We show how Newtonian cosmological simulations can be employed to investigate the non-linear evolution of two particle species in a relativistic context. We discuss the application for massive neutrinos and other multi-species systems such as Cold Dark Matter (CDM) plus baryons or Warm Dark Matter (WDM). We propose a method that allows us to perform simulations including massive neutrinos and general relativistic effects at almost the same computational cost as ordinary CDM only N-body simulations, employing tailor-made initial conditions and a dictionary for the interpretation of the simulation output.

Initial conditions for cosmological N-body simulations are usually calculated by rescaling the present day linear power spectrum obtained from an Einstein-Boltzmann solver to the initial time employing the scale-independent matter growth function. For the baseline Lambda-CDM model, this has been shown to be consistent with General Relativity (GR) even in the presence of relativistic species such as photons. We show that this approach is not feasible in cosmologies with massive neutrinos and present an alternative method employing the Newtonian motion gauge framework.