  We develop a higher-order perturbation theory for large-scale structure
formation involving a free-streaming hot or warm dark matter species. We focus
on the case of mixed cold dark matter and massive neutrinos, although our
approach is applicable also to a single warm dark matter species. In order to
capture the suppressed growth of neutrino density perturbations on small
scales, we account for the full momentum dependence of the phase space
distribution using the Vlaslov equation, and derive from it a formal
closed-form nonlinear equation for the neutrino density. Using a systematic
perturbative expansion of this equation we compute high-order corrections to
the neutrino density contrast without the explicit need to track the perturbed
neutrino momentum distribution. We calculate the leading-order total matter
bispectrum for several neutrino masses. Using our result as a benchmark, we
test the accuracy of the fluid approximation and a linear approximation used in
perturbative and N-body analyses, as well as a new hybrid approach that
combines the exact linear evolution with the nonlinear structure of the fluid
equations. Aiming at $\lesssim1%$ accuracy, we find that the total matter
bispectrum with a low neutrino mass m=0.046 eV can be reproduced by all but the
fluid approximation, while for larger neutrino masses m=0.46-0.93 eV only the
hybrid approach has the desired accuracy on a large range of scales. This
result serves as a cautionary note that approximate nonlinear models of
neutrino clustering that reproduce the gross features of some observables may
not suffice for precision calculations, nor are they guaranteed to apply to
other observables. All of the approximation schemes fail to reproduce the
bispectrum of the neutrino density perturbations at better than 20% accuracy
across all scales, indicating that an exact treatment of nonlinear neutrino
perturbations is necessary.
