  The non-trivial phase-space distribution of relic neutrinos is responsible
for the erasure of primordial density perturbations on small scales, which is
one of the main cosmological signatures of neutrino mass. In this paper, we
present a new code, FastDF, for generating 1%-accurate particle realisations of
the neutrino phase-space distribution using relativistic perturbation theory.
We use the geodesic equation to derive equations of motion for massive
particles moving in a weakly perturbed spacetime and integrate particles
accordingly. We demonstrate how to combine geodesic-based initial conditions
with the $\delta f$ method to minimise shot noise and clarify the definition of
the neutrino momentum, finding that large errors result if the wrong
parametrisation is used. Compared to standard Lagrangian methods with ad-hoc
thermal motions, FastDF achieves substantial improvements in accuracy. We
outline the approximation schemes used to speed up the code and to ensure
symplectic integration that preserves phase-space density. Finally, we discuss
implications for neutrino particles in cosmological N-body simulations. In
particular, we argue that particle methods can accurately describe the neutrino
distribution from $z=10^9$, when neutrinos are linear and ultra-relativistic,
down to $z=0$, when they are nonlinear and non-relativistic. FastDF can be used
to set up accurate initial conditions (ICs) for N-body simulations and has been
integrated into the higher-order IC code monofonIC.

Cosmology and Nongalactic Astrophysics