# Recovering the Missing Large-Scale Density Modes in 21cm Intensity Map from
the Scalar-Type Clustering Fossils

Zhenyuan Wang (王震远)  Donghui Jeong

###### Abstract

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

##  1 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 [22, 10]. Furthermore,
the frequency-dependence response of the instrument can lead to a foreground
wedge contamination in the Fourier space for
$k_{\parallel}<k_{\perp}\tan\psi$, where the angle $\psi$ 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].

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 [44]. In a statistically homogeneous universe, the two-point
correlation function
$\left<\delta^{*}({\bm{k}})\delta({\bm{k}}^{\prime})\right>$ in Fourier space
only has diagonal components proportional to
$\delta^{D}({\bm{k}}-{\bm{k}}^{\prime})$: 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. [44], we construct a quadratic clustering-fossil estimator for
the long modes based on the off-diagonal correlation between short modes.
Unlike Ref. [44], which pursues inflationary spectator fields using the non-
Gaussian coupling in the early universe [45, 46]; 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. [47, 48] 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) [49]. 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].

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) [50, 51, 52]. 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.

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 short-wavelength 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,

$\displaystyle f({\bm{k}})=\int{\rm
d}^{3}xf(\bm{x})e^{-i{\bm{k}}\cdot\bm{x}},$ (1.1) $\displaystyle
f(\bm{x})=\int\frac{{\rm
d}^{3}k}{(2\pi)^{3}}f({\bm{k}})e^{i{\bm{k}}\cdot\bm{x}}.$ (1.2)

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,

$\displaystyle{\bm{k}}_{1\cdots n}=\sum_{i=1}^{n}{\bm{k}}_{i}.$ (1.3)

##  2 Estimating long modes from the clustering-fossil estimator

We begin the section by motivating the clustering-fossil method [44] 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.

###  2.1 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 $\delta(\bm{x})\equiv n(\bm{x})/\bar{n}-1$, where $n(\bm{x})$ is the
density at location $\bm{x}$ and $\bar{n}$ is the cosmic mean density. The
two-point correlation function is defined in terms of the density contrast as

$\displaystyle\xi(\bm{r})=\left<\delta(\bm{x})\delta(\bm{x}+\bm{r})\right>\,.$
(2.1)

Here, the bracket $\left<...\right>$ represents the average over the
statistical ensemble of the cosmic density field. Note that $\xi(\bm{r})$
depends only on the separation vector $\bm{r}$ between two positions $\bm{x}$
and $\bm{x}+\bm{r}$ but is independent of the position $\bm{x}$ where we make
the measurement. This property manifests the statistical homogeneity.
Furthermore, the statistical isotropy would restrict $\xi(\bm{r})=\xi(r)$.
Equivalently, the power spectrum, the Fourier transformation of the two-point
correlation function, is defined as

$\displaystyle\left<\delta({\bm{k}}_{1})\delta({\bm{k}}_{2})\right>=(2\pi)^{3}\delta^{D}({\bm{k}}_{12})P({\bm{k}}_{1}),$
(2.2)

