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%. 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$.


INTRODUCTION
The neutrino content of the Universe, Ω = /(93 eV ℎ 2 ), becomes a powerful probe for cosmology once the implied neutrino masses are confronted with data from neutrino oscillations (Esteban et al. 2020) and the kinematics of -decay (Aker et al. 2021). A non-zero detection of Ω would be consequential for fundamental physics. It would confirm that a background of relic neutrinos survived until the epoch of structure formation, provide insight into the origin of neutrino mass, and constrain the search for dark matter and dark sectors. Oscillation experiments provide a lower bound of > 0.058 eV, while cosmology provides upper bounds of < 0.15 eV or better assuming ΛCDM (Palanque- Delabrouille et al. 2020;Choudhury & Hannestad 2020;Porredon et al. 2021;Di Valentino et al. 2021), with ongoing and future surveys promising significant further improvement. Planck and future cosmic microwave background experiments, together with largescale structure surveys like DESI, Euclid, and Vera Rubin, could achieve sensitivities in the 0.01 -0.02 eV range (Hamann et al. 2012;Abazajian et al. 2015;Brinckmann et al. 2019;Chudaykin & Ivanov 2019). Such small shifts in neutrino mass correspond to tiny 0.5% -1.5% effects on the power spectrum of matter fluctuations on 0.1 Mpc −1 to 1 Mpc −1 scales, requiring theoretical predictions that are at least as accurate.
With this goal in mind, many groups have studied the effects of massive neutrinos on large-scale structure. At early times and on large enough scales, perturbation theory is the method of choice for this purpose. Cosmological perturbation theory (Bernardeau et al. 2002) is essential for providing analytical insight and a necessary complement to more expensive numerical simulations. The effects of neutrinos on the nonlinear matter power spectrum were first calculated at one-loop by Saito et al. (2008) and Wong (2008). Subsequent work has dealt more realistically with the neutrino phasespace distribution (Shoji & Komatsu 2010;Dupuy & Bernardeau 2014;Blas et al. 2014;Führer & Wong 2015;Levi & Vlah 2016;Chen et al. 2021), which parallels similar efforts on the numerical simulations side. Other advances were made by including neutrinos in the effective field theory of large-scale structure (Senatore & Zaldarriaga 2017;Colas et al. 2020) and using time renormalisation group perturbation theory (Lesgourgues et al. 2009;Upadhye 2019), which improved agreement with -body simulations. More closely related to this work, Wright et al. (2017) extended the hybrid COLA simulation method to cases with massive neutrinos using second-order Lagrangian perturbation theory (2LPT) and Aviles & Banerjee (2020) incorporated nonlinear neutrino effects in Lagrangian perturbation theory up to third order (3LPT). On the numerical simulations side, where higher-order LPT has been used to great effect to produce accurate initial conditions (ICs) for conventional simulations without massive neutrinos (Scoccimarro 1998;Sirko 2005;Crocce et al. 2006), neutrino effects have not been included and higher-order LPT is therefore rarely used for neutrino simulations (but see Brandbyge et al. 2008;Yèche et al. 2017). In this work, we propose a novel scheme for generating LPT ICs for neutrino simulations based on all-order recursive solutions in the small-scale limit. We also generate ICs based on a full calculation of scale-dependent neutrino effects in 2LPT, dealing with framelagging terms following Aviles & Banerjee (2020), and find near perfect agreement with our scheme in the final simulation product. This demonstrates that neutrino effects can be implemented beyond first order by working in the small-scale limit, paving the way for accurate neutrino simulation ICs.
-body simulations are used to solve for the nonlinear gravitational dynamics of matter on small scales, where perturbation theory fails. Cosmological simulations with ICs based on LPT were pioneered by Frenk et al. (1983); Klypin & Shandarin (1983) and Efstathiou et al. (1985). Mixed dark matter simulations with subelectronvolt neutrinos were first carried out by Brandbyge et al. (2008); Brandbyge & Hannestad (2009) ;Viel et al. (2010). We refer the reader to  for a review of neutrino simulation methods. As with perturbation theory, the accuracy of modern surveys places stringent demands on simulations, popularly expressed as a requirement for 1% accurate calculations of the matter power spectrum (Schneider et al. 2016). A major source of uncertainty concerns the interface between perturbation theory and simulation, in the form of ICs, and associated transients (Scoccimarro 1998). We may distinguish two major sources of uncertainty related to the choice of ICs (Efstathiou et al. 1985;Michaux et al. 2021). The first arise from discrepancies between the ICs and the actual nonlinear solution at the initial time. When the solution is calculated perturbatively at order , this uncertainty can be understood as the truncation error introduced by neglecting terms of order + 1 and greater. The second source of uncertainty relates to discreteness effects that build up over time as the continuous fluid equations are solved by means of a discrete particle representation (Marcos et al. 2006;Garrison et al. 2016). There is a tension between these two, as early starts minimize truncation errors but entail larger discreteness errors, while late starts do the opposite. For example, the first-order solution of Zel'Dovich (1970) has the largest possible truncation error, driving practitioners to start simulations early when higherorder corrections are small. However, such simulations manifest a greater dependence on particle resolution due to discreteness errors. While such errors can be corrected (Garrison et al. 2016), this reasoning provides strong motivation for using higher-order ICs at late times (Michaux et al. 2021).
Neutrinos complicate this picture in two ways. First, neutrinos introduce an additional length scale into the problem. Due to their large thermal velocities, neutrinos free stream out of potential wells, otherwise stated in terms of a suppression of clustering on scales smaller than a typical free-streaming length (Lesgourgues & Pastor 2006). This in turn causes a scale-and time-dependent suppression of dark matter and baryon clustering that must be accounted for in the initial conditions. Zennaro et al. (2017) showed how to incorporate such scale-dependence in a first order back-scaling procedure, but a consistent framework for higher-order ICs has thus far been lacking. We note that after we submitted our paper to the journal, Heuschling et al. (2022) presented a recipe for second-order neutrino ICs. Like us, they use a back-scaled transfer function for the cold dark matter and baryon species. The second complication is that late-time observables are more strongly correlated with the initial conditions and less determined by the internal structure of halos, when clustering is suppressed on small scales. This means that simulations with different neutrino masses are affected by errors to different degrees, contaminating ratios such as the suppression of the matter power spectrum. We will show that such ratios are more robust than absolute statistics, but still depend on the choice of initial conditions on small scales.
The paper is organized as follows. We begin by summarizing our recipe for generating higher-order ICs for neutrino simulations in section 2. The second part of the paper is concerned with a derivation of the higher-order solutions necessary for ICs, starting with the set-up in section 3, limiting solutions at all orders in section 4.1, and the full second-order solution in section 4.2. The final third of the paper contains a systematic analysis of higher-order ICs in section 5. Finally, we conclude in section 6. Throughout this paper, we use a default neutrino mass sum of = 0.3 eV to showcase our results, except where indicated otherwise.

