Reconstructing the Long-wavelength Matter Density Fluctuation Modes from the Scalar-Type Clustering Fossils

Revealing the large-scale structure from the 21cm intensity mapping surveys is only possible after the foreground cleaning. However, most current cleaning techniques relying on the smoothness of the foreground spectrum lead to a severe side effect of removing the large-scale structure signal along the line of sight. On the other hand, the clustering fossil, a coherent variation of the small-scale clustering over large scales, allows us to recover the long-wavelength density modes from the off-diagonal correlation between short-wavelength modes. In this paper, we study the requirements for an unbiased and optimal clustering-fossil estimator and show that (A) the estimator is unbiased only when using an accurate bispectrum model for the long-short-short mode coupling and (B) including the connected four-point correlation functions is essential for characterizing the noise power spectrum of the estimated long mode. The clustering fossil estimator based upon the leading-order bispectrum yields an unbiased estimation of the long-wavelength ($k\lesssim 0.01~[h/{\rm Mpc}]$) modes with the cross-correlation coefficient of $0.7$ at redshifts $z=0$ to $3$.


Introduction
The intensity mapping of the 21cm hyper-fine transition line (21cmIM) can be a powerful probe of the large-scale structure of the Universe [1,2].With the 21cm line's low opacity and spectroscopic nature, the 21cmIM will provide a three-dimensional map of neutral hydrogen distribution over the unprecedently large volume.The ongoing and future 21cmIM projects, such as Tianlai [3], CHIME [4], HIRAX [5], BINGO [6], and SKA [7], are designed to observe the large-scale structure at redshifts z = 0.5 − 6, to measure, for example, the Hubble expansion rate and the angular diameter distance from the baryon acoustic oscillation (BAO) [8].
One of the main challenges of these 21cmIM experiments is the contamination from the extremely loud foreground from synchrotron, free-free, and thermal dust emissions.For example, the dominating galactic synchrotron emission is more than five orders of magnitude louder compared to the targeted large-scale structure signal in 21cmIM [9,10].The standard foreground cleaning method [11][12][13][14][15][16][17][18][19][20][21] takes advantage of the fact that synchrotron radiation is smooth in the frequency domain, while the 21cmIM signals show genuine cosmological density fluctuations along the radial direction.After the foreground cleaning, however, the long-wavelength Fourier modes parallel to the line of sight become inaccessible [10,22].Furthermore, the frequency-dependence response of the instrument can lead to a foreground wedge contamination in the Fourier space for k ∥ < k ⊥ tan ψ, where the angle ψ increases with the field of view [23][24][25][26][27][28][29].
These long-wavelength modes plagued by foreground cleaning are usually in the linear regime where the observed galaxy clustering statistical can be modeled by the Kaiser formula [30,31], or near horizon scales where general relativistic corrections are important (see Ref. [32] for a review).For both cases, the theoretical model is well understood in the framework of linear perturbation theory, and the accurate measurement of these long-wavelength modes will be translated to the measurement of the growth rate of the cosmic density field and metric perturbations as well as the local-type primordial non-Gaussianities [33].Measuring these parameters, in turn, will inform us of the nature of dark energy and the physics of the primordial universe, respectively.In particular, given the large volume that 21cmIM covers, the large number of long-wavelength modes could provide an unprecedented tight constraint on cosmological parameters [34][35][36][37].Therefore, to fully exploit the information from 21cmIM, it is crucial to recover the contaminated long-wavelength Fourier modes (hereafter, "long mode," for short).
To recover the contaminated long modes, Ref. [38] have proposed an innovative technique called cosmic-tide reconstruction.The basic idea of the reconstruction is to exploit the non-Gaussianity of the density field coming from the nonlinear gravitational evolution that couples Fourier modes of different wavelengths.Namely, knowing the details of this coupling should allow us to infer the long-wavelength Fourier modes from the short-wavelength Fourier modes (hereafter, "short mode," for short).Based on this idea, Ref. [39] developes a quadratic estimator for the long modes based on the anisotropic modulation of short modes' power spectrum by the large-scale tidal field.The coupling coefficient in the estimator is determined by the leading-order tidal interaction [40].Applying this quadratic estimator to the N -body simulation results shows that the cross-correlation coefficient between the reconstructed long mode and the original long mode can reach 0.9 at k < 0.1 h/Mpc.This implies that the phase of the reconstructed long mode is highly correlated to the true long mode.Ref. [41] has further confirmed this result.Additional insights into tide reconstruction from tracers in real and redshift space are provided by [42] and [43].Also see [44] for its application to recovering missing radial long modes in 21cm intensity mapping surveys, and [45] in the context of lensing reconstruction of line intensity mapping.
While following the spirit of cosmic-tide reconstruction, this paper investigates a different mathematical framework, called clustering fossil, for reconstructing the long modes.The clustering fossil utilizes the three-point non-Gaussian correlation between long mode and two short modes, which generates, in Fourier space, the off-diagonal two-point correlators between two short modes [46].In a statistically homogeneous universe, the two-point correlation function ⟨δ * (k)δ(k ′ )⟩ in Fourier space only has diagonal components proportional to δ D (k − k ′ ): the power spectrum.In the presence of long modes, however, the non-Gaussian coupling between short modes and long modes makes locally measured two-point correlators differ from place to place.The statistical homogeneity is broken up locally, and the non-zero off-diagonal two-point correlators among short modes contain information about the long mode coupling with them.As we show in Section 2.2, this off-diagonal correlation is proportional to the amplitude of the long mode and the coupling coefficient, which can be read from the squeezed-limit bispectrum.
Following Ref. [46], we construct a quadratic clustering-fossil estimator for the long modes based on the off-diagonal correlation between short modes.Unlike Ref. [46], which pursues inflationary spectator fields using the non-Gaussian coupling in the early universe [47,48]; however, we focus on the non-Gaussianities caused by nonlinear matter clustering.In particular, applying the method to the case of 21cmIM, we treat the contaminated long mode as a scalar-type fossil field.Along this line, Refs.[49,50] have applied the clustering-fossil estimator to reconstruct the long mode based on the leading-order (tree-level) bispectrum in standard cosmological perturbation theory (SPT) [51].Analyzing a suite of N -body simulations, they have recovered the long mode, with morphological features similar to the ground truth.They found, however, that the power spectrum of the recovered long mode is biased, which is consistent with the findings of Ref. [39].A further study [52] has applied this technique to the biased tracer fields to recover the matter density long mode and forecast the improved constraint on primordial non-Gaussianity through multi-tracer method, where the galaxy bias are also included into the estimator.
To proceed with the clustering-fossil-based reconstruction method, we need a more thorough theoretical understanding of the method's applicability and regime of validity.The goal of this paper is to scrutinize the method.We aim to provide a theoretical framework to explain the systematic bias observed in the N-body simulations and to characterize the statistical uncertainties of the recovered long mode.To solve the problem for the fully nonlinear density field, we need to include the coupling effects between long and short modes beyond the leading order and the non-Gaussian statistics of the short modes.To assess the requirement for a robust reconstruction method, we take a simpler approach here.Namely, instead of analyzing the result from full N -body simulations, we test the reconstruction method against a controlled sample density field only containing the first-and second-order density perturbations from the GridSPT (grid-based calculation of standard perturbation theory) [53][54][55].Because the GridSPT density field strictly follows SPT, the fossil estimator constructed from the leading-order bispectrum can fully capture the coupling between the long and short modes.Therefore, testing the clustering fossil estimator in this way allows us to disentangle the effect of higher-order coupling from other effects, such as the non-Gaussianity of short modes.To inject the effect from the higher-order nonlinear couplings, we use the 2LPT (second-order Lagrangian Perturbation Theory) because, while agreeing with the SPT to second order, the 2LPT density field also contains higher-order contributions.We take advantage of the efficiency of both GridSPT and 2LPT, which enable us to generate a large number of realizations to suppress the sampling variance and find any underlying systematic errors in the estimator.
By testing the estimator on well-controlled datasets, we find that the clustering-fossil estimator is unbiased only when an accurate bispectrum model is used for the coupling between long and short modes.Furthermore, to reconstruct the long-mode power spectrum, the connected four-point correlation function must be included to compute the noise power spectrum.It is important to note that the goal of this paper is to assess the potential systematics of the quadratic estimator, so we restrict our study to matter clustering in real space for clarity and transparency.To apply this method to galaxy surveys or 21cm intensity mapping experiments, we need to include the galaxy bias and redshift space distortions into account(see e.g.[52], [43]), which would make the accurate reconstruction even harder.
The rest of this paper is organized as follows.In Section 2, we review the clustering fossil technique and construct the optimal estimator of the large-scale matter density modes.In Section 3, we introduce GridSPT and 2LPT we use to generate the nonlinear shortwavelength modes and implement the fossil estimator to reconstruct the large-scale mode.
We compare the results with the theoretical prediction in Section 4. We conclude and discuss the future applications in Section 5. We provide the fossil estimator in the continuous limit in Appendix A.
Throughout this paper, we use the following convention of Fourier transformation, Note that we use the same character for a function f and distinguish the Fourier-space representation and configuration-space representation by the argument.For the compactness of the equations, we use the following convention for the sum of multiple vectors, 2 Estimating long modes from the clustering-fossil estimator We begin the section by motivating the clustering-fossil method [46] for a generic non-Gaussian density field in Section 2.1.In Section 2.2, we construct the optimal quadratic estimator for measuring the long mode from the imprint of the squeeze-limit bispectrum, or long-short-short coupling.We also present the power spectrum of the reconstructed long modes and the cross-correlation coefficient between the original and the reconstructed long modes.Finally, we present the reconstructed long modes' covariance matrix, or noise power spectrum, in Section 2.3.