where the Dirac-delta operator, $\delta^{D}({\bm{k}}_{12})$ on the right-hand
side, encodes the statistical homogeneity. That is, the statistical
homogeneity dictates that the correlation between the Fourier modes
$\delta({\bm{k}}_{1})$ and $\delta({\bm{k}}_{2})$ vanishes unless
${\bm{k}}_{1}+{\bm{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 $\delta$
as follows:

$\displaystyle\zeta^{(n)}(\bm{r},\bm{s},\cdots)=\left<\delta(\bm{x})\delta(\bm{x}+\bm{r})\delta(\bm{x}+\bm{s})\cdots\right>$
(2.3)
$\displaystyle\left<\delta({\bm{k}}_{1})\cdots\delta({\bm{k}}_{n})\right>=(2\pi)^{3}\delta^{D}({\bm{k}}_{12\cdots
n})P^{(n)}({\bm{k}}_{1},\cdots,{\bm{k}}_{n})\,.$ (2.4)

We call a density field Gaussian when the two-point correlation function is
the only non-vanishing connected $n$-point correlation function [53];
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 ($\delta_{\ell}$) can be broken
down as the multiplication of the conditional correlation functions and the
conditioning long modes:

$\displaystyle\left\langle\delta_{\ell}(\bm{x}_{1})\cdots\delta_{\ell}(\bm{x}_{n})\delta(\bm{r})\delta(\bm{s})\cdots\right\rangle=\left\langle\left.\left\langle\delta(\bm{r})\delta(\bm{s})\cdots\right\rangle\right|_{\delta_{\ell}(\bm{x}_{1}^{\prime})\cdots\delta_{\ell}(\bm{x}_{n}^{\prime})}\delta_{\ell}(\bm{x}_{1})\cdots\delta_{\ell}(\bm{x}_{n})\right\rangle\,,$
(2.5)

where $\langle\rangle|_{X}$ stands for the correlator evaluated with the
condition $X$. The expression follows from the definition of conditional
probability $P(A\cap 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 [54, 55] or the
position-dependent bispectrum [56]. In this case, we may write the conditional
correlation function as

$\displaystyle\left.\left\langle\delta(\bm{r})\delta(\bm{s})\cdots\right\rangle\right|_{\delta_{\ell}(\bm{x}^{\prime})}=C+\left.\frac{{\rm
d}}{{\rm
d}\delta_{\ell}}\left\langle\delta(\bm{r})\delta(\bm{s})\cdots\right\rangle\right|_{\delta_{\ell}=0}\delta_{\ell}(\bm{x}^{\prime})+\cdots\,,$
(2.6)

where $C$ is a function independent from $\delta_{\ell}$, and we truncate the
expansion at the linear-order response. On large scales where the long mode
$\delta_{\ell}$ is in the linear regime, we compute the original correlation
function as

$\displaystyle\left\langle\delta_{\ell}(\bm{x})\delta(\bm{r})\delta(\bm{s})\cdots\right\rangle\simeq\left.\frac{{\rm
d}}{{\rm
d}\delta_{\ell}}\left\langle\delta(\bm{r})\delta(\bm{s})\cdots\right\rangle\right|_{\delta_{\ell}=0}\xi_{\ell}(\bm{x}-\bm{x}^{\prime})\,.$
(2.7)

The general idea of the clustering fossil estimator [44] is to use Equation
(2.6) as a starting point for estimating the long-wavelength mode
$\delta_{\ell}$ for the non-Gaussian density field with underlying
$(n+1)$-point correlation function taking the form of Equation (2.7).

As shown in [44], the same argument applies when using a long mode
$h_{\ell}^{s}$ with general spin $s$. The density field is a spin-0 case.

###  2.2 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

$\displaystyle\left<\delta_{\ell}({\bm{k}})\delta({\bm{k}}_{1})\delta({\bm{k}}_{2})\right>=(2\pi)^{3}\delta^{D}({\bm{k}}_{12}+{\bm{k}})f({\bm{k}}_{1},{\bm{k}}_{2};{\bm{k}})P_{\ell}(k),$
(2.8)

with kernel $f({\bm{k}}_{1},{\bm{k}}_{2};{\bm{k}})$ encoding the details of
the coupling between the long mode $\delta_{\ell}({\bm{k}})$ and two short
modes $\delta({\bm{k}}_{1})$ and $\delta({\bm{k}}_{2})$. Hereafter, we drop
the explicit subscript $\ell$ and implicitly use ${\bm{k}}\ll{\bm{k}}_{i}$ to
indicate the long mode. We emphasize that the three wave vectors in
$f({\bm{k}}_{1},{\bm{k}}_{2};{\bm{k}})$ must satisfy
${\bm{k}}_{12}+{\bm{k}}=0$.

As we have shown in Section 2.1, Equation (2.8) is equivalent to the following
conditional two-point correlator:

$\displaystyle\left<\delta({\bm{k}}_{1})\delta({\bm{k}}_{2})\right>\big{|}_{\delta({\bm{k}})}=f({\bm{k}}_{1},{\bm{k}}_{2};{\bm{k}})\delta^{*}({\bm{k}})\delta^{\rm
K}_{{\bm{k}}_{12},-{\bm{k}}}.$ (2.9)

Here, $\delta^{\rm K}$ refers to the Kronecker delta that we use instead of
Dirac delta, after absorbing appropriate dimensionful constants inside
$f({\bm{k}}_{1},{\bm{k}}_{2};{\bm{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 $\delta({\bm{k}})$. A simple estimator turning
Equation (2.9) around would be

$\displaystyle\hat{\delta}_{r}^{\rm
naive}({\bm{k}})=\left[\sum_{-{\bm{q}}_{i}}\frac{\delta(-{\bm{q}}_{i})\delta(-{\bm{k}}+{\bm{q}}_{i})}{f(-{\bm{q}}_{i},-{\bm{k}}+{\bm{q}}_{i};{\bm{k}})}\right]^{*}=\sum_{{\bm{q}}_{i}}\frac{\delta({\bm{q}}_{i})\delta({\bm{k}}-{\bm{q}}_{i})}{f({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i};-{\bm{k}})}\,.$
(2.10)

In the second equality, we use the reality of the density field
$\left[\delta({\bm{k}})\right]^{*}=\delta(-{\bm{k}})$ and that bispectra and
$f({\bm{k}}_{1},{\bm{k}}_{2};{\bm{k}})$ must be real for even-parity density
field:
$f^{*}({\bm{k}}_{1},{\bm{k}}_{2};{\bm{k}})=f(-{\bm{k}}_{1},-{\bm{k}}_{2};-{\bm{k}})$
[57]. For a given ${\bm{k}}$, this naive estimator sums over all possible
pairs $\delta({\bm{q}}_{i})$ and $\delta({\bm{k}}-{\bm{q}}_{i})$ weighted by a
known kernel function $f({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i};-{\bm{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<|{\bm{k}}-{\bm{q}}_{i}|\leq 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 [58, 59].

Noticing that the individual contribution in Equation (2.10) is subject to
different variances, we can introduce a normalized weight $W_{i}$

$\displaystyle\hat{\delta}_{r}({\bm{k}})=\sum_{i}W_{i}\frac{\delta({\bm{q}}_{i})\delta({\bm{k}}-{\bm{q}}_{i})}{f({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i};-{\bm{k}})},$
(2.11)

to find an optimal estimator which minimizes the variance of
$\hat{\delta}_{r}({\bm{k}})$:

$\displaystyle\quad\sigma^{2}\left[\hat{\delta}_{r}({\bm{k}})\right]=\sum_{ij}W_{i}W_{j}C_{ij}.$
(2.12)

First, we compute the covariance matrix $C_{ij}$ for individual contribution
in the naive estimator [Equation (2.10)]:

$\displaystyle
C_{ij}({\bm{q}}_{i},{\bm{q}}_{j};{\bm{k}})=V^{2}\frac{P(q_{i})P(|{\bm{k}}-{\bm{q}}_{i}|)}{f^{2}({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i};-{\bm{k}})}\delta^{\rm
K}_{{\bm{q}}_{j},\,-{\bm{q}}_{i}},$ (2.13)

by ignoring all the connected parts of four-point correlators. Here $P(q)$ and
$P(|{\bm{k}}-{\bm{q}}|)$ are the measured nonlinear power spectra of the short
modes. In this case, the inverse-variance weight

$\displaystyle
W_{i}=\sum_{j}C_{ij}^{-1}\big{/}\sum_{ij}C_{ij}^{-1},\;\sigma^{2}\left[\hat{\delta}_{r}({\bm{k}})\right]=1\big{/}\sum_{j}C_{ij}^{-1}.$
(2.14)

gives the minimum variance estimator. Using the inverse-variance weight in
Equation (2.14), the optimal fossil estimator for the long mode becomes

$\displaystyle\hat{\delta}_{r}({\bm{k}})=P^{N}_{G}(k)\sum_{{\bm{q}}}\frac{f({\bm{q}},{\bm{k}}-{\bm{q}};-{\bm{k}})}{VP({\bm{q}})P({\bm{k}}-{\bm{q}})}\delta({\bm{q}})\delta({\bm{k}}-{\bm{q}}),$
(2.15)

the variance of the optimal estimator $\hat{\delta}_{r}$ becomes

$\displaystyle
P^{N}_{G}(k)=\left[\sum_{{\bm{q}}}\frac{\left|f({\bm{q}},{\bm{k}}-{\bm{q}};-{\bm{k}})\right|^{2}}{VP({\bm{q}})P({\bm{k}}-{\bm{q}})}\right]^{-1}\,.$
(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 $\delta$:

$\displaystyle{P}_{rr}(k)\equiv\frac{1}{V}\frac{1}{N_{k}}\sum_{\hat{\bm{k}}}\left<|\hat{\delta}_{r}({\bm{k}})|^{2}\right>=P_{L}(k)+P^{N}_{G}(k)\,.$
(2.17)

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_{G}^{N}(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 $\hat{\delta}_{r}({\bm{k}})$ and the original one
$\delta({\bm{k}})$ that we designate with subscript $m$:

$\displaystyle r(k)\equiv\frac{P_{rm}(k)}{\sqrt{{P}_{rr}(k)P_{mm}(k)}}$ (2.18)

Here,
$\left<\hat{\delta}_{r}({\bm{k}})\delta_{m}({\bm{k}}^{\prime})\right>=(2\pi)^{3}P_{rm}(k)\delta^{D}({\bm{k}}+{\bm{k}}^{\prime})$
is the cross power spectrum

$\displaystyle
P_{rm}(k)\equiv\frac{1}{V}\frac{1}{N_{k}}\sum_{\hat{\bm{k}}}{\rm
Re}\left<\hat{\delta}_{r}({\bm{k}})\delta^{*}({\bm{k}})\right>=\sum_{{\bm{q}}_{i}}W_{i}\frac{B({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i},-{\bm{k}})}{f({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i};-{\bm{k}})}\simeq
P_{L}(k)\,,$ (2.19)

which approaches the linear power spectrum for the squeezed ($k\to 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\to 0$ and $P_{mm}\to
P_{L}(k)$] as,

$\displaystyle r(k)\simeq\frac{1}{\sqrt{1+P^{N}_{G}(k)/P_{L}(k)}}.$ (2.20)

From Equation (2.20), we see that $r(k)\to 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 $|{\bm{k}}-{\bm{q}}|$
because most reconstruction power comes from the high-frequency
($|{\bm{q}}|\gg k$) modes. In this paper, we use the short modes
$\delta({\bm{q}}),\,\delta({\bm{k}}-{\bm{q}})$ with wavenumbers satisfying

$\displaystyle k<|{\bm{k}}-{\bm{q}}|\leq q\leq q_{\rm max},$ (2.21)

where $q_{\rm max}$ is the wavenumber of the smallest scale used for
reconstruction. Ideally, we want to make $q_{\rm 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_{\rm max}$ in the
quasi-linear scale where perturbation theory accurately models the
nonlinearities [60, 61]. 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, ${\bm{q}}^{\prime}$=${\bm{k}}-{\bm{q}}$ is identical to
${\bm{q}}$), which leads to a factor of two difference between Equation (2.16)
here and the expressions shown in previous studies [44, 47, 48, 62]. 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, all 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].
We now present the full noise power spectrum, including the connected four-
point function.

###  2.3 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

$\displaystyle C_{ij}({\bm{q}}_{i},{\bm{q}}_{j};{\bm{k}})$ $\displaystyle=$
$\displaystyle\left<\frac{\delta({\bm{q}}_{i})\delta({\bm{k}}-{\bm{q}}_{i})}{f({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i};-{\bm{k}})}\frac{\delta({\bm{q}}_{j})\delta(-{\bm{k}}-{\bm{q}}_{j})}{f({\bm{q}}_{j},-{\bm{k}}-{\bm{q}}_{j};{\bm{k}})}\right>\bigg{|}_{\delta({\bm{k}}),{\delta(-{\bm{k}})}}$
(2.22)
$\displaystyle-\left<\frac{\delta({\bm{q}}_{i})\delta({\bm{k}}-{\bm{q}}_{i})}{f({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i};-{\bm{k}})}\right>\bigg{|}_{\delta({\bm{k}})}\left<\frac{\delta({\bm{q}}_{j})\delta(-{\bm{k}}-{\bm{q}}_{j})}{f({\bm{q}}_{j},-{\bm{k}}-{\bm{q}}_{j};{\bm{k}})}\right>\bigg{|}_{\delta(-{\bm{k}})}$
$\displaystyle=$ $\displaystyle
V^{2}\frac{P(q_{i})P(|{\bm{k}}-{\bm{q}}_{i}|)}{|f({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i};-{\bm{k}})|^{2}}\left(\delta^{\rm
K}_{{\bm{q}}_{j},\,-{\bm{q}}_{i}}+\delta^{\rm
K}_{{\bm{q}}_{j},\,-({\bm{k}}-{\bm{q}}_{i})}\right)$
$\displaystyle+V\frac{T\left({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i},{\bm{q}}_{j},-{\bm{k}}-{\bm{q}}_{j}\right)}{f({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i};-{\bm{k}})f({\bm{q}}_{j},-{\bm{k}}-{\bm{q}}_{j};{\bm{k}})}\,,$

where
$T\left({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i},{\bm{q}}_{j},-{\bm{k}}-{\bm{q}}_{j}\right)$
is the connected four-point function, or trispectrum, defined as follows:

$\displaystyle
T\left({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i},{\bm{q}}_{j},-{\bm{k}}-{\bm{q}}_{j}\right)=$
$\displaystyle\left<\delta({\bm{q}}_{i})\delta({\bm{k}}-{\bm{q}}_{i})\delta({\bm{q}}_{j})\delta(-{\bm{k}}-{\bm{q}}_{j})\right>$
$\displaystyle-\left<\delta({\bm{q}}_{i})\delta({\bm{k}}-{\bm{q}}_{i})\right>\left<\delta({\bm{q}}_{j})\delta(-{\bm{k}}-{\bm{q}}_{j})\right>$
$\displaystyle-\left<\delta({\bm{q}}_{i})\delta({\bm{q}}_{j})\right>\left<\delta({\bm{k}}-{\bm{q}}_{i})\delta(-{\bm{k}}-{\bm{q}}_{j})\right>$
$\displaystyle-\left<\delta({\bm{q}}_{i})\delta(-{\bm{k}}-{\bm{q}}_{j})\right>\left<\delta({\bm{k}}-{\bm{q}}_{i})\delta({\bm{q}}_{j})\right>\,.$
(2.23)

After removing the duplicated contributions by imposing the inequality in
Equation (2.21), the covariance matrix becomes

$\displaystyle
C_{ij}({\bm{q}}_{i},{\bm{q}}_{j};{\bm{k}})=V^{2}\frac{P(q_{i})P(|{\bm{k}}-{\bm{q}}_{i}|)}{|f({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i};-{\bm{k}})|^{2}}\delta^{\rm
K}_{{\bm{q}}_{j},\,-{\bm{q}}_{i}}+V\frac{T\left({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i},{\bm{q}}_{j},-{\bm{k}}-{\bm{q}}_{j}\right)}{f({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i};-{\bm{k}})f({\bm{q}}_{j},-{\bm{k}}-{\bm{q}}_{j};{\bm{k}})}\,.$

From the covariance matrix, we compute the full noise power spectrum for the
optimal quadratic estimator, Equation (2.15):

$\displaystyle P_{\rm
full}^{N}(k)=P^{N}_{G}(k)+\bigg{\\{}\left[P^{N}_{G}(k)\right]^{2}\sum_{ij}\frac{T\left({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i},{\bm{q}}_{j},-{\bm{k}}-{\bm{q}}_{j}\right)f({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i})f({\bm{q}}_{j},-{\bm{k}}-{\bm{q}}_{j})}{V^{2}P(q_{i})P(q_{j})P(|{\bm{k}}-{\bm{q}}_{i}|)P(|-{\bm{k}}-{\bm{q}}_{j}|)}\bigg{\\}}\,.$
(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 [63, 64], we expect that the noise
power spectrum $P^{N}_{\rm full}(k)$ would saturate beyond some $q_{\rm max}$.
Later, we verify the statements by comparing numerically the two noise power
spectra $P^{N}_{\rm full}(k)$ and $P^{N}_{G}(k)$ for the two nonlinear density
realizations: 2LPT [65] and GridSPT [50].

##  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) [49]. 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 [47, 48, 62], which tested
the quadratic estimator based upon the second-order SPT (Equation (3.4) below)
against a few $N$-body simulation results. 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_{\rm
max}$, because (A) the squeezed-limit bispectrum deviates from the second-
order SPT prediction, and (B) the trispectrum contribution from small scales
to noise power spectrum becomes important. We cannot distinguish the 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.

We have performed the comparison study at two redshifts ($z=1$ and $z=0$). For
all studies, we use the simulation box of $V=(1000~{}{\rm Mpc}/h)^{3}$, and we
adopt the following cosmological parameters [24]: $\Omega_{b0}h^{2}=0.022307$,
$\Omega_{c0}h^{2}=0.11865$, $h=0.6778$, $n_{s}=0.9672$, $\Omega_{\nu
0}h^{2}=0.000638$, $\Omega_{\Lambda 0}=0.69179$, ${\cal A}_{s}=2.147\times
10^{-9}$ and $\sigma_{8}=0.8166$. We use CAMB [66] to generate the input
linear power spectrum of GridSPT and 2LPT.

###  3.1 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({\bm{k}}_{1},{\bm{k}}_{2};{\bm{k}})$ in Equation (2.15). In SPT, the
leading-order matter bispectrum from nonlinear gravity is given as

$\displaystyle
B(k_{1},k_{2},k_{3})=2F_{2}({\bm{k}}_{1},{\bm{k}}_{2})P_{L}(k_{1})P_{L}(k_{2})+(2~{}{\rm
cyclic})\,,$ (3.1)

with the second-order SPT kernel

$\displaystyle
F_{2}({\bm{k}}_{1},{\bm{k}}_{2})=\frac{7}{5}+\frac{2}{7}\left(\frac{{\bm{k}}_{1}\cdot{\bm{k}}_{2}}{k_{1}k_{2}}\right)^{2}+\frac{1}{2}\frac{{\bm{k}}_{1}\cdot{\bm{k}}_{2}}{k_{1}k_{2}}\left(\frac{k_{1}}{k_{2}}+\frac{k_{2}}{k_{1}}\right)$
(3.2)

and the linear matter power spectrum $P_{L}(k)$. Note that the kernel
$F_{2}({\bm{k}}_{1},{\bm{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:

$\displaystyle B(k,k_{1},k_{2})\xrightarrow{k\to
0}\left[2F_{2}({\bm{k}}_{1},{\bm{k}})P_{L}(k_{1})+2F_{2}({\bm{k}}_{2},{\bm{k}})P_{L}(k_{2})\right]P_{L}(k)\,,$
(3.3)

and we can read off the following fossil-kernel expression from Equation (2.8)
as:

$\displaystyle
f({\bm{k}}_{1},{\bm{k}}_{2};{\bm{k}})=2F_{2}({\bm{k}}_{1},{\bm{k}})P_{L}(k_{1})+2F_{2}({\bm{k}}_{2},{\bm{k}})P_{L}(k_{2}).$
(3.4)

Note that the third term
$2F_{2}({\bm{k}}_{1},{\bm{k}}_{2})P_{L}(k_{1})P_{L}(k_{2})$ vanishes in the
squeezed limit because $F({\bm{q}},-{\bm{q}})=0$. For a finite ${\bm{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

$\displaystyle P_{rm}(k)$ $\displaystyle=$
$\displaystyle\sum_{i}W_{i}\frac{B({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i},-{\bm{k}})}{f({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i};-{\bm{k}})}$
(3.5) $\displaystyle=$ $\displaystyle
P_{L}(k)+2P^{N}(k)\sum_{i}\frac{F_{2}({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i})P_{L}(q_{i})P_{L}(|{\bm{k}}-{\bm{q}}_{i}|)}{f({\bm{q}}_{i},{\bm{k}}-{\bm{q}}_{i};-{\bm{k}})}\,.$

###  3.2 Nonlinear density fields

####  3.2.1 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 [67, 68, 69, 70, 71, 72, 73, 60] approximates the
evolution as a pressure-less fluid with the following evolution equations for
the density contrast
$\delta({\bm{x}},\tau)\equiv\rho_{m}({\bm{x}},\tau)/\bar{\rho}_{m}(\tau)-1$
and the peculiar velocity ${\bm{v}}({\bm{x}},\tau)$:

$\displaystyle\dot{\delta}+\nabla\cdot\left[(1+\delta){\bm{v}}\right]=0\,,$
(3.6)
$\displaystyle\dot{\bm{v}}+({\bm{v}}\cdot\nabla){\bm{v}}+\frac{\dot{a}}{a}{\bm{v}}=-\nabla\Phi\,,$
(3.7)

along with the Poisson equation:

$\nabla^{2}\Phi=4\pi G\bar{\rho}_{\mathrm{m}}\,a^{2}\delta\,.$ (3.8)

Here, dot represents the conformal-time derivative, $d\tau=dt/a$ with $a(t)$
being the scale factor and $t$ being the cosmic time, $\nabla$ is comoving-
coordinate derivative, $\bar{\rho}_{\mathrm{m}}$ is the mean matter density,
and $\Phi$ 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 $\nabla\cdot(\delta{\bm{v}})$ and $({\bm{v}}\cdot\nabla){\bm{v}}$. In
SPT, we solve Equations(3.6)–(3.8) by expanding the nonlinear density contrast
$\delta$ and velocity-gradient field $\theta\equiv-\nabla\cdot{\bm{v}}/(aHf)$
as

$\delta({\bm{x}},\tau)=\sum_{n}[D(\tau)]^{n}\delta^{(n)}({\bm{x}})\,,~{}~{}\theta({\bm{x}},\tau)=\sum_{n}[D(\tau)]^{n}\theta^{(n)}({\bm{x}})\,,$
(3.9)

where $D(\tau)$ is the linear growth factor.

For a given realization of the linear density field on regular grid points,
the GridSPT [50] provides a way to compute the matter density field $\delta$
and the velocity field ${\bm{v}}$ of LSS perturbatively by solving the fluid
equations [Equations(3.6)–(3.8)], which becomes the recursion relation as:

$\displaystyle\left(\begin{array}[]{c}\delta^{(n)}(\bm{x})\\\
\theta^{(n)}(\bm{x})\end{array}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{(2n+3)(n-1)}\left(\begin{array}[]{ll}2n+1{}{}&1\\\
3{}{}&n\end{array}\right)\sum_{m=1}^{n-1}\left(\begin{array}[]{c}\nabla\cdot\left(\delta^{(m)}{\bm{v}}^{(n-m)}\right)\\\
\nabla^{2}\left({\bm{v}}^{(n-m)}\cdot{\bm{v}}^{(m)}\right)\end{array}\right)\,.$
(3.16)

From a given set of linear density field, $\delta_{1}=\theta_{1}=\delta_{L}$,
we can use Equation (3.16) to compute the nonlinear density $\delta^{(n)}$ and
velocity-gradient $\theta^{(n)}$ fields order by order manner. Using the FFT
(Fast Fourier Transform) to evaluate the $\nabla$ operators on the right-hand
side, the GridSPT enables us to quickly generate the $n$th 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_{\rm 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
[50].

####  3.2.2 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 $\bm{\Psi}({\bm{q}},\tau)$ from the regular Lagrangian
position ${\bm{q}}$, which makes the Eulerian position $\bm{x}$ as

$\bm{x}({\bm{q}},\tau)={\bm{q}}+\bm{\Psi}({\bm{q}},\tau)\,.$ (3.17)

We can obtain the LPT solutions by solving the equation of motion in an
expanding universe

$\ddot{\bm{x}}+\frac{\dot{a}}{a}\dot{\bm{x}}=-\nabla_{\bm{x}}\Phi\,$ (3.18)

for an irrotational displacement field $\nabla\times\bm{\Psi}=0$
perturbatively. The second-order solution for the displacement field is given
as (see, for example, appendix E of [53] for a full derivation)

$\bm{\Psi}({\bm{q}},\tau)=-\nabla_{{\bm{q}}}\phi^{(1)}({\bm{q}},\tau)+\nabla_{{\bm{q}}}\phi^{(2)}({\bm{q}},\tau)\,,$
(3.19)

where the linear Lagrangian potential is related to the linear density
contrast as

$\nabla_{{\bm{q}}}^{2}\phi^{(1)}({\bm{q}},\tau)=\delta^{(1)}(\bm{x},\tau)\,,$
(3.20)

and the second-order Lagrangian potential is related to $\phi^{(1)}$ as

$\nabla_{{\bm{q}}}^{2}\phi^{(2)}({\bm{q}},\tau)=-\frac{D_{2}(\tau)}{D(\tau)^{2}}\sum_{i>j}\left\\{\phi^{(1)}_{,ii}({\bm{q}},\tau)\phi^{(1)}_{,jj}({\bm{q}},\tau)-\left[\phi^{(1)}_{,ij}({\bm{q}},\tau)\right]^{2}\right\\}\,.$
(3.21)

Here, $D_{2}(\tau)$ is the solution of the following differential equation:

$\ddot{D}_{2}(\tau)+\frac{\dot{a}}{a}\dot{D}_{2}(\tau)-\frac{3}{2}\left(\frac{\dot{a}}{a}\right)^{2}\left[\Omega_{m}(\tau)D_{2}(\tau)-D^{2}(\tau)\right]=0\,,$
(3.22)

and $D_{2}=-3/7D^{2}(\tau)$ 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.17) and Equation (3.19). 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. [65] and [74]). 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_{\rm grid}=256^{3}$ grids and preserve the density modes of
wavenumber $k<k_{\rm Nyq}/2$ to avoid the aliasing effect [53].

####  3.2.3 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 [75] while too weak
in 2LPT [50] 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 [75]. As for the $q_{\rm
max}(z)$, the smallest scale to be included in the reconstruction at each
redshift, we use the result of Ref. [61], which has measured the maximum
wavenumber below which the tree-level bispectrum accurately model the
nonlinear bispectrum from a suite of $N$-body simulations.

##  4 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.

###  4.1 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.

####  4.1.1 The reconstructed long modes vs. the original modes

Figure 1: 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_{{\rm max}}=0.4\,h/$Mpc. As indicated by the
color bar, both are normalized to their own variance.

Figure 2: 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_{\rm 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.

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_{\rm
max}=$0.4 $h/$Mpc for reconstruction, and 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_{\rm 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_{\rm max}$. That is because there are fewer and
fewer short modes for reconstruction as $k$ approaches $q_{\rm 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_{\rm
max}$, with an especially large improvement from $q_{\rm max}=0.1\;h/$Mpc to
0.2 $h/$Mpc. With $q_{\rm 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_{\rm max}$, the cross-correlation
coefficient at $z=0$ is larger than $z=1$ for all four $q_{\rm max}$ cases. As
derived in Equation (2.20), the cross-correlation coefficient is close to
unity when $P_{G}^{N}/P_{L}\ll 1$. Since the noise power spectrum $P_{G}^{N}$
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 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.

####  4.1.2 Signal-to-noise ratio

Figure 3: The 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_{\rm 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).

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_{\rm 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_{\rm max}$.

The signal-to-noise ratio of the long-mode reconstruction significantly
improves as increasing $q_{\rm 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.3 demonstrates this point explicitly.

####  4.1.3 The power spectrum of the recovered long modes

Figure 4: The ensemble mean of the power spectrum from 1,000 2LPT simulations
at $z=1$, using $q_{{\rm max}}=0.1\,h$/Mpc. We highlight the importance of the
fossil kernel by comparing the fossil kernel
$f({\bm{k}}_{1},{\bm{k}}_{2};{\bm{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.

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_{\rm 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.

The left panel of Figure 4 shows the reconstruction with the tree-level fossil
kernel $f({\bm{k}}_{1},{\bm{k}}_{2};{\bm{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_{\rm max}=0.1~{}h/{\rm 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_{\rm 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_{G}^{N}$
estimated using the Gaussian approximation reduces them closer to $P_{mm}(k)$,
but still $P_{rr}-P_{G}^{N}$ disagrees with $P_{mm}$. As we shall show in
Section 4.2, this is due to the trispectrum contribution to the noise power
spectrum.

###  4.2 The ideal toy: reconstruction from the second-order GridSPT

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.

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 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 GridSPT realizations with $q_{{\rm 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}_{\rm full}$ according to Equation (2.25).

Therefore, it is essential to model the correct noise power spectrum,
including the off-diagonal 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.

###  4.3 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.

####  4.3.1 The reach of the tree-level fossil estimator

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_{{\rm
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
[61]. The error bars indicate the standard deviation of the mean measured from
100 GridSPT realizations.

Figure 7: The ensemble mean of the power spectrum measured from 1,000 GridSPT
realizations at $z=0$, $1$, $2$, and $3$, using the $q_{{\rm max}}$ determined
as one-loop bispectrum models the squeezed-limit with 2 % accuracy [61]. The
blue triangles are the recovered power spectrum subtracting the noise
calculated from the Gaussian approximation. The orange diamonds are the cross-
power spectrum between the recovered and the original field. Symbols are
identical to Figure 4. For all cases, the bottom panel shows the ratio between
the power spectra and the true matter power spectrum $P_{mm}$.

We apply the tree-level fossil estimator in fourth-order GridSPT, using
$q_{\rm 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 [61]. 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.

####  4.3.2 The limit of the tree-level fossil estimator

Figure 8: The ratio between the neglected nonlinear correction to the
reconstructed matter power spectrum and the cosmic variance of the recovered
auto power spectrum at $z=0,1,2,3$. The results are measured in 100 fourth-
order GridSPT realizations. We use the same $q_{\rm max}$ values as in Figure
6.

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
$\sigma(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.

##  5 Conclusion & Discussion

Here, we investigate the reconstruction of long-wavelength density modes based
on the off-diagonal 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. As a result, the noise power
spectrum would receive the contribution from the connected four-point
correlation function. Furthermore, the contributions from connected terms
would become increasingly important if we increased the largest wavenumber of
short modes for reconstruction.

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 a 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
[76, 77] enables us to measure the squeezed-limit bispectrum in mild-sized
$N$-body simulations accurately, which can potentially extend $q_{\rm 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. 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 and can constrain the local-type
primordial non-Gaussianity during the inflation without cosmic variance [78,
79]. 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.

##  Appendix A Fossil Estimator in the Continuous Limit

To get a fast estimation of the noise power spectrum of the reconstructed
modes, we can write down the fossil estimator in the continuous limit.

$\displaystyle\delta_{r}(k)$ $\displaystyle=$ $\displaystyle
P^{N}_{G}({\bm{k}})\int_{q}\frac{\delta({\bm{q}})\delta({\bm{k}}-{\bm{q}})f({\bm{q}},{\bm{k}}-{\bm{q}})}{P(q)P(|{\bm{k}}-{\bm{q}}|)}$
$\displaystyle=$ $\displaystyle P^{N}_{G}(k)\int_{k}^{q_{\rm max}}\frac{{\rm
d}q}{(2\pi)^{2}}q^{2}\int{\rm
d}\mu\frac{\delta({\bm{q}})\delta({\bm{k}}-{\bm{q}})f({\bm{q}},{\bm{k}}-{\bm{q}})}{P(q)P(|{\bm{k}}-{\bm{q}}|)}.$

In the continuous limit, the noise power spectrum under the Gaussian
assumption is

$\displaystyle P^{N}_{G}(k)$ $\displaystyle=$
$\displaystyle\left[\int_{k}^{q_{\rm max}}\frac{{\rm
d}q}{(2\pi)^{2}}q^{2}\int{\rm
d}\mu\frac{\left|f({\bm{q}},{\bm{k}}-{\bm{q}})\right|^{2}}{P(q)P(|{\bm{k}}-{\bm{q}}|)}\right]^{-1}$
(A.2)

Here we have used the relation between the integration and the grid-based
discrete summation over the Fourier space

$\displaystyle\frac{1}{V}\sum=\int\frac{{\rm d}^{3}q}{(2\pi)^{3}}.$ (A.3)

The range of $\mu\equiv\hat{{\bm{k}}}\cdot\hat{{\bm{q}}}$ is constrained by
$q\geq|{\bm{k}}-{\bm{q}}|>k$, which is

$\displaystyle\frac{k}{2q}\leq\mu\leq{\rm min}\\{1,\,\frac{q}{2k}\\}.$ (A.4)

## Acknowledgments

We want to thank Caryl Gronwall and Michael Sigel for useful comments on the
early draft of this paper. ZW wants to acknowledge Hong-ming Zhu and Oliver
Philcox for useful discussions and Yuanheng Wang for helpful support during
the completion of this paper. This research made use of the Roar Supercomputer
at Penn State University. ZW and DJ acknowledge support from the National
Science Foundation Grant No. AST-2307026 at Penn State University. DJ is
supported by KIAS Individual Grant PG088301 at Korea Institute for Advanced
Study and was supported by NASA ATP program (80NSSC18K1103) at Penn State
University.

## References

  * [1] T.-C. Chang, U.-L. Pen, J.B. Peterson and P. McDonald, _Baryon Acoustic Oscillation Intensity Mapping of Dark Energy_ , _Phys. Rev. Lett._ 100 (2008) 091303 [0709.3672].
  * [2] A. Loeb and J.S.B. Wyithe, _Possibility of Precise Measurement of the Cosmological Power Spectrum with a Dedicated Survey of 21cm Emission after Reionization_ , _Phys. Rev. Lett._ 100 (2008) 161301 [0801.1677].
  * [3] X. Chen, _The Tianlai Project: a 21CM Cosmology Experiment_ , in _International Journal of Modern Physics Conference Series_ , vol. 12 of _International Journal of Modern Physics Conference Series_ , pp. 256–263, March, 2012, DOI [1212.6278].
  * [4] K. Bandura, G.E. Addison, M. Amiri, J.R. Bond, D. Campbell-Wilson, L. Connor et al., _Canadian Hydrogen Intensity Mapping Experiment (CHIME) pathfinder_ , in _Ground-based and Airborne Telescopes V_ , L.M. Stepp, R. Gilmozzi and H.J. Hall, eds., vol. 9145 of _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , p. 914522, July, 2014, DOI [1406.2288].
  * [5] L.B. Newburgh, K. Bandura, M.A. Bucher, T.C. Chang, H.C. Chiang, J.F. Cliche et al., _HIRAX: a probe of dark energy and radio transients_ , in _Ground-based and Airborne Telescopes VI_ , H.J. Hall, R. Gilmozzi and H.K. Marshall, eds., vol. 9906 of _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , p. 99065X, August, 2016, DOI [1607.02059].
  * [6] R.A. Battye, I.W.A. Browne, C. Dickinson, G. Heron, B. Maffei and A. Pourtsidou, _H I intensity mapping: a single dish approach_ , _Mon. Not. R. Astron. Soc._ 434 (2013) 1239 [1209.0343].
  * [7] J. Pritchard, K. Ichiki, A. Mesinger, R.B. Metcalf, A. Pourtsidou, M. Santos et al., _Cosmology from EoR/Cosmic Dawn with the SKA_ , in _Advancing Astrophysics with the Square Kilometre Array (AASKA14)_ , p. 12, April, 2015, DOI [1501.04291].
  * [8] H.-J. Seo, S. Dodelson, J. Marriner, D. Mcginnis, A. Stebbins, C. Stoughton et al., _A Ground-based 21 cm Baryon Acoustic Oscillation Survey_ , _Astrophys. J._ 721 (2010) 164 [0910.5007].
  * [9] E. Chapman and V. Jelić, _Foregrounds and their mitigation_ , _arXiv e-prints_ (2019) arXiv:1909.12369 [1909.12369].
  * [10] A. Liu and J.R. Shaw, _Data Analysis for Precision 21 cm Cosmology_ , _PASP_ 132 (2020) 062001 [1907.08211].
  * [11] M.F. Morales, J.D. Bowman and J.N. Hewitt, _Improving Foreground Subtraction in Statistical Observations of 21 cm Emission from the Epoch of Reionization_ , _Astrophys. J._ 648 (2006) 767 [astro-ph/0510027].
  * [12] X. Wang, M. Tegmark, M.G. Santos and L. Knox, _21 cm Tomography with Foregrounds_ , _Astrophys. J._ 650 (2006) 529 [astro-ph/0501081].
  * [13] J.D. Bowman, M.F. Morales and J.N. Hewitt, _Foreground Contamination in Interferometric Measurements of the Redshifted 21 cm Power Spectrum_ , _Astrophys. J._ 695 (2009) 183 [0807.3956].
  * [14] A. Liu, M. Tegmark and M. Zaldarriaga, _Will point sources spoil 21-cm tomography?_ , _Mon. Not. R. Astron. Soc._ 394 (2009) 1575 [0807.3952].
  * [15] A. Liu and M. Tegmark, _A method for 21 cm power spectrum estimation in the presence of foregrounds_ , _Phys. Rev. D_ 83 (2011) 103006 [1103.0281].
  * [16] A. Parsons, J. Pober, M. McQuinn, D. Jacobs and J. Aguirre, _A Sensitivity and Array-configuration Study for Measuring the Power Spectrum of 21 cm Emission from Reionization_ , _Astrophys. J._ 753 (2012) 81 [1103.2135].
  * [17] E. Chapman, F.B. Abdalla, G. Harker, V. Jelić, P. Labropoulos, S. Zaroubi et al., _Foreground removal using FASTICA: a showcase of LOFAR-EoR_ , _Mon. Not. R. Astron. Soc._ 423 (2012) 2518 [1201.2190].
  * [18] E. Chapman, F.B. Abdalla, J. Bobin, J.L. Starck, G. Harker, V. Jelić et al., _The scale of the problem: recovering images of reionization with Generalized Morphological Component Analysis_ , _Mon. Not. R. Astron. Soc._ 429 (2013) 165 [1209.4769].
  * [19] J.S. Dillon, A. Liu and M. Tegmark, _A fast method for power spectrum and foreground analysis for 21 cm cosmology_ , _Phys. Rev. D_ 87 (2013) 043005 [1211.2232].
  * [20] L. Wolz, C. Blake, F.B. Abdalla, C.J. Anderson, T.C. Chang, Y.C. Li et al., _Erasing the Milky Way: new cleaning technique applied to GBT intensity mapping data_ , _Mon. Not. R. Astron. Soc._ 464 (2017) 4938 [1510.05453].
  * [21] I.P. Carucci, M.O. Irfan and J. Bobin, _Recovery of 21-cm intensity maps with sparse component separation_ , _Mon. Not. R. Astron. Soc._ 499 (2020) 304 [2006.05996].
  * [22] S.R. Furlanetto, S.P. Oh and F.H. Briggs, _Cosmology at low frequencies: The 21 cm transition and the high-redshift Universe_ , _Phys. Rep._ 433 (2006) 181 [astro-ph/0608032].
  * [23] M.F. Morales, B. Hazelton, I. Sullivan and A. Beardsley, _Four Fundamental Foreground Power Spectrum Shapes for 21 cm Cosmology Observations_ , _Astrophys. J._ 752 (2012) 137 [1202.3830].
  * [24] Planck Collaboration, P.A.R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont et al., _Planck 2015 results. XIII. Cosmological parameters_ , _Astron. Astrophys._ 594 (2016) A13 [1502.01589].
  * [25] C.M. Trott, R.B. Wayth and S.J. Tingay, _The Impact of Point-source Subtraction Residuals on 21 cm Epoch of Reionization Estimation_ , _Astrophys. J._ 757 (2012) 101 [1208.0646].
  * [26] J.S. Dillon, _Power Spectrum Estimation for 21 cm Cosmology_ , in _American Astronomical Society Meeting Abstracts #224_ , vol. 224 of _American Astronomical Society Meeting Abstracts_ , p. 318.04, June, 2014\.
  * [27] A. Liu, A.R. Parsons and C.M. Trott, _Epoch of reionization window. I. Mathematical formalism_ , _Phys. Rev. D_ 90 (2014) 023018 [1404.2596].
  * [28] J.C. Pober, A. Liu, J.S. Dillon, J.E. Aguirre, J.D. Bowman, R.F. Bradley et al., _What Next-generation 21 cm Power Spectrum Measurements can Teach us About the Epoch of Reionization_ , _Astrophys. J._ 782 (2014) 66 [1310.7031].
  * [29] J.C. Pober, _The impact of foregrounds on redshift space distortion measurements with the highly redshifted 21-cm line_ , _Mon. Not. R. Astron. Soc._ 447 (2015) 1705 [1411.2050].
  * [30] N. Kaiser, _On the spatial correlations of Abell clusters._ , _Astrophys. J. Lett._ 284 (1984) L9.
  * [31] N. Kaiser, _Clustering in real space and in redshift space_ , _Mon. Not. R. Astron. Soc._ 227 (1987) 1.
  * [32] D. Jeong and F. Schmidt, _Large-scale structure observables in general relativity_ , _Classical and Quantum Gravity_ 32 (2015) 044001 [1407.7979].
  * [33] V. Desjacques, D. Jeong and F. Schmidt, _Large-scale galaxy bias_ , _Phys. Rep._ 733 (2018) 1 [1611.09787].
  * [34] M. McQuinn, O. Zahn, M. Zaldarriaga, L. Hernquist and S.R. Furlanetto, _Cosmological Parameter Estimation Using 21 cm Radiation from the Epoch of Reionization_ , _Astrophys. J._ 653 (2006) 815 [astro-ph/0512263].
  * [35] Y. Mao, M. Tegmark, M. McQuinn, M. Zaldarriaga and O. Zahn, _How accurately can 21cm tomography constrain cosmology?_ , _Phys. Rev. D_ 78 (2008) 023529 [0802.1710].
  * [36] P. Bull, P.G. Ferreira, P. Patel and M.G. Santos, _Late-time Cosmology with 21 cm Intensity Mapping Experiments_ , _Astrophys. J._ 803 (2015) 21 [1405.1452].
  * [37] D. Karagiannis, A. Slosar and M. Liguori, _Forecasts on primordial non-Gaussianity from 21 cm intensity mapping experiments_ , _JCAP_ 2020 (2020) 052 [1911.03964].
  * [38] U.-L. Pen, R. Sheth, J. Harnois-Deraps, X. Chen and Z. Li, _Cosmic Tides_ , _arXiv e-prints_ (2012) arXiv:1202.5804 [1202.5804].
  * [39] H.-M. Zhu, U.-L. Pen, Y. Yu, X. Er and X. Chen, _Cosmic tidal reconstruction_ , _Phys. Rev. D_ 93 (2016) 103504 [1511.04680].
  * [40] F. Schmidt, E. Pajer and M. Zaldarriaga, _Large-scale structure and gravitational waves. III. Tidal effects_ , _Phys. Rev. D_ 89 (2014) 083507 [1312.5616].
  * [41] N.G. Karaçaylı and N. Padmanabhan, _Anatomy of cosmic tidal reconstruction_ , _Mon. Not. R. Astron. Soc._ 486 (2019) 3864 [1904.01387].
  * [42] H.-M. Zhu, T.-X. Mao and U.-L. Pen, _Cosmic Tidal Reconstruction with Halo Fields_ , _Astrophys. J._ 929 (2022) 5 [2108.01575].
  * [43] S.-H. Zang, H.-M. Zhu, M. Schmittfull and U.-L. Pen, _Cosmic tidal reconstruction in redshift space_ , _arXiv e-prints_ (2022) arXiv:2212.04294 [2212.04294].
  * [44] D. Jeong and M. Kamionkowski, _Clustering Fossils from the Early Universe_ , _Phys. Rev. Lett._ 108 (2012) 251301 [1203.0302].
  * [45] L. Dai, D. Jeong and M. Kamionkowski, _Seeking inflation fossils in the cosmic microwave background_ , _Phys. Rev. D_ 87 (2013) 103006 [1302.1868].
  * [46] E. Dimastrogiovanni, M. Fasiello, D. Jeong and M. Kamionkowski, _Inflationary tensor fossils in large-scale structure_ , _JCAP_ 2014 (2014) 050 [1407.8204].
  * [47] P. Li, S. Dodelson and R.A.C. Croft, _Large scale structure reconstruction with short-wavelength modes_ , _Phys. Rev. D_ 101 (2020) 083510 [2001.02780].
  * [48] P. Li, R.A.C. Croft and S. Dodelson, _New Probes of Large Scale Structure_ , _arXiv e-prints_ (2020) arXiv:2007.00226 [2007.00226].
  * [49] F. Bernardeau, S. Colombi, E. Gaztañaga and R. Scoccimarro, _Large-scale structure of the Universe and cosmological perturbation theory_ , _Phys. Rep._ 367 (2002) 1 [astro-ph/0112551].
  * [50] A. Taruya, T. Nishimichi and D. Jeong, _Grid-based calculation for perturbation theory of large-scale structure_ , _Phys. Rev. D_ 98 (2018) 103532 [1807.04215].
  * [51] A. Taruya, T. Nishimichi and D. Jeong, _Covariance of the matter power spectrum including the survey window function effect: N -body simulations versus fifth-order perturbation theory on grids_ , _Phys. Rev. D_ 103 (2021) 023501 [2007.05504].
  * [52] A. Taruya, T. Nishimichi and D. Jeong, _Grid-based calculations of redshift-space matter fluctuations from perturbation theory: UV sensitivity and convergence at the field level_ , _Phys. Rev. D_ 105 (2022) 103507 [2109.06734].
  * [53] D. Jeong, _Cosmology with high ( $z>1$) redshift galaxy surveys_, Ph.D. thesis, University of Texas, Austin, August, 2010.
  * [54] C.-T. Chiang, C. Wagner, F. Schmidt and E. Komatsu, _Position-dependent power spectrum of the large-scale structure: a novel method to measure the squeezed-limit bispectrum_ , _JCAP_ 2014 (2014) 048 [1403.3411].
  * [55] C.-T. Chiang, _Position-dependent power spectrum: a new observable in the large-scale structure_ , _arXiv e-prints_ (2015) arXiv:1508.03256 [1508.03256].
  * [56] S. Adhikari, D. Jeong and S. Shandera, _Constraining primordial and gravitational mode coupling with the position-dependent bispectrum of the large-scale structure_ , _Phys. Rev. D_ 94 (2016) 083528 [1608.05139].
  * [57] D. Jeong and F. Schmidt, _Parity-odd galaxy bispectrum_ , _Phys. Rev. D_ 102 (2020) 023530.
  * [58] W. Hu and T. Okamoto, _Mass Reconstruction with Cosmic Microwave Background Polarization_ , _Astrophys. J._ 574 (2002) 566 [astro-ph/0111606].
  * [59] T. Okamoto and W. Hu, _Cosmic microwave background lensing reconstruction on the full sky_ , _Phys. Rev. D_ 67 (2003) 083002 [astro-ph/0301031].
  * [60] D. Jeong and E. Komatsu, _Perturbation Theory Reloaded: Analytical Calculation of Nonlinearity in Baryonic Oscillations in the Real-Space Matter Power Spectrum_ , _Astrophys. J._ 651 (2006) 619 [astro-ph/0604075].
  * [61] J. Tomlinson and D. Jeong, _Spherical bispectrum: a novel visualization scheme for facilitating comparisons_ , _JCAP_ 2023 (2023) 040 [2204.00668].
  * [62] O. Darwish, S. Foreman, M.M. Abidi, T. Baldauf, B.D. Sherwin and P.D. Meerburg, _Density reconstruction from biased tracers and its application to primordial non-Gaussianity_ , _Phys. Rev. D_ 104 (2021) 123520 [2007.08472].
  * [63] E. Sefusatti and R. Scoccimarro, _Galaxy bias and halo-occupation numbers from large-scale clustering_ , _Phys. Rev. D_ 71 (2005) 063001 [astro-ph/0412626].
  * [64] D. Gualdi and L. Verde, _Integrated trispectrum detection from BOSS DR12 NGC CMASS_ , _JCAP_ 2022 (2022) 050 [2201.06932].
  * [65] M. Crocce, S. Pueblas and R. Scoccimarro, _Transients from initial conditions in cosmological simulations_ , _Mon. Not. R. Astron. Soc._ 373 (2006) 369 [astro-ph/0606505].
  * [66] A. Lewis, A. Challinor and A. Lasenby, _Efficient Computation of Cosmic Microwave Background Anisotropies in Closed Friedmann-Robertson-Walker Models_ , _Astrophys. J._ 538 (2000) 473 [astro-ph/9911177].
  * [67] E.T. Vishniac, _Why weakly non-linear effects are small in a zero-pressure cosmology_ , _Mon. Not. R. Astron. Soc._ 203 (1983) 345.
  * [68] J.N. Fry, _The Galaxy correlation hierarchy in perturbation theory_ , _Astrophys. J._ 279 (1984) 499.
  * [69] M.H. Goroff, B. Grinstein, S.J. Rey and M.B. Wise, _Coupling of modes of cosmological mass density fluctuations_ , _Astrophys. J._ 311 (1986) 6.
  * [70] Y. Suto and M. Sasaki, _Quasinonlinear theory of cosmological self-gravitating systems_ , _Phys. Rev. Lett._ 66 (1991) 264.
  * [71] N. Makino, M. Sasaki and Y. Suto, _Analytic approach to the perturbative expansion of nonlinear gravitational fluctuations in cosmological density and velocity fields_ , _Phys. Rev. D_ 46 (1992) 585.
  * [72] B. Jain and E. Bertschinger, _Second-Order Power Spectrum and Nonlinear Evolution at High Redshift_ , _Astrophys. J._ 431 (1994) 495 [astro-ph/9311070].
  * [73] R. Scoccimarro and J.A. Frieman, _Loop Corrections in Nonlinear Cosmological Perturbation Theory. II. Two-Point Statistics and Self-Similarity_ , _Astrophys. J._ 473 (1996) 620 [astro-ph/9602070].
  * [74] N. McCullagh, D. Jeong and A.S. Szalay, _Toward accurate modelling of the non-linear matter bispectrum: standard perturbation theory and transients from initial conditions_ , _Mon. Not. R. Astron. Soc._ 455 (2016) 2945 [1507.07824].
  * [75] Z. Wang, D. Jeong, A. Taruya, T. Nishimichi and K. Osato, _Perturbation theory remixed: Improved nonlinearity modeling beyond standard perturbation theory_ , _Phys. Rev. D_ 107 (2023) 103534 [2209.00033].
  * [76] C. Wagner, F. Schmidt, C.T. Chiang and E. Komatsu, _Separate universe simulations._ , _Mon. Not. R. Astron. Soc._ 448 (2015) L11 [1409.6294].
  * [77] A. Barreira and F. Schmidt, _Responses in large-scale structure_ , _JCAP_ 2017 (2017) 053 [1703.09212].
  * [78] U. Seljak, _Extracting Primordial Non-Gaussianity without Cosmic Variance_ , _Phys. Rev. Lett._ 102 (2009) 021302 [0807.1770].
  * [79] P. McDonald and U. Seljak, _How to evade the sample variance limit on measurements of redshift-space distortions_ , _JCAP_ 2009 (2009) 007 [0810.0323].