Constraining Neutrino Cosmologies with Nonlinear Reconstruction

Nonlinear gravitational evolution induces strong nonlinearities in the observed cosmological density fields, leading to positive off-diagonal correlations in the power spectrum covariance. This has caused the information saturation in the power spectrum, e.g., the neutrino mass constraints from the nonlinear power spectra are lower than their linear counterparts by a factor of $\sim2$ at $z=0$. In this paper, we explore how nonlinear reconstruction methods improve the cosmological information from nonlinear cosmic fields. By applying nonlinear reconstruction to cold dark matter fields from the Quijote simulations, we find that nonlinear reconstruction can improve the constraints on cosmological parameters significantly, nearly reaching the linear theory limit. For neutrino mass, the result is only $12\%$ lower than the linear power spectrum, i.e., the theoretical best result. This makes nonlinear reconstruction an efficient and useful method to extract neutrino information from current and upcoming galaxy surveys.


INTRODUCTION
The large-scale structure of the Universe as traced by galaxy clustering provides one of the best windows on fundamental physics including neutrino masses, dark energy, inflation and gravity.The ongoing and future galaxy redshift surveys have the potential of greatly improving our constraints on cosmology and fundamental physics (e.g., DESI DESI Collaboration et al. 2016;Euclid Amendola et al. 2018;Euclid Collaboration et al. 2020;SPHEREx Doré et al. 2014;Subaru Prime Focus Spectrograph Takada et al. 2014; Vera Rubin Observatory/LSST LSST Science Collaboration et al. 2009; Roman Space Telescope/WFIRST Akeson et al. 2019;MegaMapper Schlegel et al. 2022).To fully exploit these powerful surveys, it is necessary to understand how to optimally extract information from the observed galaxy density fields.
Cosmological constraints derived from galaxy power spectrum are connected to extracting information from linear (or mildly nonlinear) modes, for which the coupling with other modes can be assumed be to subdominant and have evolved independently from the early hmzhu@nao.cas.cnstages of the Universe (Ferraro et al. 2022).However, at low redshifts and on small scales, nonlinear evolution leads to the coupling of modes on different scales and typically causes significant correlations in the smallscale power spectrum (Meiksin & White 1999;Rimes & Hamilton 2005, 2006).These positive off-diagonal correlations in the power spectrum covariance lead to the information saturation in the power spectrum (Rimes & Hamilton 2005, 2006;Villaescusa-Navarro et al. 2020;Bayer et al. 2021), i.e., the cosmological information does not increase when increasing the wavenumber in the nonlinear regime.For example, the constraints on neutrino mass from the nonlinear power spectra are lower than their linear counterparts by a factor of 2 at redshift z = 0 (Bayer et al. 2022).
The constraints on cosmological parameters including neutrino mass scale approximately linearly with the primordial physics Figure of Merit (FoM), which counts the total number of linear modes accessible to a galaxy redshift survey (see Sailer et al. 2021;Ferraro et al. 2022, for more discussions).Nonlinear structure formation erases the correlation with the initial conditions and leads to a loss of information, which degrades the FoM of low redshift surveys significantly.Therefore, recovering more linear modes from the observed nonlinear fields can im-prove the FoM and cosmological constraints greatly and thus maximize the scientific return of future surveys.
Reconstruction has been an efficient technique to restore the baryonic acoustic oscillation (BAO) feature in the linear power spectrum (Eisenstein et al. 2007) and has been widely used to analyse the galaxy survey data (e.g., SDSS Padmanabhan et al. 2012, BOSS Alam et al. 2017;eBOSS Alam et al. 2021).In standard reconstruction, the linear displacement field estimated from the observed nonlinear density under the Zel'dovich approximation is used to move galaxies, which reverses the large-scale bulk flows and partially restores the linear BAO signal (Eisenstein et al. 2007).Recently nonlinear reconstruction algorithms have been developed to further improve upon the standard reconstruction by solving the full nonlinear displacement field (Zhu et al. 2017;Schmittfull et al. 2017;Shi et al. 2018).The reconstruction methods works extraordinarily well for reducing the uncertainty on the BAO scale.Slightly more generally, the reconstruction can restore information about the linear modes that were corrupted by nonlinear structure formation and a higher primordial FoM can be obtained.It has been shown that nonlinear reconstruction reduces the off-diagonal components of power spectrum covariance prominently and increase the information content of power spectrum amplitude substantially (see e.g., Pan et al. 2017).Recently, Wang et al. (2022) have explored the improvement of cosmological constraints using the standard reconstruction algorithm, while Flöss & Meerburg (2023) have studied improving constraints using neural network based reconstruction, in particular primordial non-Gaussianity.
In this paper, we explore how nonlinear reconstruction improves the constraints on cosmological parameters, with a focus on neutrino mass.This improves upon the previous study (Pan et al. 2017), which only considered one parameter, the power spectrum amplitude.We use the Quijote simulations (Villaescusa-Navarro et al. 2020) to compute the pre-and post-reconstruction power spectrum covariance and perform a detailed Fisher analysis.We find that nonlinear reconstruction can efficiently reduce the mode-coupling effects in the nonlinear field and improve the cosmological constraints significantly, almost reaching the linear theory limit, which is an important step towards optimal extraction of neutrino mass information on small scales and has profound implications for future galaxy redshift surveys.
This paper is organized as follows.In Section 2, we describe the reconstruction algorithms and Fisher matrix.In Section 3, we present the numerical implementation of the Fisher analysis.In Section 4, we show the numerical results.We discuss and conclude in Section 5.

