  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%. However, we find that first-order Zel'dovich ICs lead to much
larger errors, even when starting as early as $z=127$, exceeding 1% at $z=0$
for $k>0.5\text{ Mpc}^{-1}$ for the power spectrum and $k>0.1\text{ 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%
agreement between 2LPT and 3LPT at $z=0$.