Position-dependent correlations induced by squeezed-limit non-Gaussianity
The standard cosmological model based on the Friedman-Lemaître-Robertson-Walker world model is spatially homogeneous and isotropic.The spatial homogeneity extends to the clustering properties of galaxies, which is often referred to as statistical homogeneity.We begin this section by defining the statistical homogeneity through n-point correlation functions.The fundamental quantity describing the galaxy clustering is the density contrast δ(x) ≡ n(x)/n − 1, where n(x) is the density at location x and n is the cosmic mean density.The two-point correlation function is defined in terms of the density contrast as Here, the bracket ⟨...⟩ represents the average over the statistical ensemble of the cosmic density field.Note that ξ(r) depends only on the separation vector r between two positions x and x + r but is independent of the position x where we make the measurement.This property manifests the statistical homogeneity.Furthermore, the statistical isotropy would restrict ξ(r) = ξ(r).Equivalently, the power spectrum, the Fourier transformation of the two-point correlation function, is defined as where the Dirac-delta operator, δ D (k 12 ) on the right-hand side, encodes the statistical homogeneity.That is, the statistical homogeneity dictates that the correlation between the Fourier modes δ(k 1 ) and δ(k 2 ) vanishes unless k 1 + k 2 = 0. We call such two-point correlators in Fourier space diagonal.We can extend this concept and write the n-point correlation function and the n-poly spectra for statistically homogeneous density contrast δ as follows: (2.4) We call a density field Gaussian when the two-point correlation function is the only non-vanishing connected n-point correlation function [56]; otherwise, the density field is non-Gaussian.
The clustering-fossil estimator exploits the fact that non-Gaussianities in density fields can yield a spatial variation of correlation functions.That is, even if the density field satisfies statistical homogeneity when taking the ensemble average as in Equation (2.3), the locally measured n-point correlation function can have an explicit position dependence through the non-Gaussian correlation.This is because the non-Gaussian correlation function involving one or more long-wavelength modes (δ ℓ ) can be broken down as the multiplication of the conditional correlation functions and the conditioning long modes: where ⟨⟩| X stands for the correlator evaluated with the condition X.The expression follows from the definition of conditional probability P (A ∩ B) = P (A|B)P (B).
The most well-studied example of clustering fossil comes from the squeezed (or soft) limit of non-Gaussian correlation functions.Namely, the squeezed-limit (n + 1)-point correlation function can modulate the n-point correlation function, inducing the position-dependent power spectrum [57,58] or the position-dependent bispectrum [59].In this case, we may write the conditional correlation function as where C is a function independent from δ ℓ , and we truncate the expansion at the linear-order response.On large scales where the long mode δ ℓ is in the linear regime, we compute the original correlation function as The general idea of the clustering fossil estimator [46] is to use Equation (2.6) as a starting point for estimating the long-wavelength mode δ ℓ for the non-Gaussian density field with underlying (n + 1)-point correlation function taking the form of Equation (2.7).As shown in [46], the same argument applies when using a long mode h s ℓ with general spin s.The density field is a spin-0 case.