METHOD
Neutrinos with finite mass suppress the small-scale power spectrum below the neutrino free-streaming scale (Bond et al. 1980;Lesgourgues & Pastor 2006).The suppression is proportional to the sum of neutrino masses, thus allowing us to measure or constrain the total neutrino mass M ν .Detecting the minimum neutrino mass from the effects of neutrinos on the total matter δ m and cold dark matter (CDM)+baryon perturbations δ cb is one of the key science goals of the near-term and future galaxy surveys.
For galaxy clustering in the presence of massive neutrinos, the now-standard approximation that galaxies trace the CDM+baryon density field δ cb has been shown to be an excellent approximation well into the quasi-linear regime (Castorina et al. 2015).The galaxy overdensity is given by δ g = F [δ cb ] and at linear order δ g = bδ cb , where b is the linear galaxy bias.Recently, Bayer et al. (2022) demonstrated that all the information of neutrinos in the CDM field is in the linear power spectrum set by the initial conditions for k ≲ 1 hMpc −1 using the phasematched N -body simulations.Thus, galaxy clustering, which is based on the δ cb field will have negligible constraining power on neutrino masses beyond that which exists in the linear density field over the same range of scales (Bayer et al. 2022).They have found similar conclusions for other parameters such as the total matter density Ω m and primordial power spectrum amplitude A s .While the late-time linear power spectrum is not an observable, it provides a bound on the error of neutrino masses available from galaxy clustering.However, the nonlinear power spectrum does not reach this bound, because of the positive off-diagonal components in its covariance which degrades the constraints substantially.The resulting constraining power of the nonlinear power spectrum for δ cb is a factor of 2 lower than the linear one (Bayer et al. 2022).
Therefore, there is still room for improving neutrino mass constraints from galaxy clustering and potentially a factor of 2 better constraints can be achieved.The goal of this paper is to investigate whether nonlinear reconstruction can recover the information from the nonlinear δ cb field to the linear theory limit.
In this section, we review the reconstruction algorithms including standard and nonlinear reconstruction and describe the Fisher analysis of cosmological information.

Reconstruction
Reconstruction is an efficient way to reduce nonlinearities induced by the large-scale bulk flows.The standard reconstruction has been first proposed to reconstruct the linear BAO signal (Eisenstein et al. 2007).The algorithm devised by Eisenstein et al. (2007) consists of the following steps: • Smooth the nonlinear density field with a Gaussian kernel W R (k) of scale R, to filter out small-scale modes.
• Compute the shift s from the smoothed density field using the Zel'dovich approximation, (1) • Move the galaxies by s to compute the "displaced" density field, δ d (x).
• Shift an spatially uniform distribution of particles by s to form the shifted density field, δ s (x).
• The reconstructed density field is defined as The Zel'dovich approximation only captures the linear displacement on large scales, limiting the efficacy of standard reconstruction to relatively large scales, k ∼ 0.2 hMpc −1 .Recent new reconstruction algorithms have been presented to solve the nonlinear displacement field (Zhu et al. 2017;Schmittfull et al. 2017;Shi et al. 2018), which improve upon the linear standard reconstruction with a substantially improved small-scale performance.While these methods are different in the detailed numerical implementations, the underlying physical principles are the same, i.e., solving the full displacement from the mass conservation relation.In the following analysis, we use the algorithm of Schmittfull et al. (2017), solving the displacement field iteratively using the fast Fourier transform, which is efficient and easy to apply.The reconstruction algorithm devised by Schmittfull et al. (2017) consists of the following steps: • Iteratively displace objects as follows: 1. From the catalog, compute the density field on the grid.
2. Smooth the density field with the Gaussian window W R (k), with smoothing scale in the nth iteration step and truncate smallscale modes by setting δ where k Ny is the Nyquist frequency.