N-BODY INITIAL CONDITIONS
We begin by outlining our approach for setting up for 3-fluid ICs with cold dark matter (c), baryons (b), and neutrinos ( ). Initially, we deal with a single cold fluid, described in terms of the the massweighted density contrast and velocity, In a final step, the cold fluid is separated into two components with distinct transfer functions. Our approach is based on a growing mode solution of the LPT equations in the small-scale limit, motivated by the hierarchy between the neutrino free-streaming scale and the nonlinear scale, fs nl , at the redshifts relevant for ICs. In section 5, we confirm that this is an excellent approximation suited for precision simulations. The recipe boils down to the following steps: (i) Compute a back-scaled transfer function cb ( ) (ii) Compute particle displacements via Eqs. (3-11) (iii) Compute particle velocities via Eqs. (12-14) (iv) Perturb particle masses and velocities via Eqs. (15)(16)(17)(18)(19) These steps can be performed using a modified version of the -IC code (Michaux et al. 2021), which we have made publicly available 1 . We briefly discuss the steps in order and then deal with possible extensions in section 2.5 and 2.6.

Transfer functions and back-scaling
In this paper, we follow the commonly used back-scaling approach. This approach begins by choosing a pivot redshift, typically = 0, where the simulation should reproduce linear theory on the largest scales. This is necessary because conventional -body codes solve Newtonian equations and therefore fail to capture the large-scale general relativistic dynamics in which matter and radiation are coupled through the Einstein-Boltzmann equations. We note that there exist alternative solutions to this problem Brandbyge et al. 2017;Fidler & Kleinjohann 2019;Tram et al. 2019;Partmann et al. 2020) as well as fully relativistic -body codes (Adamek et al. 2017;Barrera-Hinojosa & Li 2020a,b), which can avoid it altogether. In the back-scaling procedure, one uses a linear Einstein-Boltzmann code such as (Lesgourgues 2011) or (Lewis & Challinor 2011) to calculate the density transfer functions for each fluid species at pivot , which are then scaled back to the starting redshift of the simulation using the exact linear dynamics of the Newtonian code. For ΛCDM without neutrinos, this amounts to rescaling the dark matter transfer function by a constant growth factor ratio ( )/ ( pivot ).
Adding massive neutrinos makes the linear solution scaledependent, precluding a simple rescaling factor. Nevertheless, the same philosophy can be applied by solving the Newtonian dynamics of an -body code with massive neutrinos at linear order. Following Zennaro et al. (2017), we do this using a first-order Newtonian fluid approximation (Shoji & Komatsu 2010;Blas et al. 2014), but see also Heuschling et al. (2022) for a relativistic formulation. This back-scaling method for neutrino cosmologies was first implemented in the code. To streamline the procedure for the end-user and to reduce the potential for human error, we built a lightweight back-scaling library that interfaces directly with and the initial conditions generator IC. The final result of these steps is a rescaled density transfer function cb ( ) = cb ( , )/ cb ( , pivot ) · cb ( , pivot ) for a cold dark matter-baryon fluid (cb), where the growth factor ratio is computed with and the transfer function with .

Displacements
The displacement field, = − , relates the particle position to the corresponding Lagrangian coordinate . To determine , we first obtain the linear potential by solving ∇ 2 (1) ( ) = cb ( ).
Unless indicated otherwise, ∇ = ∇ . We observe that (1) is not the gravitational potential, which also includes a neutrino contribution, but a notation that reflects the fact that we are solving for the displacements of cb fluid particles. Our fast approximate 3LPT (Buchert 1994;Bouchet et al. 1995;Melott et al. 1995) scheme for the displacement field in the presence of massive neutrinos has the simple form where are scale-independent factors that capture the absence of neutrino perturbations in the small-scale limit, = / , and ( ) have the same form in terms of (1) as in ΛCDM. In the notation of Michaux et al. (2021), these are given by with higher-order potentials given by , − where commas denote partial derivatives and we sum over repeated indices. In section 4.1, we show that can be expressed in terms of the neutrino fraction, = Ω /Ω m . The correction, as it turns out, is small and approximately linear in : 1 + 2 5(2 + 3) . (11) For a minimal neutrino mass sum of = 0.06 eV, one finds 2 − 1 = 5 × 10 −4 . For our fiducial mass sum of = 0.3 eV, it is 0.3%. At = 1 eV, the effect is about one per cent. The third-order correction 3 is larger, but since (3) is suppressed by another power of the growth factor, the overall impact is smaller.

Velocities
The velocity field is cb = d /d . Given a satisfactory scheme for computing the displacement field, the time derivative can be evaluated numerically. This is our preferred method, since it requires no additional approximations. However, a faster method that avoids calculating higher order terms more than once is to use the asymptotic logarithmic growth rate to convert displacements to velocities, setting cb = ∞ (1) + 2 2 (2) (13) By construction, this gives the correct particle velocities on small scales. To recover also the correct behaviour on horizon scales, we add a large-scale correction ( ) cb given by where cb is the dimensionless energy flux transfer function computed with . We verified that the resulting simulated power spectrum agrees with linear theory to better than 0.1% at the pivot redshift of = 0 on large scales. However, this approximation neglects possible nonlinear effects of scale-dependent growth on particle velocities. Another alternative is to rescale the velocities by the scale-dependent growth rate (Zennaro et al. 2017), which faces a similar problem beyond linear order.

