Higher order initial conditions with massive neutrinos

The discovery that neutrinos have mass has important consequences for cosmology. The main effect of massive neutrinos is to suppress the growth of cosmic structure on small scales. Such growth can be accurately modelled using cosmological N-body simulations, but doing so requires accurate initial conditions (ICs). There is a trade-off, especially with first-order ICs, between truncation errors for late starts and discreteness and relativistic errors for early starts. Errors can be minimized by starting simulations at late times using higher order ICs. In this paper, we show that neutrino effects can be absorbed into scale-independent coefficients in higher order Lagrangian perturbation theory (LPT). This clears the way for the use of higher order ICs for massive neutrino simulations. We demonstrate that going to higher order substantially improves the accuracy of simulations. To match the sensitivity of surveys like DESI and Euclid, errors in the matter power spectrum should be well below $1{{\ \rm per\ cent}}$. However, we find that first-order Zel'dovich ICs lead to much larger errors, even when starting as early as z = 127, exceeding $1{{\ \rm per\ cent}}$ at z = 0 for k > 0.5 Mpc^-1 for the power spectrum and k > 0.1 Mpc^-1 for the equilateral bispectrum in our simulations. Ratios of power spectra with different neutrino masses are more robust than absolute statistics, but still depend on the choice of ICs. For all statistics considered, we obtain $1{{\ \rm per\ cent}}$ agreement between 2LPT and 3LPT at z = 0.