Compute the negative Zel
(3) 4. Move the object by the displacement s and update the object's position.
• Compute the total displacement of each object χ( qend ) = qend − x start , where x start is the starting position before the first iteration step and qend is the end position of each object after iteration (x start is the Eulerian position of each object and qend is the estimated Lagrangian position).
• Paint the estimated displacement χ( qend ) defined on the estimated Lagrangian position qend to a regular grid and truncate the small-scale modes beyond k Ny .
• The reconstructed density field from nonlinear reconstruction is An improved estimation of the small-scale displacement allows us to recover more small-scale linear modes.The cross-correlation coefficient of the reconstructed field with the linear initial conditions is higher than 0.5 for k ≲ 1 hMpc −1 (Zhu et al. 2017;Schmittfull et al. 2017;Shi et al. 2018).Note that such an improvement can also be achieved by combining the standard reconstruction and convolutional neural networks (see Shallue & Eisenstein 2022;Flöss & Meerburg 2023).

Fisher analysis
To quantify the cosmological information in the power spectrum of pre-and post-reconstruction fields, we use the Fisher information matrix (Fisher 1925).The Fisher matrix is defined as where P i = P (k i ) is the power spectrum of the ith bin and C ij is the power spectrum covariance matrix, and θ a represents the cosmological parameters (Tegmark 1997).
We have neglected the dependence of covariance on the cosmological parameters to avoid underestimating the errors (Carron 2013).
The inverse of the Fisher matrix gives the covariance matrix of the parameters and the expected errors on the parameters are equal to the square root of the diagonal elements of which represents the minimum error achievable.

SIMULATIONS
To compute the covariance matrix of the power spectrum and the response of power spectrum with respect to the cosmological parameters, we use the Quijote simulations (Villaescusa-Navarro et al. 2020), run with the TreePM code Gadget-III, an improved version of Gadget-II (Springel 2005).Each simulation evolves 512 3 CDM particles (plus 512 3 neutrino particles for nonzero neutrino mass) in a cosmological volume of 1 (h −1 Gpc) 3 .The Quijote simulations have been widely used to quantify the information content of nonlinear statistics (e.g., Villaescusa-Navarro et al. 2020;Hahn et al. 2020;Massara et al. 2021;Bayer et al. 2021;Paillas et al. 2022).More details of the Quijote simulations can be found in Villaescusa-Navarro et al. (2020).
We consider six cosmological parameters, Ω m , Ω b , h, n s , σ 8 , and M ν .We use 5, 000 simulations of the fiducial cosmology to compute the covariance matrix.We use 1, 000 simulations to compute the derivatives with respect to each parameter.The detailed information of simulations we used in the analysis is shown in Table 1.We have performed a detailed convergence test and verified that our results are converged, i.e., the constraints do not change much if the power spectrum covariance and derivatives are evaluated with less simulations (See Appendix A for more details).
We use the snapshot of CDM particles from the simulation output at redshift z = 0.Here we assume baryons trace CDM.We compute the nonlinear δ cb field on a 512 3 grid using the cloud-in-cell interpolation scheme.For standard reconstruction, we use a Gaussian smoothing of scale R = 10 h −1 Mpc, W R (k) = e −(kR) 2 /2 following Eisenstein et al. (2007).For nonlinear reconstruction, we use R init = 10 h −1 Mpc, ϵ R = 0.5, and R min = 1 h −1 Mpc.We use eight iteration steps in the reconstruction (Schmittfull et al. 2017).

