  We perform detailed comparison of the semi-analytic halo model predictions
with measurements in numerical simulations of the two and three point
correlation functions (3PCF), as well as power spectrum and bispectrum. We
discuss the accuracy and self-consistency of the halo model description of
gravitational clustering in the non-linear regime and constrain halo model
parameters. We exploit the recently proposed multipole expansion of three point
statistics that expresses rotation invariance in the most natural way. This not
only offers technical advantages by reducing the integrals required for the
halo model predictions, but amounts to a convenient way of compressing the
information contained in the 3PCF. We find that, with an appropriate choice of
the halo boundary and mass function cut-off, halo model predictions are in good
agreement with the bispectrum measured in numerical simulations. However, the
halo model predicts less than the observed configuration dependence of the 3PCF
on ~ Mpc scales. This effect is mainly due to quadrupole moment deficit,
possibly related to the assumption of spherical halo geometry. Our analysis
shows that using its harmonic decomposition, the full configuration dependence
of the 3PCF in the non-linear regime can be compressed into just a few numbers,
the lowest multipoles. Moreover, these multipoles are closely related to the
highest signal to noise eigenmodes of the 3PCF. Therefore this estimator may
simplify future analyses aimed at constraining cosmological and halo model
parameters from observational data.
