  We use a suite of high resolution state-of-the-art N-body Dark Matter
simulations of chameleon $f(R)$ gravity to study the higher order volume
averaged correlation functions $\bar{\xi_n}$ together with the hierarchical
$n$-th order correlation amplitudes $S_n=\bar{\xi}_n/\bar{\xi}_2^{n-1}$ and
density distribution functions (PDF). We show that under the non-linear
modifications of gravity the hierarchical scaling of the reduced cumulants is
preserved. This is however characterised by significant changes of both the
$\bar{\xi_n}$ and $S_n$'s values and their scale dependence with respect to
General Relativity gravity (GR). In addition, we measure a significant increase
of the non linear $\sigma_8$ parameter reaching $14, 5$ and 0.5% in excess of
the GR value for the three flavours of our $f(R)$ models. We further note that
the values of the reduced cumulants up to order $n=9$ are significantly
increased in $f(R)$ gravity for all our models at small scales $R\simlt
30\hmpc$. In contrast the values of the hierarchical amplitudes, $S_n$'s, are
smaller in $f(R)$ indicating that the modified gravity density distribution
functions are deviating from the GR case. Furthermore we find that the redshift
evolution of relative deviations of the $f(R)$ hierarchical correlation
amplitudes is fastest at high and moderate redshifts $1\leq z \leq4$. The
growth of these deviations significantly slows down in the low redshift
universe. We also compute the PDFs and show that for scales below $\sim
20\hmpc$ they are significantly shifted in $f(R)$ gravity towards the low
densities. Finally we discuss the implications of our theoretical predictions
for measurements of the hierarchical clustering in galaxy redshifts surveys,
including the important problems of the galaxy biasing and redshifts space
distortions.