Covariance matrix
We estimate the covariance matrix using 5, 000 simulations of the fiducial cosmology as where ⟨•⟩ denotes the ensemble average over simulations.
The covariance matrix for fields δ cb , δ r and δ nr are estimated numerically using the Quijote simulations.
For the linear density field (a Gaussian random field), the power spectrum covariance matrix is diagonal, where N ki is the number of Fourier modes in the corresponding k i bin.This can be directly computed us-ing the linear power spectrum at z = 0 from the linear Boltzmann code CAMB (Lewis et al. 2000).
Since the covariance matrix is estimated from simulations, it is itself a random variable with errors.The matrix inversions lead to biased constraints on the parameters.To account for this, we apply a correction to the covariance matrix (Percival et al. 2022;Paillas et al. 2022) where n sim is the number of simulations to calculate the covariance matrix, n θ is the number of parameters in Fisher analysis, n bins is the number of power spectrum bins.In our case, we have 5, 000, n θ = 6, and n bins depends on the maximum wavenumber k max included in the analysis.We have checked that this correction factor is less than 1% for k max < 0.5 hMpc −1 .

Derivatives
For the parameters Ω m , Ω b , h, n s and σ 8 , we use a central difference scheme to compute the derivatives, where P (k) is the measured power spectrum from simulations and only the considered parameter is varied and other parameters are fixed.We have used 500 pairs of simulations to calculate each derivative using Eq. ( 10).
For neutrino mass M ν , the second term in the numerator of Eq. ( 10) corresponds to a cosmology with negative neutrino mass, since the fiducial value of the neutrino mass is at the boundary of the allowed region, i.e., M ν = 0. Therefore, we use a first-order forward difference scheme to compute the derivative, where the fiducial neutrino mass M ν = 0. We use 1,000 simulations to compute the derivative with respect to the neutrino mass.

RESULTS
In this section, we present the reconstruction results, including the cross-correlation coefficient, power spectrum covariance matrix, and cosmological parameter constraints.

Cross-correlation coefficient
To quantify the reconstruction performance, we use the cross-correlation coefficient, where P δ A δ A (k) and P δ B δ B (k) are the power spectrum of fields δ A (k) and δ B (k), and P δ A δ B (k) is the cross-power spectrum between the two fields.
In Figure 1, we plot the cross-correlation coefficient The cross-correlation coefficient with the linear initial conditions for the nonlinear δ cb field, the reconstructed fields using standard reconstruction δr and nonlinear reconstruction δnr.While both reconstruction algorithms improve the correlation with the linear initial conditions substantially, the nonlinear method has a better performance on smaller scales, k ≳ 0.2 hMpc −1 .
with the linear initial conditions δ L for the nonlinear δ cb field, the reconstructed fields using standard reconstruction δ r and nonlinear reconstruction δ nr .We note that the both reconstruction algorithms improve the correlation with the linear density field substantially, es-pecially at small scales.While being competitive with nonlinear reconstruction at large scales k ≲ 0.2 hMpc −1 , the standard reconstruction has a lower correlation with the linear field δ L at smaller scales k ≳ 0.2 hMpc −1 , which needs to be improved by solving the small-scale displacement or using the convolutional neural network (e.g.Zhu et al. 2017;Schmittfull et al. 2017;Shi et al. 2018;Shallue & Eisenstein 2022).Nevertheless, we expect that both reconstruction methods can recover more small-scale linear modes, i.e., a much higher primordial physics FoM, and therefore improve the cosmological constraints.
Note that reconstruction algorithms based on solving the full nonlinear displacement have a similar performance with correlation coefficient r(k) ≳ 0.5 at k ≲ 1 hMpc −1 (Zhu et al. 2017;Schmittfull et al. 2017;Shi et al. 2018).To reduce the computational cost, we have used the iterative algorithm of Schmittfull et al. (2017).

