  Many-body densities and correlation functions are of paramount importance for
understanding quantum many-body physics. Here, we present a method to compute
them; our approach is general and based on the action of bosonic or fermionic
annihilation field operators on the many-body wavefunction. We analyze $N = 6$
quasi-one-dimensional harmonically-trapped bosons with weak to strong contact
interaction strength up to the Tonks-Girardeau limit with infinite repulsion
using the MultiConfigurational Time-Dependent Hartree method for
indistinguishable particles (MCTDH-X). We compare our MCTDH-X solutions to the
analytical ones in the infinite repulsion regime as well as to the so-called
correlated pair wavefunction approach and find a good agreement. Since
numerical approximations are not bound to the cases where analytical solutions
are known, we thus demonstrate a general method to investigate high-order
reduced density matrices and correlation functions in systems for which
analytical solutions are unknown. We trace the build-up of correlation features
in the crossover from weak interactions to the Tonks-Girardeau limit and find
that the higher-order correlation functions and densities resemble those in the
Tonks-Girardeau limit for way smaller interactions than anticipated from just
the one-body density.
