  We present a novel method for including the impact of massive neutrinos in
cold dark matter N-body simulations. Our approach is compatible with widely
employed Newtonian N-body codes and relies on only three simple modifications.
First, we use commonly employed backscaling initial conditions, based on the
cold dark matter plus baryon power spectrum instead of the total matter power
spectrum. Second, the accurate Hubble rate is employed in both the backscaling
and the evolution of particles in the N-body code. Finally, we shift the final
particle positions in a post-processing step to account for the integrated
effect of neutrinos on the particles in the simulation. However, we show that
the first two modifications already capture most of the relevant neutrino
physics for a large range of observationally interesting redshifts and scales.
The output of the simulations are the cold dark matter and baryon distributions
and can be analysed using standard methods. All modifications are simple to
implement and do not generate any computational overhead. By implementing our
methods in the N-body codes GADGET and gevolution, we show that any
state-of-the-art Newtonian N-body code can be utilised out of the box. Our
method is also compatible with higher order Lagrangian perturbation theory
initial conditions and accurate for combined neutrino masses of up to at least
0.3 eV. Being formulated in relativistic gauge theory, in addition to including
the impact of massive neutrinos, our method further includes relativistic
corrections relevant on the large scales for free.