Full shape reconstructed power spectrum
Figure 2 shows the measured power spectra divided by the linear initial power spectrum linearly scaled to the redshift z = 0 of the simulation.The original nonlinear power and linear power agree within 5% at k < 0.15 hMpc −1 , and the reconstructed power spectrum agrees with the linear power spectrum within 5% at k < 0.16 hMpc −1 .Nonlinear reconstruction slightly extends the k range where the linear theory is valid.As we have noticed in Fig. 1, the cross-correlation coefficient r becomes smaller than unity on smaller scales, we would have to add corrections to the linear power spectrum to model the power spectrum after reconstruction.Power spectra from the simulations at z = 0, divided by the linear initial power spectrum linearly scaled to z = 0.The reconstructed power spectrum agrees with the linear power spectrum within 5% at k < 0.16 hMpc −1 , while the original nonlinear power and linear power agree within 5% at k < 0.15 hMpc −1 .Compared to the nonlinear density without reconstruction, nonlinear reconstruction slightly improves the agreement with the linear linear power spectrum.However, the best-fit theoretical model agrees with the reconstructed power spectrum within 1% at k ≲ 0.24 hMpc −1 and 3% at k ≲ 0.5 hMpc −1 .
The fact that nonlinear reconstruction misses smallscale broadband power relative to the linear power is not necessarily a problem.All that is needed to estimate cosmological parameters is a theoretical model describing the reconstructed power spectrum as a function of cosmological parameters.Following the literature on modeling the reconstructed power spectrum (see e.g.Noh et al. 2009;Padmanabhan et al. 2009;White 2015), we use a simple theoretical model to fit the power spectrum after reconstruction, where the propagator C(k) = e −k 2 Σ 2 /2 describes the damping of the linear power spectrum, P δ L is the linear initial power spectrum linearly scaled to redshift z = 0, and the noise power spectrum P N quantifies the nonlinearities after reconstruction.We then fit the two parameters Σ and P N by minimizing the sum of squares where the hat denotes the reconstructed power spectrum estimated from simulations and k i denotes the power spectrum in the ith bin.Here, we directly use the measured linear initial power spectrum Pδ L when computing the theoretical model, to reduce the impact of variations due to the cosmic variance.The maximum wavenumber included in the fitting is k max = 0.8 hMpc −1 .We use equal weights for all k bins to avoid over-fitting small scales.The best-fit parameters are Σ = 1.67 h −1 Mpc and P N = 38 [(h −1 Mpc) 3 ].The best-fit theoretical model agrees with the reconstructed power spectrum within 1% at k ≲ 0.24 hMpc −1 and 3% at k ≲ 0.5 hMpc −1 , which makes it possible to perform cosmological parameter inference at the percent level using the reconstructed power spectrum.However, when applying to galaxy surveys, more complexities have to be considered such as galaxy biasing and redshift space distortions, which will be explored in the future.On smaller scales, k ≳ 0.5 hMpc −1 , the model predicts a higher noise level than the data.This is partly due to that the simple white noise term misses the scaledependence of the residual nonlinearities.Additional improvements on small scales may be possible by using the effective field theory of large-scale structure.On small scales we would have to add more counterterms to this simple model to match the reconstructed power spectrum (e.g., Schmittfull et al. 2017).We leave this for future work.
In the following analysis, we use the reconstructed power spectrum in the simulation to compute the Fisher information.

Covariance matrix
We apply two reconstruction algorithms to the 5, 000 simulations of fiducial cosmology and compute the covariance matrix of the power spectrum.
In Figure 3, we plot the correlation matrix, i.e., the normalized covariance matrix, defined as C ij / C ii C jj , where C ij is the covariance matrix.For the nonlinear δ cb field, the correlation matrix is almost diagonal on large scales, k ≲ 0.1 hMpc −1 , and exhibits significant offdiagonal correlations on smaller scales, k ≳ 0.1 hMpc −1 , which leads to the information saturation of the nonlinear power spectrum (e.g., Rimes & Hamilton 2005, 2006;Villaescusa-Navarro et al. 2020;Bayer et al. 2021).The standard reconstruction reduces the off-diagonal elements and the covariance matrix is nearly diagonal at k ≲ 0.2 hMpc −1 , in agreement with the cross-correlation coefficient observed in Figure 1.However, the covariance matrix of nonlinear reconstruction is diagonal well into the nonlinear regime, k ≲ 0.8 hMpc −1 , where the coupling with other modes are subdominant.This nearly diagonal covariance matrix for the δ nr field indicates more linear information available than the power spectra of the δ cb and δ r fields.
With enough simulations, we have not observed large negative off-diagonal elements before and after recon- )) Figure 3. Correlation matrix for the power spectrum of the nonlinear δ cb field, and the reconstruction fields δr and δnr using standard and nonlinear reconstruction methods, respectively.The covariance matrix is estimated using 5, 000 simulations.
Reconstruction reduces the off-diagonal correlation substantially.The density modes are nearly independent at k ≲ 0.8 hMpc −1 for the nonlinear reconstruction method.
struction.There only exists a few negative off-diagonal elements with very small absolute values ∼ 10 −2 , mainly due to the random statistical fluctuations.Therefore, our results should be more robust than the previous work which uses less realizations to estimate the covariance matrix (Pan et al. 2017).We have confirmed that our results have converged, i.e., the parameter constraints do not change much if less simulations are used.See Appendix A for more details of the convergence tests.