Additional steps for 3-fluid ICs
The steps above are sufficient for simulations with neutrinos and a single cold fluid. To separate this cold fluid into baryon and CDM components with distinct transfer functions, we follow the approach of Hahn et al. (2021). In short, the component densities are related to the mass-weighted average via 2 where the difference variable, bc = b − c , is constant at first order. The velocity difference too is conserved and vanishes at all orders: bc = b − c = 0. These results, derived for ΛCDM without massive neutrinos , carry over to the neutrino case, essentially due to the fact that the neutrino contribution cancels in the difference equations (Appendix A). The transfer function 2 We remind the reader that = Ω /Ω cb for ∈ {c, b} even as = Ω /Ω m = Ω /(Ω cb +Ω ) = 1− cb . Furthermore, bc ≠ cb and bc ≠ cb . difference, bc ( ) = b ( ) − c ( ), is computed with at the pivot redshift and, since it is conserved, is not scaled back.
After assigning displacements and velocities to both particle species using the mass-weighted average fields, the density difference is implemented by setting the masses to with¯the mean particle mass for type ∈ {c, b}. Perturbing the masses, rather than the displacements, was found by Hahn et al. (2021) to limit discreteness errors. By construction, Newtonian simulations with initial conditions set up using the above procedure, reproduce the expected evolution of two cold fluids with a shared velocity field and a relative density contrast that is approximately conserved. However, like the largescale velocity correction (14), a further modification is needed to bring the dynamics back into agreement with at first order: where ∞ ( ) is the small-scale growth factor at the starting redshift and of the dimensionless energy flux transfer functions is computed with at the pivot redshift.

Neutrino particles
Massive neutrinos can be included in -body codes using a variety of methods. The most common approach is to solve for the neutrino perturbations self-consistently by including them as a separatebody particle species (Brandbyge et al. 2008;Viel et al. 2010). Initial conditions are then also needed for these neutrino particles.
Capturing the full neutrino phase-space distribution is non-trivial even in linear theory and it is therefore not sufficient to compute only the first two moments, as is done for baryons and CDM. Accurate neutrino particle initial conditions can be generated by integrating geodesics from high redshift (Ma & Bertschinger 1994;Adamek et al. 2017), where the perturbed phase-space distribution can be expressed analytically (Ma & Bertschinger 1995), but care must be taken that the equations of motion remain valid in the ultrarelativistic régime (Elbers, in prep.). This procedure can be carried out efficiently using our F DF code. We stress that the focus of this paper is on dark matter and baryon ICs and the results apply regardless of whether the neutrino implementation uses particles.

Scale-dependent effects
Finally, we verified the approximations above by performing a full calculation of scale-dependent effects on the second-order displacement field. This is done by replacing (7) with a convolution of two copies of the first-order potential (1) ( ), modulated by kernels (2) ( 1 , 2 ) and (2) ( 1 , 2 ), computed in section 4.2. This numerical calculation is expensive, but we will show in section 5 that simulations with ICs based on the full calculation agree extremely well with those based on the approximate scheme described above.
The reason for this is the hierarchy of scales, fs nl , which implies that higher-order corrections are important only on scales where neutrinos do not cluster, at least at redshifts that are relevant for ICs. Since the overall impact of the third-order correction factor,

THEORETICAL SET-UP
We now proceed with the set-up of a 3-fluid model, which is solved in section 4. We consider three fluids indexed by ∈ {c, b, } for cold dark matter, baryons, and neutrinos. Throughout, we will treat baryons like dark matter particles and denote the mass-weighted CDM-baryon fluid by subscript cb. Let ( ) be the density, ( ) the peculiar velocity flow, and ( ) the stress tensor. We also write = /¯− 1 for the density contrast.

Euler equations
Taking moments of the Boltzmann equation yields the Euler fluid equations (Bernardeau et al. 2002) where is conformal time, = / 2 is the Hubble constant (given explicitly below) and Φ the Newtonian potential. While the neutrino distribution function and its higher-order moments are complicated, the stress tensor can be neglected for the cold dark matter and baryon fluids on the scales of interest, c = b = 0. Taking the mass-weighted average of the cold dark matter and baryon equations, we obtain at all orders (see Appendix A) The potential is given by Poisson's equation, in terms of the total matter density, m = cb cb + , which includes a massive neutrino contribution. To complete the system, we assume the linear response approximation for the neutrino density: where lin ( ) refers to the density transfer function of ∈ { , cb} computed in relativistic linear perturbation theory with . The total matter density contrast is then where we have introduced the convenient notation = lin /( cb lin cb ) for the linear theory ratio. The linear response approximation is accurate while neutrinos and dark matter remain in phase, which is a reasonable assumption at the early times considered here (see below). Inserting this in (24) yields is written in terms of present-day values. We look for a growing solution of the form cb ( , ) = cb ( , ) cb ( , 0 ). Linearising (22-24), we find 2 cb + cb = 0 (1 + ) cb .
In contrast to the ΛCDM case, this equation is scale-dependent due to the appearance of ( ). To proceed, we will take the limit → ∞. Since lim →∞ ( ) = 0, we simply obtain We denote the solution of (29) by ∞ to indicate that this is the small-scale solution. At this point, an equally valid description could be given in the large-scale limit or indeed for an arbitrary pivot scale. We deliberately choose the small-scale limit for two reasons. First, most simulations are not large enough to realize the largescale limit. Second, we are interested in nonlinear corrections to the initial conditions which are negligible on large scales.