The Optimal Clustering Fossil estimator
Let us focus on the quadratic clustering-fossil estimator using the squeezed-limit three-point correlation function, or its Fourier transform, bispectrum.The Fourier-space counterpart of Equation (2.7), when considering the three-point correlator, becomes with kernel f (k 1 , k 2 ; k) encoding the details of the coupling between the long mode δ ℓ (k) and two short modes δ(k 1 ) and δ(k 2 ).Hereafter, we drop the explicit subscript ℓ and implicitly use k ≪ k i to indicate the long mode.We emphasize that the three wave vectors in f (k 1 , k 2 ; k) must satisfy k 12 + k = 0.As we have shown in Section 2.1, Equation (2.8) is equivalent to the following conditional two-point correlator: Here, δ K refers to the Kronecker delta that we use instead of Dirac delta, after absorbing appropriate dimensionful constants inside f (k 1 , k 2 ; k), and the superscript * is the complex conjugate.The locally broken statistical homogeneity is now manifested by the off-diagonal correlator in Equation (2.9).Note that Equation (2.8), more generally Equation (2.7), assumes that the long mode is in the linear regime where we can neglect the higher order terms in Equation (2.7); so the fossil estimator is unbiased only in that limit.A generic unbiased fossil estimator would require additional corrections to include higher-order effects.In Section 4.3.2 of this paper, however, we show that the nonlinear part of the long-mode power spectrum is less than a few percent of the cosmic variance, which ensures that the nonlinear correction is unimportant at the scale within with we implement the estimator.
Let us translate the fossil equation, Equation (2.9), to the quadratic estimator for the long mode δ(k).A simple estimator turning Equation (2.9) around would be δnaive In the second equality, we use the reality of the density field [δ(k)] * = δ(−k) and that bispectra and f (k 1 , k 2 ; k) must be real for even-parity density field: . For a given k, this naive estimator sums over all possible pairs δ(q i ) and δ(k − q i ) weighted by a known kernel function f (q i , k − q i ; −k) determined by the bispectrum as in Equation (2.8).In order to exclude the duplication of counting the same quadratic contribution twice, we demand the condition k < |k − q i | ≤ q i .Since the estimator is based on the multiplication of two short modes, it is called a quadratic estimator.A similar quadratic estimator has already been used widely in the CMB lensing reconstruction [61,62].
Noticing that the individual contribution in Equation (2.10) is subject to different variances, we can introduce a normalized weight to find an optimal estimator which minimizes the variance of δr (k): (2.12) First, we compute the covariance matrix C ij for individual contribution in the naive estimator [Equation (2.10)]: by ignoring all the connected parts of four-point correlators.Here P (q) and P (|k − q|) are the measured nonlinear power spectra of the short modes.In this case, the inverse-variance weight gives the minimum variance estimator.Using the inverse-variance weight in Equation (2.14), the optimal fossil estimator for the long mode becomes the variance of the optimal estimator δr becomes . (2.16) The minimum variance P N G (k) turns out to be the noise term of the recovered long mode's auto power spectrum, when taking only the Gaussian (disconnected) four-point correlator of the density contrast δ: Here, N k is the number of Fourier modes for a fixed k.Equation (2.17) suggests that we have to subtract the noise power spectrum P N G (k) from the auto power spectrum, in order to recover an unbiased long-mode power spectrum.
Let us turn into the cross-correlation coefficient r(k) between the recovered long mode δr (k) and the original one δ(k) that we designate with subscript m: r(k) ≡ P rm (k) Here, δr (k)δ m (k ′ ) = (2π) 3 P rm (k)δ D (k + k ′ ) is the cross power spectrum which approaches the linear power spectrum for the squeezed (k → 0) limit.In Section 3.1, we calculate the correction term to this approximation when using the full bispectrum expression.The leading contribution of the correction term comes from the correlation between the second-order of the long mode and two first-order short modes, which vanishes in the squeezed limit.
Combining all results, we calculate the leading-order expression for the cross-correlation coefficients in the squeezed limit [k → 0 and P mm → P L (k)] as, From Equation (2.20), we see that r(k) → 1 if and only if when the noise power spectrum is much suppressed compared to the linear power spectrum.The noise power spectrum [Equation (2.16)] becomes smaller as we include more and more short modes in the quadratic estimator.Meanwhile, the noise power spectrum is insensitive to the lower limit of q and |k − q| because most reconstruction power comes from the high-frequency (|q| ≫ k) modes.
In this paper, we use the short modes δ(q), δ(k − q) with wavenumbers satisfying where q max is the wavenumber of the smallest scale used for reconstruction.Ideally, we want to make q max as large as possible to suppress the noise P N G .In practice, however, we are limited by strong nonlinearities on small scales, so we can only keep q max in the quasi-linear scale where perturbation theory accurately models the nonlinearities [63,64].As we shall show below, using the short modes beyond the quasi-linear regime leads to the systematic bias in the reconstructed long modes.
As stated earlier, the summations in the equations above excludes the duplications (that is, q ′ =k − q is identical to q), which leads to a factor of two difference between Equation (2.16) here and the expressions shown in previous studies [46,49,50,52].The restriction in the summation, however, cancels the factor-of-two difference, so the expressions are identical.
However, just like our derivation in this section, most previous studies to date have ignored the connected four-point correlation function in deriving Equation (2.13), so the noise power spectrum is underestimated (see Figure 3 below).We now present the full noise power spectrum, including the connected four-point function.Readers with interest can find related discussions of connected four-point correlation function in the noise power spectrum of quadratic estimators in [45,52] for long-mode reconstruction and [65] for CIB lensing.