Cosmological constraints
In Figure 4, we show the marginalized 1σ constraints on cosmological parameters, Ω m , Ω b , h, n s , σ 8 , and M ν , as a function of k max .On large scales where the cosmic density field evolved linearly from the initial conditions, all fields, δ L , δ cb , δ r , and δ nr , have a similar constraint on the parameters at k ≲ 0.1 hMpc −1 .The constraints from the nonlinear δ cb start to saturate at k max = 0.2 hMpc −1 , especially for σ 8 and M ν , for which the dominant information comes from the power spectrum amplitude.For standard reconstruction, the constraints are much better than the nonlinear δ cb field.The information from δ r continue to increase for k ≳ 0.2 hMpc −1 and saturates at k max ∼ 0.3 hMpc −1 .This is due to the limited validity of the Zel'dovich approximation adopted in the standard reconstruction.However, the constraining power of the δ nr field is even better than δ r and the constraints obtained using nonlinear reconstruction continue to improve until k max = 0.5 hMpc −1 , nearly reaching the linear theory constraints.The improvement at k ≳ 0.3 hMpc −1 upon the standard reconstruction is mainly because a better estimation of small-scale displacement using nonlinear methods removes more shift nonlinearities, which is not captured by the Zel'dovich approximation.
We notice that the constraints on neutrino mass M ν with reconstruction are almost the same as the linear field at k ≲ 0.25 hMpc −1 and the constraint still nearly reaches the linear theory limit at k ≳ 0.25 hMpc −1 for the nonlinear reconstruction method.This result is remarkable, since the δ cb field and derived quantities, e.g., the halo field, contain little information about neutrino mass beyond that which exists in the linear power spectrum for k ≲ 1 hMpc −1 , demonstrated by using phasematched simulations (Bayer et al. 2022).Therefore, the linear power spectrum provides a bound on the error that can be obtained using δ cb and it is exciting to see that nonlinear reconstruction can recover nearly all the information regarding neutrino mass from the nonlinear δ cb field.
The M ν −σ 8 degeneracy has limited the neutrino mass constraint on small scales (Villaescusa-Navarro et al. 2020;Bayer et al. 2021), which could not be reduced using power spectrum alone (See Figure 8 of Appendix B).Combining other nonlinear statistics such as halo mass function and void size function can tighten the constraint on neutrino mass by breaking the degeneracy between different parameters (e.g.Bayer et al. 2021), which will be subject of further inquiry.
In Table 2, we list the marginalized 1σ errors for k max = 0.5 hMpc −1 .We find that nonlinear reconstruction can achieve a significant level of improvement on all 6 parameters and the improvements to be a factor of 2. 72, 2.65, 2.64, 2.65, 2.49, 2.54 for Ω m , Ω b , h, n s , σ 8 and M ν , respectively.Thus, the constraining power of the reconstructed power spectrum P δnr , i.e., the signal to noise ratio of neutrino mass measurement ∝ 1/σ(M ν ), is . The marginalized 1σ error for each cosmological parameters as a function of kmax.The parameter constraints from the nonlinear δ cb start to saturate at kmax = 0.2 hMpc −1 .The information from standard reconstruction continues to increase until kmax ∼ 0.3 hMpc −1 , limited by the linear assumption when estimating the displacement field.The constraining power of the δnr field is much better than δ cb and δr and continues to improve until kmax ∼ 0.5 hMpc −1 .

