  Cosmology places the strongest current limits on the sum of neutrino masses.
Future observations will further improve the sensitivity and this will require
accurate cosmological simulations to quantify possible systematic uncertainties
and to make predictions for nonlinear scales, where much information resides.
However, shot noise arising from neutrino thermal motions limits the accuracy
of simulations. In this paper, we introduce a new method for simulating
large-scale structure formation with neutrinos that accurately resolves the
neutrinos down to small scales and significantly reduces the shot noise. The
method works by tracking perturbations to the neutrino phase-space distribution
with particles and reduces shot noise in the power spectrum by a factor of
$\mathcal{O}\left(10^2\right)$ at $z=0$ for minimal neutrino masses and
significantly more at higher redshifts, without neglecting the back-reaction
caused by neutrino clustering. We prove that the method is part of a family of
optimal methods that minimize shot noise subject to a maximum deviation from
the nonlinear solution. Compared to other methods we find permille level
agreement in the matter power spectrum and percent level agreement in the
large-scale neutrino bias, but large differences in the neutrino component on
small scales. A basic version of the method can easily be implemented in
existing N-body codes and makes it possible to run neutrino simulations with
significantly reduced particle load. Further gains are possible by constructing
background models based on perturbation theory. A major advantage of this
technique is that it works well for all masses, enabling a consistent
exploration of the full neutrino parameter space.