The Full Noise Power Spectrum
In Section 2.2, we use the Gaussian approximation for when calculating the covariance matrix of naive quadratic terms estimators in Equation (2.10).In this section, we supplement that calculation by deriving the full covariance matrix, including the connected four-point function.We also present the improved noise power spectrum by using the full covariance matrix but still adopting the inverse-variant weight using the Gaussian approximation.
We first compute the full covariance matrix of each term contributing to the naive estimator Equation (2.10) as where T (q i , k − q i , q j , −k − q j ) is the connected four-point function, or trispectrum, defined as follows: After removing the duplicated contributions by imposing the inequality in Equation (2.21), the covariance matrix becomes (2.24) From the covariance matrix, we compute the full noise power spectrum for the optimal quadratic estimator, Equation (2.15): (2.25) Note that the estimator Equation (2.15) is optimal under the Gaussian assumption but may not remain optimal when the trispectrum contributes significantly to the covariance matrix.Since the relative contribution from the trispectrum rises toward small scales [66,67], we expect that the noise power spectrum P N full (k) would saturate beyond some q max .Later, we verify the statements by comparing numerically the two noise power spectra P N full (k) and P N G (k) for the two nonlinear density realizations: 2LPT [68] and GridSPT [53].
3 Leading-order clustering-fossil estimator and nonlinear density fields The main goal of this paper is to assess the requirement for an accurate and unbiased clustering-fossil estimator in Equation (2.15).As for the simplest but realistic nonlinear density field, as shown in Section 3.1, we study the clustering-fossil estimator using the second-order density perturbations as defined in the SPT (Standard Perturbation Theory) [51].Specifically, we apply the quadratic clustering-fossil estimator to measure the long modes from the squeezed-limit bispectrum that encodes an off-diagonal power spectrum of the short modes.
Our approach is different from the previous studies [49,50], which tested the quadratic estimator based upon the second-order SPT (Equation (3.4) below) against a few N -body simulation results, or tested the second-order EFT against dark matter halo fields [52].In these studies, they have drawn affirmative conclusions about the possibility of reconstructing large-scale modes.At the same time, the quadratic estimator must fail beyond some q max , because the squeezed-limit bispectrum deviates from the tree-level SPT or EFT prediction.Also, the noise power spectrum of the reconstructed lone mode must be underestimated because the trispectrum contribution becomes important on small scales.Neglecting the trispectrum contribution would bias the power spectrum of reconstructed long modes.We cannot distinguish these two effects in N -body simulations.Also, we need multiple realizations to draw a robust statistical conclusion, but N -body simulations are too expensive for this purpose.
In contrast, our method based on the SPT allows a more systematic and in-depth study of the method because we can generate nonlinear density fields with a well-controlled nonlinear order.In this paper, we adopt a novel grid-based standard perturbation theory (GridSPT) method to generate the density fields up to second-order (Section 3.2.1).We then study the higher-order contribution's effect on the fossil estimator by comparing the estimated long mode from GridSPT realizations with the result from the second-order Lagrangian Perturbation Theory (2LPT) realizations (Section 3.2.2) and fourth-order GridSPT (Section 3.2.3).Both GridSPT and 2LPT are fast enough to generate each realization in less than one minute.

The leading-order matter bispectrum from standard perturbation theory
Constructing the quadratic clustering-fossil estimator requires the squeeze-limit bispectrum to set the fossil kernel f (k 1 , k 2 ; k) in Equation (2.15).In SPT, the leading-order matter bispectrum from nonlinear gravity is given as with the second-order SPT kernel and the linear matter power spectrum P L (k).Note that the kernel F 2 (k 1 , k 2 ) only depends on the three wavenumbers (k 1 , k 2 , k 3 ) because of the statistical homogeneity and isotropy.The clustering-fossil estimator facilitates the squeezed limit of the bispectrum: and we can read off the following fossil-kernel expression from Equation (2.8) as: Note that the third term 2F 2 (k 1 , k 2 )P L (k 1 )P L (k 2 ) vanishes in the squeezed limit because F (q, −q) = 0.For a finite k, this term indeed contributes to the statistics of the recovered long mode.For example, when computing the cross power spectrum between the recovered long modes (using Equation (2.11)) and the original one, the term yields the correction to the squeezed-limit expression in Equation (2.19) as