Marginalized Fisher Constraints Probes
Ωm Table 2. Marginalized 1σ errors of cosmological parameters for kmax = 0.5 hMpc −1 .We list the constraints obtained from the power spectrum of δ cb , δr, δnr, and δL, and the multiplicative improvement in the constraints with nonlinear reconstruction, compared to those from the power spectrum of δ cb field alone.The signal to noise ratio of neutrino mass measurement, ∝ 1/σ(Mν ), is only 12% lower than the linear power spectrum, i.e., the theoretical upper limit achievable with the δ cb field.
only 12% lower than the theoretical upper limit achievable with the δ cb field.We notice that the exact values of σ(M ν ) are slightly different between our results and Bayer et al. (2022).We have confirmed that this is due to the different fiducial cosmology adopted to compute the derivatives to neutrinos.Our fiducial cosmology is for a Universe with massless neutrinos, while Bayer et al. (2022) have used a fiducial neutrino mass M ν = 0.05 eV when computing the linear theory constraints.Nevertheless, we find that the relative improvements are not sensitive to these numerical details and thus this does not impact our major conclusions presented above, i.e., nonlinear reconstruction can recover almost all the neutrino information, nearly reaching the theoretical upper limit.3. Multiplicative improvement using the reconstructed power spectrum and the theoretical model for kmax = 0.5 hMpc −1 .

Test with the theoretical model
In the above analysis, we have used the reconstructed power spectrum from simulations for computing the covariance matrix and derivatives to cosmological parameters.This requires an accurate modeling of the reconstructed power spectrum, which might be achieved with the emulator constructed from simulations (see e.g., Wang et al. 2023, for a recent exploration).
When using the perturbation theory approach, usually we need to use more parameters to describe the nonlinear effects in the reconstructed power spectrum, such as the simple model we introduced earlier.The constraints on cosmological parameters will degrade after marginalizing over these nuisance parameters due to possible degeneracy between parameters.We can employ the simplified model in Eq. ( 13) to calculate the Fisher information.In this approach, we utilize the reconstructed power spectrum covariance matrix in the simulations, while using Eq. ( 13) for derivative computations.In this way, we include the damping scale Σ and noise power spectrum P N in the Fisher analysis.
In Table 3, we list the multiplicative improvement using the reconstructed power spectrum from simulations and the theoretical model for k max = 0.5 hMpc −1 .We see that the improvement of the constraint on neutrino mass is ∼ 20% lower than using the power spectrum from simulation.However, nonlinear reconstruction can still achieve a substantial improvement, by a factor of 2, compared to the original nonlinear power spectrum.The constraints on other parameters are not as good as the neutrino mass, which is due to the degeneracy between the cosmological and nuisance parameters, especially the spectrum index n s , since both the spectrum index and the damping scale change the power spectrum shape.