Asymptotic form
We can find an analytic 3 solution to (29) if the contribution of radiation to the Hubble rate is neglected. We will return to this point further below. For now, let us assume that In this case, the growing mode can be expressed in terms of the hypergeometric function as (see Appendix B) with Λ = Ω Λ /Ω m and = √︁ 1 + 24(1 − )/4 − 1/4. This is normalized such that lim →0 ∞ / = 1. Taking = 0, we recover the ΛCDM solution with = 1 (Rampf et al. 2015). Taking instead Λ → 0, we recover the solution during matter domination (MD) which agrees with Bond et al. (1980). For ΛCDM without massive neutrinos, accurate nonlinear predictions can be made by substituting the growth factor for the scale factor, → , in solutions obtained for the Einstein-de Sitter model. This is facilitated by using the growth factor as time variable (e.g. Matsubara 2015; Rampf et al. 2015Rampf et al. , 2021. Here, we will pursue a similar strategy and make a change of time variables to ∞ . Defining the quantity and the new velocity variable cb = ∞ , the fluid equations can be rewritten as where the rescaled potential = Φ/( 0 ∞ ) is given in terms of a convolution, denoted by * , of cb and the linear response (1 + ).
Although written in terms of ∞ , this is completely general.  Given suitable boundary conditions, Eqs. (34-36) are analytic at ∞ = 0. In particular, we require that ini m = ini cb = 0. This agrees with our use of growing mode solutions for particle displacements, ↦ → + , where the unperturbed particle grid represents a uniform density field. The scaling, 2 ∝ −3 , of the Hubble rate at early times ensures that such mass transport problems are well-posed (Brenier et al. 2003;Rampf et al. 2015). This scaling does not hold in the presence of radiation, a problem that already occurs in ΛCDM on account of the cosmic microwave background radiation, but is certainly made worse by the inclusion of massive neutrinos, which scale like radiation in the relativistic régime. Therefore, we need to start the integration at a time when the relativistic contribution of neutrinos to the Hubble expansion can be neglected. Note that we make this assumption to ensure a consistent mathematical framework for the higher-order LPT solutions. However, it is not needed for the linear transfer functions, the back-scaling procedure or in the -body code itself. In each of those cases, we do take the relativistic neutrino contribution into account.
Before proceeding, let us give the following convenient expression for ∞ in the limit Λ → 0: Both numerator and denominator scale approximately as (1− ) 1/2 . The numerator is simply the exponent of the growing mode in (32), while the dependence of the denominator can be traced to the appearance of 0 on the right-hand side of (29). The resulting smallness of ∞ − 1 explains why neutrino corrections at th order are small relative to ∞ : the lack of neutrino clustering is largely compensated by slower growth of the linear solution. In the next section, we will validate the assumptions made up to this point.

Validity of assumptions
Central to the approach of section 4 is the linear response approximation (25) for the nonlinear neutrino density, ( ). This approximation is very accurate at early times, but underestimates neutrino clustering on small scales and neglects the phase shift between neutrinos and dark matter that builds up at late times (see Fig. 6 in Elbers et al. 2021). The top panel of Fig. 1 shows the nonlinear neutrino density contrast, computed from a simulation with neutrino particles, relative to the linear neutrino response evaluated at = 0.60 Mpc −1 . The neutrino mass is = 0.3 eV. The figure suggests that the approximation is valid at this scale up to ≈ 1.5, when perturbation theory has presumably already broken down. Hence, approximation (25) is well-suited for our application at much higher redshifts.
A second approximation is that we neglect the contribution of the relativistic tail of the neutrino distribution to the Hubble rate in (30). We reiterate that this approximation is only made for the calculation of the higher-order kernels and not in any of the calculations at first order. The middle panel of Fig. 1 shows that this approximation is accurate to better than 1% for > 0.01, for our default neutrino mass of = 0.3 eV. In particular, at the fiducial starting redshift of = 31, the error is 0.3%. We are helped in this regard by our preference for late starts. Finally, we assume that ∞ is constant in section 4.1. The bottom panel of Fig. 1 shows that this is an excellent approximation, except at late times during Λ-domination. The figure suggests that there is a window where all assumptions are valid, potentially allowing us to push to even later starts, with the breakdown of LPT likely being the limiting factor.

LAGRANGIAN APPROACH
In the Lagrangian approach to gravitational instability (Zel'Dovich 1970;Buchert 1989;Moutarde et al. 1991;Bouchet et al. 1992;Gramann 1993;Buchert & Ehlers 1993;Bouchet et al. 1995), the objective is to describe fluid particle trajectories in terms of a displacement field . We use the Helmholtz decomposition, writing the Laplacian of a smooth vector field as What remains is to solve for the longitudinal and transverse derivatives. The displacement is related to the Eulerian density, cb , through the mass conservation equation where ( ) is the determinant of the Jacobian of the coordinate transformation, = / , given by Let ( / ∞ ) L = ∞ + cb · ∇ be the Lagrangian derivative. The Lagrangian form of the Euler equation (34) can be written as where we used cb = ( / ∞ ) L and introduced the linear operator For the large majority of configurations, the system attains the approximate value. The shaded region indicates the range of scales for which the power spectrum of · (2) is at least 0.01% of that of · (1) .
Using (36) and taking the divergence and curl of (42), we find that the evolution of the displacement is governed by To facilitate a fully Lagrangian description, we define the framelagging terms ( Frame-lagging terms arise from mapping the Eulerian neutrino response to Lagrangian coordinates. We give explicit expressions up to second order in Appendix C. Transforming the derivatives on the left-hand side of (44) and (45) using = ( / ) = −1 and using the Monge-Ampère equation (40), we write these equations in Lagrangian coordinates as It will be the task in the following sections to find perturbative solutions for . We perform an expansion in displacements, writing where ( ) is of order (1) .

Limiting solutions
Having set up the Lagrangian equations for the neutrino-cb fluid model, we are now in a position to look for approximate solutions.
The aim is to find expressions for the displacement on large and small scales. In the small-scale limit, neutrinos do not cluster and only contribute to the background expansion as encoded by ∞ .
Meanwhile, in the large-scale limit, neutrinos cluster like cold dark matter and one recovers behaviour analogous to ΛCDM. In both cases, we can find simple solutions in the form of LPT recursion relations (Rampf 2012;Zheligovsky & Frisch 2014;Matsubara 2015;Rampf et al. 2015;Schmidt 2021). These limiting solutions will be used as initial conditions for the numerical integration of the general problem and provide the basis for the recipe of section 2.
In this section, we assume that ∞ = constant, which is exact during matter domination (Eq. (37)), and a very good approximation in general (Fig. 1). On large scales, we also have 1 + ( ) = 1 + / cb 4 and on small scales 1 + ( ) = 1. Hence, if all modes involved in the problem are either large or small, we can approximate the convolution with the neutrino response as multiplication by a constant = 1 + ( ). In such cases, the frame-lagging terms also vanish, as will be confirmed in section 4.2. Given these assumptions, (47) reduces to Using the identities −1 = (1/2) and = (1/6) , we rewrite (50) and (48) as Hence, using = + , and substituting the expansion (49), we obtain equations for the longitudinal and transverse parts at order in terms of perturbations of orders 1 + 2 = (for ≥ 2) and 1 + 2 + 3 = (for ≥ 3): , , The first-order equations separate. The longitudinal equation (53) has the particular time-dependent solution (1) = ∞ with = 1 4 √︁ 4 + 3 ∞ (8 + 3 ∞ − 4) − 3 4 ∞ + 1 2 , while the transverse equation (54) has constant and decaying solutions. Identifying the fastest growing solutions order by order, we find that ( ) ∝ ∞ . In particular, we find that the fastest growing solution at second order grows as Reinserting = 1 + ( ), we obtain a useful approximation of the magnitude of neutrino effects on the second-order coefficient, relative to the ΛCDM value of 3/7. This is shown by the dashed line in Fig. 2 for a model with = 0.3 eV at = 31. We stress that this approximation neglects the non-trivial coupling with the neutrino 4 This is not strictly true, since > cb on the largest scales due to the relativistic tail of the neutrino distribution. We ignore this small effect in the current section and in Fig. 2, but take it into account in section 4.2.
response in the general case. As we will see in the next section, the second-order solution can be described in full by two kernels, (2) ( 1 , 2 ) and (2) ( 1 , 2 ). For most configurations on the 6D Fourier space lattice that we use to generate -body ICs, both 1 and 2 are large and the result is close to the estimate of Eq. (55). However, for cases with one mode large and one mode small or for squeezed configurations with = || 1 + 2 || 1 ≈ 2 , the value may depart from this estimate, as shown by the histogram in Fig. 2. Nevertheless, the figure demonstrates that the large-and smallscale limits provide reasonable bounds on the effect at intermediate scales. Overall, the magnitude of the effect is O 10 −3 , in line with the estimate given in section 2.2 for this mass. The figure also demonstrates that the ΛCDM value of 3/7 is only reached for < 10 −3 Mpc −1 , while the second-order potential is important for > 10 −1 Mpc −1 , reflecting the hierarchy between the neutrino freestreaming scale and the nonlinear scale, fs nl , that motivates the approach of section 2.
Using ( ) ∝ ∞ , we derive recursion relations for the fastest growing solution at order ≥ 2: , , For the purposes of higher-order ICs, we are primarily interested in deriving corrections to the ΛCDM coefficients in the small-scale limit with = = 1. Reading off coefficients from (56), we find that these can be conveniently expressed in terms of Proceeding as in Appendix D, we obtain the 3LPT form given in section 2.2. Combining Eqs. (58) and (37) yields an accurate approximation of in terms of : with = √︁ 1 + 24(1 − ). For = 2, the above expression agrees with that given by Wright et al. (2017). The next section is dedicated to relaxing the assumptions on ∞ and ( ), finding the general solution at second order.

General solution
For the general solution, we need to deal with the frame-lagging terms ( ). Here, we will follow the approach of Aviles & Banerjee (2020). We are interested in solutions at second order. The transverse equation (48) only has non-trivial solutions for ≥ 3. Therefore, we concentrate on the longitudinal part. We repeat (47) for convenience: Using (41) and −1 = ∞ =0 [( − ) ] , we can write this up to second order in the displacement: where the second-order frame-lagging terms, (2) , are given in Appendix C. At first order, the displacement admits a growing solution (1) ∝ (1) with a growth factor that satisfies This is simply a reformulation of the Eulerian equation for the first-order growth factor (28). Using the expansion (49) in (61) and collecting second-order terms then yields , , In Fourier space, each of the quadratic terms in (63), including the second-order frame-lagging term, is a convolution of derivatives of (1) ( 1 ) and (1) ( 2 ). Expressing the displacements in terms of potentials as and identifying terms, we thus obtain where ∫ 1 , 2 = ∫ d 1 d 2 (2 ) −6 (3) ( 1 + 2 − ) and 12 = 1 · 2 and = (1) ( ) for = 1, 2. Notice the similarity of this equation with Eq. (7). The difference is that the two terms now have distinct scale-and time-dependent coefficients satisfying where the functions and are given by for = || 1 + 2 ||. The terms in square brackets correspond to the frame-lagging terms. In the small-scale limit with , 1 , 2 fs , we have = = 0. Hence, (2) = ( 2) and (65) factorizes as in Eq. (7). Similarly, in the large-scale limit with , 1 , 2 fs , we obtain again the approximate form described in section 4.1 with = ≈ / cb . In both limits, the frame-lagging terms drop out, as anticipated. Intermediate configurations will deviate from the asymptotic solutions, as was already discussed in section 4.1 and shown in Fig. 2.
For the numerical solution, we begin the integration at a time when the non-relativistic neutrino fraction is 50%. For the fiducial neutrino mass, (2) ( 1 , 2 ) and (2) ( 1 , 2 ). In the asymptotic approximation (black), we use Eqs. (4) and (13), but truncate third-order terms. In the ΛCDM approximation (red), we additionally set 2 = 1. The vertical dotted line is the particle Nyquist frequency the second-order kernels, using the approximate model of Eq. (55) as initial conditions. The results, projected onto the -axis, are shown in Fig. 2. When generating 2LPT particle initial conditions, we begin by generating a realisation of the back-scaled first-order potential, (1) . We then perform the convolution integral of Eq. (65) explicitly, interpolating from tables of (2) , ( , 1 , 2 ). To ensure completion in a reasonable time frame, we impose cut-offs at 1 ≤ cut and 2 ≤ cut . We performed convergence tests to ensure that the results are independent of the cut-off scale, finding that a cutoff at cut = 1 Mpc −1 was more than adequate for the resolutions considered in this paper. Table 2. Description of the simulations. The listed particle mass, , refers to the cb particles. The neutrino fraction is listed as = Ω /(Ω cb + Ω ).
All simulations used the same random phases in an = 800 Mpc cube.

RESULTS
We will now discuss the power spectra, bispectra, and halo mass functions of massive neutrino simulations with different ICs. We introduce our simulation suite in section 5.1. We then consider the impact of different approximation schemes for the second-order kernels in section 5.2 and follow it up with a comparison of Zel'dovich (ZA), 2LPT, and 3LPT ICs at various starting redshifts in section 5.3. Finally, we consider the impact of ICs on the suppression of the power spectrum as a function of neutrino mass in section 5.4.

Simulations
We use the cosmological hydrodynamics code (Schaller et al. 2016(Schaller et al. , 2018, which uses task-based parallelism, asynchronous communication, fast neighbour finding, and vectorised operations to achieve significant speed-ups. The code uses the Fast Multipole Method (FMM) for short-range gravitational forces and the Particle Mesh method for long-range forces. Neutrinos are modelled as a separate particle species. We employ the method to suppress the effects of shot noise (Elbers et al. 2021) and generate neutrino particle initial conditions by integrating geodesics from high redshift using our F DF code. Additionally, we use fixed initial conditions to limit cosmic variance (Angulo & Pontzen 2016). Apart from the neutrino mass, we use cosmological parameters based primarily on Year 3 results from the Dark Energy Survey (Porredon et al. 2021) and Planck 2018 (Aghanim et al. 2020). Our choice of parameters is (ℎ, Ω m , Ω b , , ) = (0.681, 0.306, 0.0486, 2.09937 × 10 −9 , 0.967), with different choices for the neutrino density Ω . The parameters used by the gravity solver are listed in table 1 and an overview of the simulations is given in table 2.
There is a subtle point regarding comparisons between simulations with and without massive neutrinos. Codes like employ a multipole acceptance criterion to determine when the multipole approximation is sufficiently accurate to be used without further refinement. The adaptive criterion used for the runs in this paper is based on error analysis of forces on test particles. This means that the accuracy of the -body calculation depends on the number of particles contained in any given volume. When comparing two runs with equal numbers of dark matter particles, one with neutrinos and the other without, all other things being equal, forces will be calculated more accurately in the run with neutrinos. To account for this difference, we included an equal number of massless 'spectator' neutrino particles in the = 0 runs, with velocities corresponding to = 0.05 eV neutrinos. These particles contribute no forces and only affect the -body simulation through the multipole acceptance criterion, ensuring that the accuracy of the massless runs is comparable to that of the massive neutrino runs. Such massless runs are considered in section 5.4.

Validation of approximate treatment
To validate our approach, we compare three different implementations of 2LPT, based on the following models: (i) The asymptotic model of section 2 (ii) A model with ΛCDM coefficients (iii) A reference model with scale-dependent effects The first order displacements and velocities are identical in each of the approaches, obtained from the back-scaled linear power spectrum at = 0. In the asymptotic scheme, we use Eqs. (4) and (13), but truncate the 3LPT terms. In the ΛCDM approximation, we additionally set 2 = 1, which corresponds to neglecting neutrino effects at second order. Finally, we compare these two approximate methods with a reference run that relied on a numerical calculation of the scale-dependent 2LPT kernels, (2) ( 1 , 2 ) and (2) ( 1 , 2 ).
With respect to Fig. 2, the asymptotic approximation corresponds to using the small-scale limit, the ΛCDM approximation corresponds to the large-scale limit, and the reference run corresponds to the underlying histogram. We use simulations with side length = 800 Mpc and cb = 1200 3 particles. Fig. 3 shows the impact of these approximations on the power spectrum of the evolved CDM & baryon density field. The differences are most evident at = 3 (bottom panel). On the largest scales, < 0.05 Mpc −1 , nonlinear corrections are small and all simulations agree to machine precision. For > 0.05 Mpc −1 , the ΛCDM simulation systematically underestimates clustering with a maximum error of 0.04% at = 4 Mpc −1 . For the asymptotic run, the error is two orders of magnitude smaller over the same scales. Between = 31 and = 3, the evolution is virtually identical in the asymptotic and reference runs, but we begin to see some noise in the ratio on the smallest scales at = 1 (middle panel). These perturbations continue to grow until = 0 (top panel), where we find a scatter of 2×10 −4 for > 1 Mpc −1 in both the asymptotic/reference and ΛCDM/reference ratios. It is hard to attribute this noise to any particular run as the power spectrum on these scales is increasingly determined by the internal structure of poorly resolved halos. On larger scales, < 1 Mpc −1 , the asymptotic run performs extremely well with errors below 10 −5 , while the systematic deficit in the ΛCDM run persists.
These results demonstrate that, at second order, the effect of the suppressed neutrino perturbations can be absorbed into a scaleindependent factor 2 and that further scale-dependent neutrino effects are negligible as far as initial conditions are concerned. We expect that this continues to hold for third-order corrections, which are confined to even smaller scales. Including the correction factor 2 is clearly superior to simply using the ΛCDM coefficient. However, we also observe that this higher-order neutrino effect is below 0.1%, and therefore beyond the sensitivity of current experiments. Hence, we conclude that for most purposes both the ΛCDM approximation and the asymptotic approximation are justified.

Choice of LPT order and starting time
We are now in a position to study the effects of LPT order and starting time on massive neutrino simulations, using the asymptotic approximation. Fig. 4 shows the late-time power spectrum for simulations with = 800 Mpc and cb = 1200 3 particles, comparing in the first instance Zel'dovich (solid red) and 2LPT (solid black) with 3LPT (dotted gray) as a baseline. All three runs were started at = 31. The most striking observation is that the differences are much larger than those shown in Fig. 3. This means that using higherorder LPT in some fashion is more important than getting the details right. Next, we find per cent agreement between 2LPT and 3LPT over the entire range of scales probed for ≤ 1 and approximately a 1% error at = 3 for > 2 Mpc −1 . We also find that the Zel'dovich approximation performs very poorly with errors of (4, 7, 15)% for > 1 Mpc −1 at = (0, 1, 3). This well-known fact (Crocce et al. 2006) has motivated practitioners to start Zel'dovich simulations at higher redshifts, when truncation errors are smaller. We demonstrate this with Zel'dovich runs started at = 63 (dashed, red) and = 127 (dotted, red). While the agreement with the higherorder runs improves, we still find per cent agreement only up to = 0.4 Mpc −1 . Moreover, starting earlier introduces inaccuracies of a different sort. To see this, we repeat the exercise at a lower resolution with cb = 600 3 particles. The resulting power spectra at = 0 are shown in Fig. 5, with Zel'dovich runs compared against 3LPT in the top panel. We observe that for runs started at = 31 (red), the error is almost independent of resolution. However, for earlier starts at = 63 (black) and = 127 (blue), the lower resolution runs increasingly underestimate the power spectrum on small scales. This shows that while truncation errors decrease, resolution effects increase as simulations are started earlier. The pattern reverses for 2LPT (bottom panel), with earlier starts performing worse than later starts. This can easily be explained by the fact that truncation errors are much smaller for 2LPT, such that the effect of increasing discreteness errors dominates. We confirm the finding of Michaux et al. (2021) that the size of discreteness errors is independent of LPT order. This demonstrates that, at fixed resolution and LPT order, starting earlier does not guarantee convergence onto the higher-order solution. As was the case for truncation errors, discreteness errors are much larger at = 1, 3 (not shown).
We also consider three-point statistics, which are sensitive to transients from initial conditions (Crocce et al. 2006) and an interesting probe of neutrino physics (Chiang et al. 2018;Ruggeri et al. 2018;Hahn et al. 2020). For the equilateral bispectrum, ( ) = ( 1 , 2 , 3 ) with = 1 = 2 = 3 , shown in Fig. 6 at late times, the same pattern is broadly repeated as for the power spectrum. However, errors are approximately twice as large as for the power spectrum. In detail, we again find per cent agreement between 2LPT and 3LPT for ≤ 1 with larger errors on small scales at = 3. For the Zel'dovich runs, we find significant errors compared to 3LPT, even when starting at = 127, with per cent agreement only up to = 0.1 Mpc −1 at = 0, and not even there for ≥ 1. Finally, we compare the halo mass function at = 0. Halos are identified with VELOCI (Elahi et al. 2019) using a 6D friends-of-friends algorithm applied to the cb particles. Spherical overdensity masses are computed within spheres for which the density equals 200 times the mean CDM & baryon density¯c b . The reason for using¯c b instead of the total mass density¯m is that it is this cold density field that produces universal and unbiased results in halo model calculations (Ichiki & Takada 2012;Castorina et al. 2014;Massara et al. 2014). The results are shown Fig. 7. We once again find per cent agreement between 2LPT and 3LPT over the entire mass range, but large errors for the Zel'dovich runs. There is an interesting pattern in the Zel'dovich error as the starting time is varied. For late starts (solid red), the simulation agrees well at the low-mass end but underestimates the number of very massive, 10 15 , halos by more than 7%. This can be understood in terms of the deficit of power seen also in Fig. 4, resulting in a suppressed growth of large structures. Meanwhile, for early starts (dotted and dashed red), the agreement at the high-mass end improves like the small-scale power spectrum. However, the number of low-mass halos decreases by a similar factor, likely due to discreteness errors. This seems to be broadly consistent with the ΛCDM results of Michaux et al. (2021), but not with Nishimichi et al. (2019) who find little dependence on starting time at = 0.

Dependence on neutrino mass
Thus far, we have focused on a single neutrino mass of = 0.3 eV. However, it is of great interest to determine the effect of initial conditions on the suppression of the power spectrum for different neutrino masses. We consider three cases:  (iii) degenerate = 0.30 eV neutrinos ( = 0.023).
In each case, we adjust Ω cdm to keep the total matter density Ω m fixed. We primarily use lower resolution simulations with cb = 600 3 particles in an = 800 Mpc cube. First, we consider the effect of LPT order. In Fig. 8, we show the suppression of the CDM & baryon power spectrum relative to the massless case, comparing ZA/ZA (solid), 2LPT/2LPT (dashed), and 3LPT/3LPT (shaded). Evidently, it is crucial to compare like with like simulation, keeping the LPT order and starting redshift the same. Not doing so introduces large errors in the ratio, as might be expected from the fixed neutrino mass results discussed above. We illustrate this by including a dotted line for the ZA/2LPT ratio, which is clearly off the mark. However, even when comparing like with like, we find a residual error that is proportional to the neutrino mass, rises with , and peaks around the turn-over of the spoon. This feature is most clearly visible at = 1 for ZA, with a maximum error of 0.05 . The effect is already present in the initial conditions and can be explained by a mass-dependent suppression of nonlinear terms. As virialized structures grow, both the turn-over of the spoon and the peak of the error move to larger scales. At = 0, the error is 0.025 around = 0.3 Mpc −1 . On smaller scales, we see a scatter of order 0.5%, treading outside the scale-dependent error bars that correspond to a ±0.005 eV shift in . For 2LPT, both the systematic effect and the noise are greatly suppressed, resulting in 0.1%-level agreement with 3LPT even at early times.
Next, we consider the effect of the starting time of the simulation. In Fig. 9, we show the suppression for simulations with 2LPT ICs started at = 127 (solid), = 63 (dashed), = 31 (shaded). Once again, we compare like with like simulations. Even so, we find a small residual effect with earlier starts overestimating the suppression. The differences between = 31 and = 63 are minimal for both neutrino masses. However, starting at = 127 results in (0.1, 0.2) errors at = (0, 1) for > 1 Mpc −1 . These errors once again exceed the threshold for a ±0.005 eV shift in . Based on the discussion above, and given that we are using 2LPT, we expect that truncation errors are small at both redshifts. This suggests that the differences are caused by resolution effects, which grow in importance with the starting redshift. To test this, we repeated some of the simulations at a higher resolution with cb = 1200 3 particles, starting at = 127 and = 31. The ratio is shown as a dotted line in the bottom panels of Fig. 9. The agreement between the early and late starts improves to 0.1% up to = 10 Mpc −1 at = 0, comparable to the low-resolution = 63 start. However, the suppression is still slightly overestimated at = 1. One possible alternative explanation is that errors could be introduced by the back-scaling procedure (section 2.1). To test this hypothesis, we repeated some of the simulations with "forward" ICs, as in Elbers et al. (2021). We found nearly identical results for these runs, ruling out this explanation. Another possibility is that the errors could be the result of shot noise, since we use a particle-based implementation of neutrino perturbations. However, this is unlikely as the differences already appear at high redshift when shot noise is highly suppressed due to our use of the method. Finally, one might expect differences due to relativistic effects that are increasingly important for earlier starts. Once again, this is unlikely since relativistic effects would appear on the largest scales, where the differences shown in Fig. 9 are minimal. Since the error decreases for the higher resolution runs, discreteness effects likely account for the majority of the difference, with massive neutrino simulations being more sensitive to such errors, due to the suppressed growth of structure. Late starts can be utilized to minimize the effect of particle resolution, as shown in Fig. 5.

DISCUSSION
We have investigated the use of higher-order Lagrangian initial conditions (ICs) for cosmological simulations with massive neutrinos. We solved the fluid equations for a neutrino-CDM-baryon model with approximate time-dependence in the large-and small-scale limits, finding that higher-order neutrino effects can be described by scale-independent coefficients that are easy to implement in existing IC codes. To validate our approach, we constructed ICs based on a rigorous treatment of the scale-dependent neutrino response in 2LPT, obtaining agreement with our scheme to better than one part in 10 5 up to = 1 Mpc −1 in the power spectrum of the evolved CDM and baryon perturbations at late times.
Compared to these small differences, we find that the truncation error associated with using the first-order Zel'dovich approximation is much larger. For our fiducial model with = 0.3 eV and a starting redshift of = 31, the error is 4% in the power spectrum and 7% in the equilateral bispectrum around = 0.5 Mpc −1 at = 0. Ratios of statistics from simulations with different neutrino masses can be calculated much more robustly, provided that the LPT order and starting redshift are the same. Nevertheless, even such ratios have a residual dependence on the ICs. For instance, Zel'dovich ICs introduce a mass-dependent error in the suppression of the power spectrum that grows with wavenumber and redshift , peaking around the turn-over of the spoon. We also find that the starting time of the simulation has an impact on the suppression over a wide range of scales and redshifts. Simulations started at = 127 overestimate the suppression of the power spectrum on small scales, compared to later starts. While simulations can be started at higher redshifts to reduce truncation errors, this also increases the importance of particle resolution and relativistic effects. To minimize errors from initial conditions and particle resolution, simulations can be started at late times using higher-order ICs. A major target of cosmological surveys is to measure the sum of neutrino masses. Assuming the minimum value allowed under the normal mass ordering, = 0.06 eV, cosmology could provide a 3 detection and rule out the inverted mass ordering at 2 by reaching a sensitivity of 0.02 eV, which is in reach of future cosmic microwave background and large-scale structure experiments (Hamann et al. 2012;Abazajian et al. 2015;Brinckmann et al. 2019;Chudaykin & Ivanov 2019). This corresponds to detecting 1% effects on the matter power spectrum on 0.1 Mpc −1 < < 1 Mpc −1 scales. We should therefore aim for neutrino simulations with errors that are well below 1% on these scales. While Zel'dovich ICs fall short of this mark, our findings suggest that 2LPT is sufficiently accurate for most applications. Higher-order statistics at high redshift seem to be the notable exception, which could be relevant for Lyman-forest simulations. The accuracy of neutrino simulations depends on many factors: the accuracy of the linear transfer functions and back-scaling procedure (Lesgourgues & Tram 2011;Zennaro et al. 2017), the implementation of neutrino perturbations (e.g. Bird et al. 2018;Elbers et al. 2021), neutrino initial conditions (Elbers, in prep.), and dark matter and baryon initial conditions (this paper). It has now been demonstrated that each of these factors can be controlled to within 1%. The remaining uncertainty is likely dominated by the choice of gravity solver. Achieving 1% agreement between different -body codes is non-trivial even in the absence of neutrinos (Schneider et al. 2016;Garrison et al. 2019;Grove et al. 2021). Fortunately, the accuracy of -body codes should not in the first place be expected to deteriorate in the presence of neutrinos. In fact, the accuracy could even improve for particle-based implementations due to 'spectator' effects (section 5.1). A systematic comparison of neutrino simulations with different codes and identical initial conditions could establish whether this is indeed the case. Such explorations would improve our ability to simulate nonlinear clustering in Universes with massive neutrinos, allowing us to meet the demands of the next generation of surveys. = 0.30 eV (black). The bottom panels show the suppression relative to runs with the same resolution, but started at = 31. The shaded areas represent a ±0.005 eV shift (light) and a constant 0.1% error (dark) where this is smaller. All simulations used 2LPT initial conditions.

APPENDIX A: DIFFERENCE AND SUM EQUATIONS
As in (34-36), the component fluid equations (20-21) can be rewritten using ∞ as time variable and = / ∞ as velocity: for ∈ {c, b} with = Φ/( 0 ∞ ) and ∞ defined in (33). The initial conditions at ∞ = 0 must be c = b = −∇ for (A1) not to diverge. Taking the difference of (A1) for = b and = c gives where bc = b − c . Notice that the neutrino contribution contained in ∇ has dropped out. Consequently, we obtain results analogous to the ΛCDM case without massive neutrinos ).