  We use a set of N-body simulations employing a modified gravity (MG) model
with Vainshtein screening to study matter and halo hierarchical clustering. As
test-case scenarios we consider two normal branch Dvali-Gabadadze-Porrati
(nDGP) gravity models with mild and strong growth rate enhancement. We study
higher-order correlation functions $\xi_n(R)$ up to $n=9$ and associated
hierarchical amplitudes $S_n(R)\equiv\xi_n(R)/\sigma(R)^{2n-2}$. We find that
the matter PDFs are strongly affected by the fifth-force on scales up to
$50h^{-1}$Mpc, and the deviations from GR are maximised at $z=0$. For reduced
cumulants $S_n$, we find that at small scales $R\leq10h^{-1}$Mpc the MG is
characterised by lower values, with the deviation growing from $7\%$ in the
reduced skewness up to even $40\%$ in $S_5$. To study the halo clustering we
use a simple abundance matching and divide haloes into thee fixed number
density samples. The halo two-point functions are weakly affected, with a
relative boost of the order of a few percent appearing only at the smallest
pair separations ($r\leq 5h^{-1}$Mpc). In contrast, we find a strong MG signal
in $S_n(R)$'s, which are enhanced compared to GR. The strong model exhibits a
$>3\sigma$ level signal at various scales for all halo samples and in all
cumulants. In this context, we find that the reduced kurtosis to be an
especially promising cosmological probe of MG. Even the mild nDGP model leaves
a $3\sigma$ imprint at small scales $R\leq3h^{-1}$Mpc, while the stronger model
deviates from a GR-signature at nearly all scales with a significance of
$>5\sigma$. Since the signal is persistent in all halo samples and over a range
of scales, we advocate that the reduced kurtosis estimated from galaxy
catalogues can potentially constitute a strong MG-model discriminatory as well
as GR self-consistency test.