DISCUSSION AND CONCLUSION
In this paper, we apply nonlinear reconstruction to the CDM density fields from simulations and find the parameter constraints are much tighter than those from the power spectrum of δ cb directly, with the improvement to be a factor of ∼ 2 for each parameter.A salient conclusion is that the constraint on neutrino mass can nearly reach that of the linear power spectrum, with a signal-to-noise ratio only 12% lower than the theoretical best result, making reconstruction an efficient way to extract cosmological information from the nonlinear galaxy survey data.
We have assumed baryons trace CDM when using the CDM field as an approximation of the CDM+baryon distribution.This is a valid approximation on large and intermediate scales, k ≲ 0.5 hMpc −1 , where the baryonic effects are subdominant (see e.g.Chisari et al. 2019;Villaescusa-Navarro et al. 2021).The analysis with δ cb corresponds to an unbiased tracer with negligible shot noise, δ g = bδ cb with b = 1.Therefore, our constraints on M ν should be interpreted as maximum information content in principle.Higher order biasing contributes tens of percent of the total galaxy power spectrum at k ∼ 0.5 hMpc −1 (see e.g., Schmittfull et al. 2019) and the importance of shot noise depends on the specific surveys.In future work, we plan to extend the analysis to simulated galaxy samples, which are more representative of real data.
In redshift space, more information regarding neutrinos is available, since the redshift space distortions is induced by the peculiar velocity, which is sourced by the total matter field and thus is sensitive to the neutrino overdensity δ ν (Bayer et al. 2022).Therefore, it is possible to obtain better constraints in redshift space, which we plan to explore in future.
To obtain an unbiased inference of cosmological parameters, it is necessary to have an accurate power spectrum model for reconstruction.The modeling of standard reconstruction has been explored using both the Lagrangian perturbation theory (Noh et al. 2009;Padmanabhan et al. 2009;White 2015;Seo et al. 2016;Chen et al. 2019) and standard perturbation theory (Hikage et al. 2017(Hikage et al. , 2020)).Recently, Ota et al. (2021Ota et al. ( , 2022) ) have explored the modeling of nonlinear displacement using Lagrangian perturbation theory.The simulationbased method such as power spectrum emulator might be needed on nonlinear scales, where the perturbation theory methods are no longer available.
While for the current galaxy samples such as BOSS (Alam et al. 2017) and eBOSS (Alam et al. 2021), the nonlinear reconstruction methods have a similar efficiency as standard reconstruction because of the higher shot noise.The near term and upcoming surveys such as DESI (DESI Collaboration et al. 2016), Euclid (Amendola et al. 2018;Euclid Collaboration et al. 2020), MegaMapper (Schlegel et al. 2022), etc, will generally have a much higher galaxy number density and therefore nonlinear reconstruction can achieve a better performance with the lower shot noise.We expect that reconstruction will be a promising method to extract cosmological information from future galaxy surveys.
We acknowledge Adrian Bayer and Gong-Bo Zhao for enlightening discussions about the In this appendix, we present the convergence tests of our results.In Figure 5, we show the marginalized 1σ errors on cosmological parameters as a function of the number of simulations used to estimate the covariance matrix for nonlinear reconstruction.The results are presented as the ratio of the errors computed using N cov simulations and 5000 simulations.The derivatives are computed using 500 pairs of simulations.From left to right, we show the convergence of marginalized constraints for different k max .We see that the results converge well for all scales from 0.05 hMpc −1 to 1 hMpc −1 , with a small error at percent level.In Fig. 6, we show the derivatives of power spectrum with respect to cosmological parameters.The derivatives for the power spectrum show good convergence for 500 pairs of simulations.
Figure 7 shows the marginalized 1σ errors with nonlinear reconstruction on the cosmological parameters as a function of the simulation pairs to compute the derivatives.The results are presented as the ratio between the errors calculated with N Deri pairs and 500 pairs of simulations.The covariance matrix is computed using 5000 simulations.For the derivatives, a good convergence is achieved for k max from 0.06 hMpc −1 to 0.5 hMpc −1 .For k max = 0.05 hMpc −1 , the fractional error is about 20% when using 500 simulation pairs.There are not enough modes on these large scales to get converged results with the simulations available.A higher accuracy power spectrum estimation is required to obtain reliable results at k max = 0.05 hMpc −1 , but it is fine since the dominant constraining power comes from small-scale modes.For k max = 0.06 hMpc −1 , the results become much more stable.The fractional error is under 2% when using 500 simulation pairs to obtain partial derivatives.Therefore, we present the Fisher constraints from k max = 0.06 hMpc −1 to k max = 0.5 hMpc −1 as our default results.

B. CONFIDENCE INTERVALS
In Figure 8, We show the two-dimensional 68% and 95% confidence intervals from the Fisher analysis for k max = 0.5 hMpc −1 .The constraints are generally much tighter after reconstruction.However, it does not break the strong degeneracy between σ 8 and M ν using the reconstructed power spectrum alone.By introducing other nonlinear statistics, a much tighter constraint can be achieved (Bayer et al. 2021).
Figure2.Power spectra from the simulations at z = 0, divided by the linear initial power spectrum linearly scaled to z = 0.The reconstructed power spectrum agrees with the linear power spectrum within 5% at k < 0.16 hMpc −1 , while the original nonlinear power and linear power agree within 5% at k < 0.15 hMpc −1 .Compared to the nonlinear density without reconstruction, nonlinear reconstruction slightly improves the agreement with the linear linear power spectrum.However, the best-fit theoretical model agrees with the reconstructed power spectrum within 1% at k ≲ 0.24 hMpc −1 and 3% at k ≲ 0.5 hMpc −1 .

Figure 5 .
Figure 5. Marginalized 1σ error as a function of the number of simulations used to calculate the covariance matrix for nonlinear reconstruction.From left to right, we show the results for different kmax.For Ncov = 5000, a tight convergence is achieved.

Figure 6 .Figure 7 .
Figure6.The derivatives of the power spectrum to cosmological parameters.The derivatives are well converged using 500 pairs of simulations.

Table 1 .
The details of the Quijote simulations used in this work.The fiducial cosmology contains 5, 000 simulations, which are used to estimate the power spectrum covariance matrix.Other simulations with different cosmological parameters are used to compute the partial derivatives of the power spectrum to cosmological parameters, where only the considered parameter is varied and other parameters are fixed.The initial conditions are generated at z = 127 using the second-order Lagrangian perturbation theory, except for the neutrino simulations where the Zel'dovich approximation is used.