Fully relativistic treatment of light neutrinos in $N$-body simulations

Cosmological $N$-body simulations are typically purely run with particles using Newtonian equations of motion. However, such simulations can be made fully consistent with general relativity using a well-defined prescription. Here, we extend the formalism previously developed for $Lambda$CDM cosmologies with massless neutrinos to include the effects of massive, but light neutrinos. We have implemented the method in two different $N$-body codes, CONCEPT and PKDGRAV, and demonstrate that they produce consistent results. We furthermore show that we can recover all appropriate limits, including the full GR solution in linear perturbation theory at the per mille level of precision.