Second-order GridSPT
The Grid-based standard perturbation theory (GridSPT) is a method of computing the n-th order SPT density and velocity fields from the recursion relation.Standard Perturbation Theory [63,[70][71][72][73][74][75][76] approximates the evolution as a pressure-less fluid with the following evolution equations for the density contrast δ(x, τ ) ≡ ρ m (x, τ )/ρ m (τ ) − 1 and the peculiar velocity v(x, τ ): along with the Poisson equation: Here, dot represents the conformal-time derivative, dτ = dt/a with a(t) being the scale factor and t being the cosmic time, ∇ is comoving-coordinate derivative, ρm is the mean matter density, and Φ is the peculiar gravitational potential.Note that we omit the spacetime coordinate in equations to avoid clutter.The set of equations describes the non-relativistic-matter (cold-dark matter and baryon) fluid on scales larger than the baryonic Jeans scale.Note that the nonlinearities in Equations(3.6)-(3.7)comes from the second-order terms such as ∇ • (δv) and (v • ∇)v.In SPT, we solve Equations(3.6)-(3.8)by expanding the nonlinear density contrast δ and velocitygradient field θ ≡ −∇ • v/(aHf ) as where D(τ ) is the linear growth factor.For a given realization of the linear density field on regular grid points, the GridSPT [53] provides a way to compute the matter density field δ and the velocity field v of LSS perturbatively by solving the fluid equations [Equations(3.6)-(3.8)],which becomes the recursion relation as: From a given set of linear density field, δ 1 = θ 1 = δ L , we can use Equation (3.10) to compute the nonlinear density δ (n) and velocity-gradient θ (n) fields order by order manner.Using the FFT (Fast Fourier Transform) to evaluate the ∇ operators on the right-hand side, the GridSPT enables us to quickly generate the nth order quantities.We use the second-order GridSPT, which contains first-and second-order density contrast, to test the fossil estimator.The squeezed bispectrum in the second-order GridSPT strictly follows the tree-level bispectrum in SPT so the constructed estimator can fully capture the coupling between long and short modes.Therefore, the second-order GridSPT is an ideal tool to conduct a proof-of-concept study of the fossil estimator without any uncontrolled systematic error.
We compute the GridSPT density fields on the N grid = 512 3 grids.We adopt the empirical cutoff k 1 = 1.0 h/Mpc for the first-order density contrast and k 2 = 1.33 h/Mpc for higher-order density contrasts used in [53].

Second-order Lagrangian perturbation theory (2LPT)
Lagrangian perturbation theory (LPT) is an alternative to the SPT in modeling the nonlinear density fields.The fundamental object in LPT is the displacement field Ψ(q, τ ) from the regular Lagrangian position q, which makes the Eulerian position x as x(q, τ ) = q + Ψ(q, τ ) . (3.11) We can obtain the LPT solutions by solving the equation of motion in an expanding universe ẍ + ȧ a ẋ = −∇ x Φ (3.12) for an irrotational displacement field ∇×Ψ = 0 perturbatively.The second-order solution for the displacement field is given as (see, for example, appendix E of [56] for a full derivation) where the linear Lagrangian potential is related to the linear density contrast as and the second-order Lagrangian potential is related to ϕ (1) as ,ii (q, τ )ϕ ,jj (q, τ ) − ϕ ,ij (q, τ ) (3.15) Here, D 2 (τ ) is the solution of the following differential equation: and D 2 = −3/7D 2 (τ ) for the Einstein de-Sitter (spatially flat and matter-dominated) universe.
The 2LPT (seond-order LPT) prescription is to displace regularly spaced particles using Equation (3.11) and Equation (3.13).Since the nonlinearities are modeled by particle displacement, the resulting density contrast, while agreeing with the SPT prediction to second order, contains a myriad of higher-order nonlinear contributions (see, e.g., Refs.[68] and [77]).The bispectrum of the 2LPT density field, therefore, deviates from the tree-level predictions in Equation (3.3), particularly on small scales.By applying the same fossil estimator using the kernel given in Equation (3.4) to the 2LPT density fields, we can test the behavior of the fossil estimator when ignoring the higher-order nonlinear couplings in data.
We implement the 2LPT by using the displacement of 512 3 Lagrangian particles in the cubic box of volume (1000 Mpc/h) 3 .We then estimate the density with N grid = 256 3 grids and preserve the density modes of wavenumber k < k Nyq /2 to avoid the aliasing effect [56].

Fourth-order GridSPT
While we test the effect of neglecting higher-order nonlinear couplings by comparing the second-order GridSPT and 2LPT, the density fields in these two toy cases are far from reality on small scales.For example, the nonlinearity in the density field is too strong in second-order GridSPT [78] while too weak in 2LPT [53] compared to N -body simulations.To overcome this, we also test the performance of the tree-level fossil estimator in a more realistic density field using the fourth-order GridSPT.The power spectrum and bispectrum in the fourth-order GridSPT can be accurately modeled by the complete one-loop power spectrum and one-loop bispectrum in SPT, which are closer to the N -body result than either second-order GridSPT or 2LPT [78].As for the q max (z), the smallest scale to be included in the reconstruction at each redshift, we use the result of Ref. [64], which has measured the maximum wavenumber below which the tree-level bispectrum accurately model the nonlinear bispectrum from a suite of N -body simulations.

Results
We present the results of the analysis in reconstructing the long modes by applying the tree-level fossil estimator to the following three nonlinear density fields: second-order and fourth-order GridSPT and 2LPT.

The tree-level fossil estimator on the 2LPT density field
First, we applied the tree-level fossil estimator to the 2LPT density fields to reconstruct the long modes.

The reconstructed long modes vs. the original modes
In Figure 1, we compare the reconstructed (left panel) and true (right panel) large-scale density field.We show the projected density distribution at z = 1 in the x-y plane of thickness 2 Mpc/h.We use short modes up to q max =0.4 h/Mpc for reconstruction, and 2LPT, L = 1000 Mpc/h, z = 0 The cross-correlation coefficient, r(k) defined in Equation (2.18), between the reconstructed and true long modes in 2LPT simulations at z = 0 (left) and z = 1 (right).Each line shows the result from different maximum wavenumber of short modes for reconstruction, q max = 0.1, 0.2, 0.3, and 0.4 h/Mpc.The error bars indicate the standard deviation of the mean measured from 100 2LPT realizations.
we smooth both fields with a Gaussian filter of radius R = 15 Mpc/h.In general, the reconstructed density field resembles the morphology of the original density field, but we can clearly see the visible differences between the two.This is most apparent for small-scale features in that some features in the original 2LPT field are missing in the reconstructed field, which implies that the recovered mode is more noise-dominated on smaller scales.
For a more quantitative comparison, we compute the cross-correlation coefficients between the recovered and original long modes, and the result is shown in Figure 2 for the results at redshifts z = 0 (left) and z = 1 (right).At each redshift, we show four different results corresponding to four different values of q max = 0.1, 0.2, 0.3, and 0.4 h/ Mpc.
First of all, the cross-correlation coefficient quickly drops as the long-mode wavenumber k approaches q max .That is because there are fewer and fewer short modes for reconstruction as k approaches q max .Then the recovered long modes are dominated by noise and correlate weakly with the original long modes.This is consistent with the lack of small-scale features we observe in the morphological comparison in Figure 1.
We also notice that the cross-correlation coefficients increase with q max , with an especially large improvement from q max = 0.1 h/Mpc to 0.2 h/Mpc.With q max = 0.4 h/Mpc, the cross-correlation coefficient reaches 0.95 on the large-scale end.This implies that the phase of the recovered long mode from the fossil estimator is highly correlated with the true long mode.Clearly, including more short modes leads to a more accurate reconstruction of the phase of the long modes.
Comparing the two panels, for a fixed q max , the cross-correlation coefficient at z = 0 is larger than z = 1 for all four q max cases.As derived in Equation (2.20), the cross-correlation coefficient is close to unity when P N G /P L ≪ 1.Since the noise power spectrum P N G on large scales only weakly depends on the redshift [Equation (2.16) and Equation (3.4)] while the amplitude of linear power spectrum grows in time, the cross-correlation coefficient is closer to unity at lower redshifts.In other words, a higher signal-to-noise ratio of the power spectrum leads to a larger cross-correlation coefficient at lower redshifts.This phenomenon can also be interpreted as a consequence of the stronger nonlinear coupling between the long and short For both panels, we show the linear matter power spectrum (signal) as a solid black line and the noise power spectrum with four different q max = 0.1, 0.2, 0.3, and 0.4 h/Mpc as different colors: measured from ten 2LPT realizations (dots with error bars), Gaussian estimate [Equation (2.16)] (dashed lines), and non-Gaussian estimate, full expressions in Equation (2.25) (solid lines).modes at lower redshift.That is, at higher redshifts, we can reach the same level of the cross-correlation coefficient only with an increasing dynamic range of the short modes used for the reconstruction.

Signal-to-noise ratio
To test the statistical significance of the reconstruction, we measure the noise power spectrum of the reconstructed long modes and show them in Figure 3 for two redshifts z = 0 (left) and z = 1 (right), and for four different maximum short mode wavenumber, q max .For each case, we compare the measured noise power spectrum (data points with errorbar) with the theoretical estimation from Gaussian assumption Equation (2.16) (dashed line) and the full computation [Equation (2.25)] including the trispectrum contribution (solid line).We find that while the full noise power spectrum estimate captures the measurement quite well, the Gaussian approximation always underestimates the noise.As coming from the non-Gaussian trispectrum, the discrepancy is most apparent from cases including more small-scale contributions or increasing q max .The signal-to-noise ratio of the long-mode reconstruction significantly improves as increasing q max , because more short-modes are used for the reconstruction.At the same time, including short modes in the nonlinear scales deviates the noise power spectrum from the Gaussian approximated one.Such a deviation does not merely cause an underestimation of the errorbar, because, as shown in Equation (2.17), estimating the long-mode power spectrum requires subtracting the noise power spectrum.That is, underestimation of the noise power spectrum leads to the systematic overestimation bias of the long-mode power spectrum.Section 4.1.3demonstrates this point explicitly.

The power spectrum of the recovered long modes
In Figure 4, we show the ensemble mean of the power spectra from the long modes, both auto and cross-correlation with the original mode, reconstructed from the 1,000 2LPT realizations at z = 1 (with q max = 0.1 h/Mpc) and compare them with the power spectrum P mm (k) of the original long mode.A perfect reconstruction would yield a perfect match.Figure 4.The ensemble mean of the power spectrum from 1,000 2LPT simulations at z = 1, using q max = 0.1 h/Mpc.We highlight the importance of the fossil kernel by comparing the fossil kernel f (k 1 , k 2 ; k) from the tree-level bispectrum (Left) and from the bispectrum and power spectrum measured from an average of 10,000 2LPT realizations (right).For both panels, black (sqaures), orange (diamonds), blue (triangles), and green (dots) are, respectively, the true matter power spectrum, the cross power spectrum between the recovered and the original field, and the recovered power spectrum subtracting the noise calculated from Gaussian approximation.We also show the ratio of each power spectrum to the P mm (k) in the bottom panels.
The left panel of Figure 4 shows the reconstruction with the tree-level fossil kernel f (k 1 , k 2 ; k) in Equation (3.4) coming from the tree-level bispectrum, and the right panel shows the results with the measured fossil kernel.We measured the fossil kernel by taking the ratio between the average of the bispectrum and the average long-mode power spectrum using the 10,000 2LPT simulations.
The cross-power spectrum P rm (k) from the tree-level fossil kernel (left panel) deviates from P mm (k), while that from the measured fossil kernel (right panel) is on top of the P mm (k).It is because the 2LPT bispectrum, even at q max = 0.1 h/Mpc at z = 1, deviates from the tree-level prediction.This result highlights the importance of having an accurate fossil kernel for the reconstruction.Equation (3.5) shows that the cross power spectrum P rm is determined by the ratio between the true bispectrum and the fossil kernel, which is equal to the linear matter power spectrum P L on large scales.The high-k deviation of P rm from P mm is due to the correlation between the second-order long mode and two first-order short modes (the second term in Equation (3.5)).
The auto power spectrum of the recovered modes P rr exceeds P mm due to the presence of the noise power spectrum (see Equation (2.17)).Therefore, P rr increase rapidly when k approaches q max .This is also consistent with the quick drop of the cross-correlation coefficient at high k as we showed in Figure 2. Subtracting the noise power spectrum P N G estimated using the Gaussian approximation reduces them closer to P mm (k), but still P rr −P N G disagrees with P mm .As we shall show in Section 4.2, this is due to the trispectrum contribution to the noise power spectrum.

The ideal toy: reconstruction from the second-order GridSPT
In Section 3.2.2,we show that using an accurate squeezed-limit bispectrum for the fossil kernel is essential to get the correct cross-correlation power spectrum.At the same time, the estimated auto-correlation of the recovered long mode deviates from that of the original field -16 - Figure 5.The reconstructed long mode power spectrum when applying the fossil estimator to the second-order GridSPT density field.The symbols are the same as Figure 4 except for the red pentagon, which shows the recovered power spectrum with the full noise subtracted.The recovered long-mode power spectrum matches the ground truth (black square) only after subtracting the noise power spectrum taking into account the non-Gaussian covariance.
when subtracting the noise power spectrum estimated using the Gaussian approximation.
To show the proof-of-concept case where we can recover both cross-correlation and auto power spectrum of the long mode, we apply the fossil estimator to the second-order GridSPT realizations.The second-order GridSPT realizations contain the nonlinear density field up to the second order so that we can estimate all relevant correlations.Figure 5 show the result of the reconstruction analysis using 1,000 second-order Grid-SPT realizations with q max = 0.1 h/Mpc.Just like in Section 3.2.2,when estimating the noise power spectrum with Gaussian approximation, the recovered power spectrum (the blue triangles) is about 30 percent higher than P mm .We reach the agreement (red pentagon symbols) only after including the full covariance matrix with the trispectrum contribution.To do so, we measure the off-diagonal covariance matrix numerically from 10,000 realizations of second-order GridSPT and calculate the full noise power spectrum P N full according to Equation (2.25).
Therefore, it is essential to model the correct noise power spectrum, including the offdiagonal terms in the covariance matrix or trispectrum.Neglecting these terms not only underestimates the uncertainties in the estimated long-mode power spectrum, but also introduces systematic bias in the power spectrum of recovered long modes.

A more realistic toy: reconstruction with the fourth-order GridSPT
Thus far, we have reconstructed the long modes from the controlled toy nonlinear density fields generated from the 2LPT (Section 3.2.2) and the second-order GridSPT (Section 4.2), from which we find that we have to use an accurate fossil kernel and to incorporate the non-Gaussian covariance (including the trispectrum) in order to achieve an unbiased reconstruction of the long modes.Although useful for the analysis, of course, none of the test nonlinear density fields is close to the real cosmic density field.In this section, we study the limitations of the tree-level fossil estimator by using the fourth-order GridSPT density field.4th GridSPT, L = 1000 Mpc/h z = 0, q max = 0.100 h/Mpc z = 1, q max = 0.133 h/Mpc z = 2, q max = 0.167 h/Mpc z = 3, q max = 0.200 h/Mpc Figure 6.The cross-correlation coefficient between the recovered field and the original field constructed using four-order GridSPT realizations.For the reconstruction, we use four different q max = 0.1, 0.133, 0.167, 0.2 h/Mpc, respectively, at z = 0, 1, 2, 3 such that the squeezed bispectrum works within 2% accuracy, as determined by [64].The error bars indicate the standard deviation of the mean measured from 100 GridSPT realizations.

The reach of the tree-level fossil estimator
We apply the tree-level fossil estimator in fourth-order GridSPT, using q max = 0.1, 0.133, 0.167, 0.2 h/Mpc, respectively, at z = 0, 1, 2, 3, which are the highest wavenumber where the tree-level bispectrum models the measured squeezed bispectrum in N -body simulations within the 2% accuracy [64].That is, the reconstruction must be unbiased within these dynamic ranges.
In Figure 6, we present the best cross-correlation coefficient that the estimator can achieve at z = 0, 1, 2, 3.For all four redshifts on large scales, the cross-correlation coefficients reach r = 0.7, which indicates the best signal-to-noise ratios of the recovered long mode based on tree-level bispectrum theory.
Figure 7 presents the measured auto-and cross-power spectra at the four redshifts.One notable feature is that the cross-power spectra match the true matter power spectra on large scales.This fact ensures that the tree-level fossil estimator is unbiased in the dynamical ranges.To enhance the signal-to-noise ratio of the reconstructed mode, we need to find a more accurate bispectrum model, e.g.one-loop SPT bispectrum, such that the estimator can safely include more short modes on smaller scales without biasing the recovered long mode.

The limit of the tree-level fossil estimator
In Section 2.2, we have proved that the fossil estimator is equivalent to the squeezed-limit bispectrum.Therefore, the fossil estimator reconstructs the linear-order density mode while the nonlinear part of the density mode is neglected.We test the scheme in this section by comparing the reconstructed power spectrum with the original density power spectrum which contains the nonlinear terms.
Figure 8 shows the difference between the reconstructed power spectrum and the original long-mode power spectrum in a unit of the cosmic variance uncertainty σ(P rr ).We use 100 fourth-order GridSPT realizations for this analysis.For all (z = 0, 1, 2, 3) redshifts.This plot shows that the nonlinear corrections are only a few percent of the cosmic variance at all four redshifts, which indicates the nonlinear correction of the recovered mode can be safely neglected with the range of q and the simulation box size we use in this study.Note that the nonlinear correction would become important if we work with either a bigger simulation box or a larger dynamical range of q, because the former raises the number of modes while the latter suppresses the noise power spectrum in P rr .In either case, the nonlinear correction could be comparable to the cosmic variance of P rr .And we have to carefully pick out the range of wavenumbers that the fossil estimator works.

Conclusion & Discussion
Here, we investigate the reconstruction of long-wavelength density modes based on the offdiagonal correlation between short-wavelength modes in the Fourier space, a phenomenon called "clustering fossil".The off-diagonal correlation can be understood as the local inhomogeneities introduced by the long mode.The coupling coefficient between short modes can be modeled from the squeezed-limit bispectrum and be used to pick up the specific nonlinear correlation that we use to estimate the long modes.
Throughout this paper, we use the tree-level bispectrum in standard perturbation theory to write down the quadratic estimator for the long mode in Equation (2.15), i.e., we only consider the coupling between the second-order short modes and the linear-order long mode.We have tested the same estimator in both second-order GridSPT and 2LPT simulations.By calculating the ensemble mean of the cross power spectrum P rm (k), we show that the estimator is unbiased in second-order GridSPT but biased in 2LPT.We manage to remove the bias of the recovered long mode in 2LPT simulations with the fossil kernel term measured from the 2LPT bispectrum.The results imply that the coupling coefficient from an accurate squeezed-limit bispectrum model guarantees an unbiased estimator.
We also aim to reconstruct the long-mode power spectrum.The auto power spectrum of the recovered mode P rr involves the four-point correlation function in Fourier space.P rr contains both the long mode power spectrum and the noise power spectrum.The noise power term is minimized if we use the inverse variance weight with the assumption that the short density modes are Gaussian, which is not true in the real universe.In fact, the noise power spectrum also receives the contribution from the connected four-point correlation function, which becomes more significant as we increase the largest wavenumber of short modes for reconstruction.Therefore, it is crucial to include the connected four-point correlation function into the noise power spectrum to accurately reconstruct the long-mode power spectrum.
We also demonstrate that the cross-correlation coefficient between the recovered and true long mode is determined by the ratio between the long-mode power spectrum and noise.To get higher cross-correlation coefficient, we can suppress the noise power spectrum by including more short modes for reconstruction.With the estimation from fourth-order GridSPT, we estimate the best cross-correlation coefficient the tree-level fossil estimator could achieve in N -body simulation is 0.7 at redshifts z = 0, 1, 2, 3 while keeping the recovered long mode unbiased.
In order to apply the fossil estimator in N -body simulations and enhance the reconstruction power, we need to include more short modes of larger wavenumbers for the reconstruction.Therefore, we need to use higher-order coupling terms beyond tree-level in the perturbation theory into the estimator.However, the estimator can still be biased because the perturbation theory fails on the fully nonlinear scale.Fortunately, the response approach [79,80] enables us to measure the squeezed-limit bispectrum in mild-sized N -body simula-tions accurately, which can potentially extend q max from quasi-linear scale to fully nonlinear scale without biasing the estimator.We leave this part to the future work.
To make this method more feasible in future galaxy surveys and 21cmIM experiments, we need to investigate the reconstruction from the biased tracer field.Using perturbative galaxy bias expansion, we can calculate the squeezed-limit galaxy-galaxy-matter bispectrum and reconstruct the matter density long mode from the galaxy short modes.We also need to include the redshift-space distortion effect and the foreground wedge in the 21cmIM [52].Our next step, then, is to understand how galaxy bias and the observational effects impact the reconstruction power.
There are many interesting applications of fossil estimators.By implementing the estimator in galaxy survey, the recovered matter long mode and observed galaxy long mode naturally form multi tracers can constrain the local-type primordial non-Gaussianity during the inflation without cosmic variance [81,82] (See [52] for more details.).Furthermore, we can construct the tensor-type fossil estimator to infer the primordial gravitational waves from the galaxy short modes, if the galaxy-galaxy-gravitational wave bispectrum can be modeled accurately.Based on the clustering fossil theory, this reconstruction method could become a powerful new statistical approach to probing the large-scale structure and increasing the scientific output of the future galaxy and 21cm surveys.

Figure 1 .
Figure1.The two-dimensional (the 2 h/Mpc slice on the x-y (z = 0) plane) morphology of the reconstructed density field (Left) and the 2LPT density fields (Right) at z = 1.We smooth all fields using the spherical Gaussian filter with the radius R = 15 h/Mpc.The maximum wavenumber of modes used for reconstruction is q max = 0.4 h/Mpc.As indicated by the color bar, both are normalized to their own variance.

Figure 3 .
Figure3.significance of the long-mode reconstruction at z = 0 (left panel) and z = 1(right panel).For both panels, we show the linear matter power spectrum (signal) as a solid black line and the noise power spectrum with four different q max = 0.1, 0.2, 0.3, and 0.4 h/Mpc as different colors: measured from ten 2LPT realizations (dots with error bars), Gaussian estimate [Equation (2.16)] (dashed lines), and non-Gaussian estimate, full expressions in Equation (2.25) (solid lines).