Anatomy of Cosmic Tidal Reconstruction

TL;DR

This paper reviews cosmic tidal reconstruction techniques for 21-cm intensity surveys, demonstrating their robustness and effectiveness in recovering large-scale signals despite foreground contamination and mode loss.

Contribution

It provides a theoretical framework and estimator construction for tidal reconstruction, analyzing its robustness and quantifying the impact of foreground wedge removal.

Findings

Reconstruction maintains high cross-correlation (r > 0.7) on large scales with mode cutoff.

Removing wedge modes reduces correlation but still retains significant signal (0.3 < r < 0.5).

Theoretical predictions align well with numerical simulations.

Abstract

21-cm intensity surveys aim to map neutral hydrogen atoms in the universe through hyper-fine emission. Unfortunately, long-wavelength (low-wavenumber) radial modes are highly contaminated by smooth astrophysical foregrounds that are six orders of magnitude brighter than the cosmological signal. This contamination also leaks into higher radial and angular wavenumber modes and forms a foreground wedge. Cosmic tidal reconstruction aims to extract the large-scale signal from anisotropic features in the local small-scale power spectrum through non-linear tidal interactions; losing small-scale modes to foreground wedge will impair its performance. In this paper, we review tidal interaction theory and estimator construction, and derive the theoretical expressions for the reconstructed spectra. We show the reconstruction is robust against peculiar velocities. Removing low line-of-sight $k$ modes, we demonstrate cross-correlation coefficient $r$ is greater than 0.7 on large scales ($k <0.1$ $h/$Mpc) even with a cutoff value $k^c_{\|}=0.1$ $h/$Mpc. Discarding wedge modes yields $0.3< r < 0.5$ and completely removes the dependency on $k^c_{\|}$. Our theoretical predictions agree with these numerical simulations.

Figures

Equations

$$\varphi (\text{x})$$

$$=\int {\text{displaylimits}}_{k<k_L}\frac{d^3\text{k}}{(2\pi {)}^3}\varphi (\text{k})e^{i\text{k}\cdot \text{x}}+\int {\text{displaylimits}}_{k>k_L}\frac{d^3\text{k}}{(2\pi {)}^3}\varphi (\text{k})e^{i\text{k}\cdot \text{x}}$$

$$={\Phi }_L(\text{x},k_L)+{\Phi }_S(\text{x},k_L)$$

$${\Phi }_L(\text{x})$$

$$=\int {\text{displaylimits}}_{k<k_L}\frac{d^3\text{k}}{(2\pi {)}^3}\varphi (\text{k})\left(1+i\text{k}\cdot \text{x}-\frac{1}{2}{k}_i{k}_j{x}^i{x}^j\right)$$

$$={\Phi }_L(\text{0})+\text{x}\cdot \nabla {\Phi }_L(\text{0})+\frac{1}{2}{x}^i{x}^j{\Phi }_{L,ij}(\text{0}).$$

$$\varphi (\text{x})={\Phi }_S+{ϵ}_t\frac{1}{2}{t}_{ij}(\tau )x^i{x}^j,$$

$$\left[\frac{d^2}{d{\tau }^2}+\mathcal{H}\frac{d}{d\tau }\right]\text{x}\equiv \mathbb{D}\text{x}=-{\nabla }_x\varphi ,$$

$$=\text{q}+\Psi (\text{q},\tau )$$

$$\Psi$$

$$\delta$$

$$=-\text{Tr}{\Psi }_{i,j}+\frac{1}{2}[(\text{Tr}{\Psi }_{i,j}{)}^2+\text{Tr}({\Psi }_{i,j}^2)],$$

$${\delta }_{\text{LG}}(\text{q})={ϵ}_s{\delta }_{1s}+{ϵ}_s{ϵ}_t\left[{\Psi }_{i,j}^{1s}{\Psi }_{j,i}^{1t}-{\Psi }_{i,i}^{1st}\right].$$

$${\delta }_{\text{EU}}(\text{x})$$

$$={\delta }_{\text{LG}}(\text{x})-\Psi \cdot \nabla {\delta }_{\text{LG}}(\text{x})$$

$${\delta }_{\text{EU}}(\text{x})$$

$${\delta }_t(\text{x})$$

$$\mathbb{D}{\Psi }_{i,i}$$

$${\Psi }_i^{1s}$$

$${\Psi }_i^{1t}$$

$${\Psi }_{i,i}^{1st}$$

$$F(\tau )$$

$$D_{1st}(\tau )$$

$${\delta }_t=t_{ij}^{(0)}\left[\alpha (\tau )\frac{{\partial }_i{\partial }_j}{{\nabla }^2}+F(\tau )x_j{\partial }_i\right]{\delta }_{1s}(\text{x},\tau ).$$

$${\tilde{P}}_{1s}(\text{k},\tau )=P_{1s}(k,\tau )+t_{ij}^{(0)}{\hat{k}}^i{\hat{k}}^{j}f(k,\tau )P_{1s}(k,\tau ),$$

$$\Delta P_{1s}$$

$$\mathcal{L}$$

$$\frac{\partial \mathcal{L}}{\partial {\gamma }_1}$$

$$\frac{\partial {\tilde{P}}_{1s}}{\partial {\gamma }_1}$$

$${\hat{\gamma }}_1=\int \frac{d^{3}k}{(2\pi {)}^3}\frac{|\tilde{\delta }(\text{k}){|}^2}{L^3}\frac{P_{1s}(k)}{{\tilde{P}}_{tot}^2(\text{k})}f(k,\tau )({\hat{k}}_1^2-{\hat{k}}_2^2).$$

$${\tilde{\delta }}^{w_i}(\text{k})=\tilde{\delta }(\text{k}){\left[\frac{P_{1s}(k)f(k)}{{\tilde{P}}_{tot}^2(\text{k})}\right]}^{1/2}i{\hat{k}}_i,$$

$${\hat{\gamma }}_1=\int \frac{d^{3}x}{L^3}\left({\tilde{\delta }}^{w_x}(\text{x}){\tilde{\delta }}^{w_x}(\text{x})-{\tilde{\delta }}^{w_y}(\text{x}){\tilde{\delta }}^{w_y}(\text{x})\right)$$

$${\hat{\gamma }}_1(\text{x})$$

$${\hat{\gamma }}_2(\text{x})$$

$${\kappa }_{3D}(\text{k})=\frac{2k^2}{3(k_x^2+k_y^2{)}^2}\left[(k_x^2-k_y^2){\hat{\gamma }}_1(\text{k})+2k_x{k}_y{\hat{\gamma }}_2(\text{k})\right]$$

$$W(k_{\perp },k_{\Vert })=\frac{⟨{\kappa }_{3D}\delta ⟩}{⟨{\kappa }_{3D}{\kappa }_{3D}⟩}.$$

$$\kappa (\text{k})={\kappa }_{3D}(\text{k})W(k_{\perp },k_{\Vert })$$

Full text

††pubyear: 2018††pagerange: Anatomy of Cosmic Tidal Reconstruction–C

1 Introduction

21-cm intensity mapping has been a promising probe for cosmological information on large-scales (Bull et al., 2015). Late-time 21-cm measurements such as CHIME (Bandura et al., 2014) and HIRAX (Newburgh et al., 2016) aim to map large volumes of the universe through red-shifted neutral hydrogen emission line in $0.8<z<2.5$ redshift range. This method faces challenges from foregrounds that are six orders of magnitude brighter than the cosmological signal. Long-wavelength modes in the line-of-sight are contaminated by these spectrally smooth foregrounds. Furthermore, instrumental imperfections cause this noise to leak into high $k_{\bot}-k_{\|}$ modes and to form a foreground wedge in Fourier space (Morales et al., 2012). Since the long-wavelength radial modes play a significant role in cross correlations of 21-cm surveys with CMB measurements and photo-z galaxy surveys (Furlanetto & Lidz, 2007; Adshead & Furlanetto, 2008; Masui et al., 2013), recovering the lost modes is important to improve the cross correlation signal.

A procedure called cosmic tidal reconstruction has been developed to extract these lost modes from small-scale signal (Pen et al., 2012; Zhu et al., 2016, 2018). Although large-scale modes are not explicitly present, they leave an anisotropic imprint on small-scale modes through non-linear tidal interactions. To better visualize the nature of this interaction, we remind the reader the Moon’s tidal field on Earth. The gravitational potential of the Moon around the Earth’s centre can be Taylor expanded up to second order. The constant term can be ignored; the first derivative is the net gravitational force between two celestial objects. The second derivative of the potential causes deformations on Earth’s oceans and is the tidal field. Similarly, the cosmological gravitational potential can be separated into short- and long-wavelength parts. The latter can be expanded up to second order around an observer. This tidal field $t_{ij}\sim\Phi_{L,ij}$ then causes deformations in the local universe and introduces anisotropies to the local small-scale correlation function and power spectrum.

After finding the theoretical expression for the local small-scale power spectrum, we now need to estimate the underlying large-scale over-density field $\delta_{L}$ given an observation $\delta$. Assigning a Gaussian probability to measuring a data set $\delta$ for a given $\delta_{L}$, we can claim that our best estimate maximizes this probability. In technical terms, optimal quadratic estimators for $\delta_{L}$ are constructed assuming the observed density field $\delta$ is Gaussian. This assumption mandates mapping $\delta$ to a Gaussian field. Choosing robust $t_{ij}$ components is the last piece of the puzzle. Galaxies have peculiar velocities on top of the Hubble expansion. The line-of-sight positions are mismeasured because of these peculiar velocities. To minimize the so called redshift space distortions, two quadrupolar distortions of the tidal field in the plane perpendicular to the line-of-sight ($\gamma_{1}=(\Phi_{L,xx}-\Phi_{L,yy})$ and $\gamma_{2}=2\Phi_{L,xy}$) are chosen. Even though both $\gamma_{1}$ and $\gamma_{2}$ are independent estimators of $\delta_{L}$ in theory, they can be combined into the three dimensional convergence field $\kappa_{3D}$ which is a better estimate for $\delta_{L}$ given that each $\gamma$ is zero in distinct $\hat{k}$ direction (Kaiser, 1992). A Wiener filter finally corrects for the bias and noise in $\kappa_{3D}$ itself.

A reconstruction algorithm first has to pass a test on ideal input by producing a highly correlated field. Then, it has to be robust against the obstacles arising from imperfections. Zhu et al. (2016) showed cosmic tidal reconstruction produced cross-correlations greater than 0.9 on scales $k\la 0.1$ $h/$Mpc at $z=0$ for full dark matter field in real space. In this paper, we run the cosmic tidal reconstruction in real and redshift spaces to test the efficiency of the algorithm. We introduce the observational challenges to assess if the reconstruction is robust. In short, our work includes redshift space distortions, testing a range for low $k_{\|}$ data loss and including the foreground wedge.

The reader may be wondering how all these steps produce a correlated signal. The theoretical expectations are concealed behind tidal interactions, estimators and multiple Fourier transforms in implementation. Because of these complications, previous works relied only on simulations. In this work, we derive the analytic expressions for cross and power spectrum. We show that the cross correlations emerge from a modified bispectrum. This modification is due to Gaussian mapping of the initial over-density field. Our theoretical model can quantify the steep decline in reconstruction efficiency due to lost modes and the minor degradation due to peculiar velocities. It also promises a Wiener filter estimate when higher order terms are considered.

In Section 2, we review the tidal interaction theory and the construction of estimators. Section 3 describes the parameters for our N-body simulations and outlines the algorithm step by step. We present our results from these simulations in Section 4. We derive analytic cross and auto spectrum expressions in Section 5; and summarize in Section 6.

2 Theory

In this section, we first review the derivation of the local small-scale power spectrum (Zhu et al., 2016; Schmidt et al., 2014) and connect it to the construction of estimators (Lu & Pen, 2008). We refer the reader to the references for a longer discussion. Moreover, Akitsu et al. (2017) and Akitsu & Takada (2018) study the impact of the large-scale tidal field in a similar but different framework.

2.1 Local Power Spectrum

The first ingredient of the cosmic tidal reconstruction is a theoretical prediction for local anisotropic features. What we call tidal interactions are the mechanism introducing these anisotropies. Zhu et al. has reviewed the traceless tidal field, and Schmidt et al. has studied the tidal interaction theory in detail using conformal Fermi Normal Coordinate frame ($\overline{\text{FNC}}$). We find it sufficient to work in Newtonian picture for our discussion. Suppressing the time dependence, we start by decomposing the gravitational potential in Fourier space to separate it into short- and long-wavelength parts.

[TABLE]

Now we consider a spherical volume of radius R around the origin such that $k_{L}R\la 1$. Suppressing $k_{L}$ we can rewrite the long-wavelength potential as

[TABLE]

We can drop the constant term in the potential. The first derivative is a net force on the local universe; for large enough scales this can be ignored as well. The trace of $\Phi_{L,ij}(\textbf{0})$ describes if the local universe lives at an under- or over-dense region. The effect of this term is an overall change to the growth rate, which is not detectable and we ignore its contribution here. So with a traceless tidal field $t_{ij}(\tau)\equiv\Phi_{L,ij}(\textbf{0})-\delta^{K}_{ij}\Phi_{L,kk}(\textbf{0})/3$, the gravitational potential $\phi$ is given by

[TABLE]

where small-scale gravitational potential $\Phi_{S}$ obeys the Poisson equation $\nabla^{2}_{x}\Phi_{S}=\frac{3}{2}\Omega_{m}(\tau)\mathcal{H}^{2}(\tau)\delta_{S}$. We will suppress subscript $S$ to simplify our notation. We added $\epsilon_{t}$ to keep track of the long-wavelength perturbations in powers of $t_{ij}$. The tidal field can be written as $t_{ij}(\tau)=T(\tau)t^{(0)}_{ij}$, where $t^{(0)}_{ij}$ is the present value, $T(\tau)=D(\tau)/a(\tau)$ and $D(\tau)$ is the linear growth function.

The equation of motion for a particle in an expanding universe is

[TABLE]

where we have defined the operator $\mathbb{D}$. We can solve this equation using Lagrangian perturbation theory up to the nearest order in non-linear coupling between long- and short- wavelength perturbations.

[TABLE]

The hypothetical particle in question still lives in the spherical volume, so $x<R$ by construction. Additionally, $1s$ term is the linear solution and $\epsilon_{t}$ signifies the contribution to the local small-scale fluctuations from the tidal field. We introduced $\epsilon_{s}$ to keep track of the coupling order.

Let us first expand the over-density field using the mass conservation relation.

[TABLE]

where $\partial^{q}_{j}\Psi_{i}\equiv\Psi_{i,j}$. We set $\Psi^{1t}_{i,i}=0$ with hindsight because the tidal field is traceless ($t_{ii}=0$). The corresponding over-density in Lagrangian coordinates is

[TABLE]

We find the over-density field in Eulerian frame $\delta_{\text{EU}}$ by demanding $\delta_{\text{LG}}(\textbf{q})=\delta_{\text{EU}}(\textbf{x})$.

[TABLE]

Then, the tidal field’s contribution to the local density fluctuations $\delta_{t}=\delta_{\text{EU}}-\delta_{1s}$ is

[TABLE]

Now, we go back to the equation of motion. We find the divergence using the chain rule $\partial^{q}_{i}=\partial^{x}_{i}+\Psi_{j,i}\partial^{x}_{j}$.

[TABLE]

We can solve for the displacement field order by order. Here, we only summarize the results.

[TABLE]

and the time evolution functions are given by

[TABLE]

These functions in terms of $z$ integral are in Appendix A.

Putting all these together and defining $\alpha(\tau)\equiv-D_{1st}(\tau)+F(\tau)$, we find $\delta_{t}$ in terms of $\delta_{1s}$ and its spatial derivatives.

[TABLE]

The local small-scale correlation function is defined with respect to a local origin, $\xi(\textbf{r})=\langle\delta(\textbf{0})\delta(\textbf{r})\rangle$ where $\delta=\delta_{1s}+\delta_{t}$. Then, the nearest order power spectrum under tidal distortions is

[TABLE]

where $f(k,\tau)=-2\alpha(\tau)-F(\tau)d\ln P_{l}(k,\tau)/d\ln k$. For clarity, we stress that $1s$ term is the linear solution and $P_{1s}\equiv P_{l}$ is the linear power spectrum. Tilde represents tidal distortions, and $\hat{k}_{i}=k_{i}/k$.

Analogous to weak lensing, we define $\gamma_{1}=(\Phi_{L,xx}-\Phi_{L,yy})$ and $\gamma_{2}=2\Phi_{L,xy}$ and obtain

[TABLE]

Note that we ignore the contribution from the $z$-components of the tidal tensor to minimize the effects of redshift space distortions on our measurements.

2.2 Estimators

To estimate the underlying long-wavelength over-density field, we assign a Gaussian probability distribution for measuring a data set $\tilde{\delta}(\textbf{k})$ that depends on parameters $\gamma_{1}$ and $\gamma_{2}$ (Lu & Pen, 2008): $P(\tilde{\delta}(\textbf{k});\gamma_{1},\gamma_{2})\propto|\mathbb{C}|^{-1/2}\exp(-\tilde{\delta}^{\dagger}\mathbb{C}^{-1}\tilde{\delta}/2)$. However, we pursue the best $\gamma$ estimate of one data set. We define the likelihood function $\mathcal{L}$ as the negative logarithm of the probability distribution over parameters $\gamma$ for a fixed $\tilde{\delta}(\textbf{k})$. We also assume a diagonal covariance matrix such that $\mathbb{C}_{ij}=\langle\tilde{\delta}(k_{i})\tilde{\delta}(k_{j})\rangle-\langle\tilde{\delta}(k_{i})\rangle\langle\tilde{\delta}(k_{j})\rangle=L^{3}\delta_{ij}\tilde{P}_{tot}(k_{i})$. Then, the likelihood function in the continuum limit is

[TABLE]

where an ensemble average gives $\tilde{P}_{tot}(\textbf{k})=\tilde{P}_{1s}(\textbf{k})+P_{N}(\textbf{k})$ and equation (25) provides the theoretical prediction for $\tilde{P}_{1s}(\textbf{k})$. The noise spectrum $P_{N}(\textbf{k})$ includes the effects of unmodelled non-linearities as well as instrumental noise.

We can construct the estimators by maximizing $\mathcal{L}$ with respect to $\gamma_{1,2}$. Let us go through the calculation for $\gamma_{1}$.

[TABLE]

Expanding the expression in parentheses to first order in estimates $\hat{\gamma}_{1,2}$ and limiting ourselves to quadratic estimators, we find

[TABLE]

We implicitly absorbed the normalization coefficient into the Wiener filter. Furthermore, we can rewrite this expression by first defining

[TABLE]

and then inverse Fourier transforming $\tilde{\delta}^{w_{i}}$. Integral over $k$ yields a Dirac delta function cancelling the other position variable.

[TABLE]

We have defined the tidal field as the second derivatives of the gravitational potential at an observer, so the tidal field and consequently $\hat{\gamma}_{i}$ are not functions of x. Nevertheless, the estimate $\hat{\gamma}_{1}$ is an average over local values as equation (32) indicates. Therefore, the estimate for $\gamma_{1}$ at a location x is the expression in parentheses. These estimators are effectively weighted derivatives. Following the same steps, we find an expression for $\hat{\gamma}_{2}$.

[TABLE]

Underlying density fluctuations and $\gamma_{i}$ are theoretically related through $k^{2}\gamma_{1}\propto(k_{x}^{2}-k_{y}^{2})\delta_{L}/2$ and $k^{2}\gamma_{2}\propto k_{x}k_{y}\delta_{L}$, so either $\gamma_{i}$ can be used to estimate $\delta_{L}$. Given $\gamma_{1}$ vanishes when $k_{x}=k_{y}$ and $\gamma_{2}$ vanishes when either $k_{x}$ or $k_{y}$ is zero, a better estimate for $\delta_{L}$ is the three dimensional convergence field (Kaiser, 1992).

[TABLE]

However, we still need to filter out noise and correct for multiplicative biases in $\kappa_{3D}$. The Wiener filter is constructed by minimizing the error $e=\langle(W\kappa_{3D}-\delta)^{2}\rangle$ with respect to $W$, which yields

[TABLE]

We compute these expectation values by direct simulations, although we also present an analytic approach in Section 5.

3 Implementation

3.1 Simulations

We run 10 simulations with $1024^{3}$ dark matter particles using GADGET-2111http://wwwmpa.mpa-garching.mpg.de/gadget/ (Springel et al., 2001; Springel, 2005) in a box of side length $L=1.5$ Gpc$/h$ with periodic boundary conditions. The boxes have the following cosmological parameters: $\Omega_{m}=0.276$, $\Omega_{b}=0.045$, $\Omega_{\Lambda}=0.724$, $h=0.7$, $n_{s}=0.961$ and $\sigma_{8}=0.811$. The simulations start from $z=49$ using 2LPT initial conditions (Scoccimarro, 1998; Jenkins, 2010) constructed with the linear power spectrum from CAMB222https://lambda.gsfc.nasa.gov/toolbox/tb_camb_form.cfm.

In the following sections, we compare the reconstruction results at $z=0$ and $z=1$ in real and redshift space without any foreground subtraction. For every other case, we run the reconstruction only at $z=1$.

3.2 Method

The tidal reconstruction algorithm presented in Zhu et al. (2016) consists of two phases. The first phase computes $\kappa_{3D}$ from noisy over-density field, while the second phase applies a Wiener filter to $\kappa_{3D}$ to obtain the clean reconstructed field $\kappa$.

In the first phase:

We interpolate the density field to a grid and smooth it with a Gaussian window function. 2. 2.

We apply Gaussian mapping333Even though the log transform is a fast and simple Gaussianization procedure, the foreground subtraction models will cause $\delta<-1$ in our simulations, which hinders evaluating $\ln(1+\delta)$. Also considering the possible numerical errors at void regions, we have chosen to map the field into a Gaussian distribution by preserving the ranking. Given that any interferometic map will have regions $\delta<-1$ with or without foreground subtraction, this mapping will be required in real data as well. while keeping the standard deviation same ($\sigma(\delta_{G})=\sigma(\delta_{R})$) to obtain the Gaussianized over-density field $\delta_{G}(\textbf{x})$ (Weinberg, 1992). 3. 3.

We construct $\delta_{G}^{w_{i}}(\textbf{x})$ using equation (31). 4. 4.

We estimate $\gamma$’s from equations (33) and (34). 5. 5.

We use equation (35) to determine the three dimensional convergence field $\kappa_{3D}$.

Note that 3 and 5 are done in Fourier space, while 4 is a configuration space operation. As is normal, all smoothing steps are done in Fourier space. We use FFTW444http://www.fftw.org for Fourier transforms on $1536^{3}$ grids and $R=1.5$ Mpc/h as our smoothing radius, enough to suppress the interpolation kernel.

To estimate $\delta_{G}^{w_{i}}$, we adopt $P_{1s}(k)$ as the linear power spectrum from CAMB and non-linear power spectrum for $\tilde{P}_{tot}(\textbf{k})$. In redshift space we use a simple Kaiser form: $\tilde{P}^{(s)}_{tot}(\textbf{k})=(1+f\mu^{2})^{2}\tilde{P}_{tot}(\textbf{k})$. The reconstruction works on $xy$ slices and does not have information in purely $z$ direction. As a result $\kappa_{3D}(k_{\bot}=0,k_{z})$ modes go to infinity and contain only noise. We set these modes to zero. Since peculiar velocities cause distortions in $x_{\|}$ direction, we expect $\gamma_{1,2}$ not to be particularly affected by redshift space distortions.

$\kappa_{3D}$ is a biased estimate of $\kappa$; it is noisy and missing a normalization constant. We construct Wiener filter from ten $\langle\kappa_{3D}\delta\rangle$ and $\langle\kappa_{3D}\kappa_{3D}\rangle$ pairs. To get better statistics, we bin the spectra in $(k_{\bot},k_{\|})$. We linearly interpolate the Wiener filter in 2D using GSL555https://www.gnu.org/software/gsl/. In the second phase, the reconstructed field $\kappa$ is the filtered three dimensional convergence field.

[TABLE]

We stress the ensemble (and angle) average in the Wiener filter expression. This filtering does not perfectly recover the true over-density field.

4 Results

Throughout this paper we plot the mean cross-correlation coefficients from N-body simulations between true over-density field $\delta$ and the reconstructed field $\kappa$, where $r_{\kappa\delta}=P_{\kappa\delta}/\sqrt{P_{\delta}P_{\kappa}}$. Error bars are computed from the standard deviation of ten simulation results.

4.1 Tests on the Reconstruction

We run the reconstruction on full dark matter field in real and redshift spaces at $z=0$ and $z=1$. On the top panel of Fig. 1, solid lines represents the results for $z=1$, whereas dashed lines represents the results for $z=0$. The results are similar for all cases, so we explicitly show data points and error bars only for $z=1$. The redshift space distortions (red circles) degrade the correlation coefficient by less than 5%. Given the reconstruction relies on non-linear coupling, its efficiency depends on redshift in a complicated way. For example, purely linear regime ($z\rightarrow\infty$) has no coupling to exploit, whereas late times has higher order coupling besides nearest order. Nevertheless, we do not see a strong improvement at higher redshift. Overall, the reconstruction keeps the correlations high ($\gtrsim 0.7$) up to $k\approx 0.1$ $h/$Mpc. However, our results are approximately $10\%$ worse than the results of Zhu et al. (2016). We attribute this difference to a combination of our lower particle density and higher smoothing radius.

The reconstruction is based on non-linear tidal interactions, so it needs small-scale modes to work efficiently. These modes are limited by the finite beam size, foreground wedge effect and the smoothing radius. Increasing the smoothing radius, we confirm the findings of Zhu et al.: The smoothing radius $R$ should not be larger than 5 Mpc$/h$. To simulate the finite beam, we remove high $k_{\bot}$ modes since they cannot be well resolved. We implement these as in the foreground case described below. Note that this cutoff makes small-scale signal weaker, but unlike Gaussian smoothing it is a real loss of information and not isotropic. The bottom panel of Fig. 1 demonstrates the reconstruction mostly uses modes with $k_{\bot}\la 1.0$ $h/$Mpc and does not work well without modes with $k_{\bot}\ga 0.3$ $h/$Mpc. Results in the redshift space are similar and are not depicted. In reality, how well an angular scale can be resolved depends on the baseline distribution for the experiment. A HIRAX-like experiment has access to a broad range of angular scales; $k_{\bot}<1$ $h/$Mpc are present, although with diminishing quantity (see White & Padmanabhan, 2017, figure 1). Moreover, the thermal noise power is a function of $k_{\bot}$. Simple estimates yield values between 150-600 Mpc${}^{3}/h^{3}$ at $z=1$ and $k_{\bot}=0.2$ $h/$Mpc (see White & Padmanabhan, 2017, Appendix A) and cross over at $k=0.6$ $h/$Mpc (Zhu et al., 2018). Since our figures and results show simply removing modes, the reality will be between the cases we have presented.

4.2 Foregrounds

In this section, we investigate if the reconstruction can recover absent modes. Foregrounds from galactic and extra-galactic synchrotron and free-free emissions are spectrally smooth and contaminate modes with small $k_{\|}$. Without any angular resolution cutoff, we now turn to the astrophysical foregrounds problem at hand. These foregrounds have been implemented using a high-pass filter in Zhu et al. (2018). However, estimator construction in Section 2.2 has motivated us to take a different approach. Instead of a high-pass filter, we adopt a sharp cutoff in $k_{\|}$ and assume $P_{N}(\textbf{k})$ to be infinity for every mode set to zero. We first remove $k_{z}<k_{\|}^{c}$ modes from the original over-density field. Since these modes are purely noise, we set $\delta_{g}^{w_{i}}(k_{z}<k_{\|}^{c})$ in equation (31) to zero explicitly. The reconstruction recovers discarded modes in configuration space. We stress that these foregrounds cost high $k_{\bot}$ modes as well and reduce the reconstruction’s efficiency.

Instrumental limitations introduce another challenge by allowing low $k_{\|}$ foregrounds to leak into higher $k_{\bot}-k_{\|}$ modes (Morales et al., 2012). This effect forms a foreground wedge and contaminates modes with $|k_{\|}|\leq m(z)k_{\bot}$, where $m(z)=D_{M}(z)H(z)/(1+z)$ and $D_{M}(z)$ is the co-moving angular diameter distance (Seo & Hirata, 2016; White & Padmanabhan, 2017). We set these modes to zero to realize the foreground wedge. This effect costs even more small-scale modes and significantly limits the reconstruction’s efficiency.

We remind the reader of our two-phase procedure. The first phase uses the foreground subtracted over-density field. The second phase reverts to using the true over-density field $\delta$ to construct the Wiener filter.

We test three values for $k^{c}_{\|}$ ($0.02$ $h/$Mpc, $0.05$ $h/$Mpc and $0.1$ $h/$Mpc) in real and redshift space at $z=1$. The correlation coefficient falls sharply after the foreground removal since we remove high $k_{\bot}$ modes in the process. However, we find the results are not truly sensitive to $k^{c}_{\|}$. The reconstruction can recover lost modes up to $k^{c}_{\|}=0.1$ $h/$Mpc. Fig. 2 displays the results in redshift space. Real space results follow the same trend with a better performance similar to Fig. 1.

The reconstruction deteriorates more dramatically in the presence of a foreground wedge as seen in Fig. 2. We show the results for $k^{c}_{\|}=0.1$ $h/$Mpc; changing $k^{c}_{\|}$ has no effect anymore since the additional information in low $k_{\bot}$ triangle is insignificant.

In Fig. 3, we show the cross correlation coefficient in the anisotropic $(k_{\bot},k_{\|})$ plane. The figure on the left is a cartoon for the foreground removal procedure, which does not depict redshift space distortions; and the removed wedge modes are represented in bright red. The middle panel shows the results for the low $k_{\|}$ removal. Modes with $k_{\|}<0.1$ $h/$Mpc are recovered up to $k_{\bot}\approx 0.2$ $h/$Mpc, but the anisotropic nature of reconstruction causes a lower $r(k)$ at $k\approx 0.2$ $h/$Mpc in Fig. 2. Without the wedge modes the reconstruction cannot recover $k_{\bot}\ga 0.1$ $h/$Mpc modes, and its performance is limited around the wedge line with $k_{\bot}\la 0.05$ $h/$Mpc. The performance of the reconstruction on $k_{\|}=0$ plane is shown in Fig. 4. The correlations are significantly higher on this plane as Fig. 3 signifies. These modes weigh more in weak lensing cross correlations.

The level to which the foreground wedge can be corrected will depend on the experiment and is still unknown. So, we prefer isolating the wedge effect from low angular resolution, whereas in reality both effects will be present. The details of foregrounds and available angular modes will depend on the instrument and should be calculated into $P_{N}$ for an exact treatment. Although it would deteriorate the reconstruction further, combining these two does not necessarily negate the recovery, since there is an overlap between two cases.

5 Theoretical Expectations

The cosmic tidal reconstruction combines distinct theoretical pieces such as tidal interactions, maximum likelihood estimators and three dimensional convergence field, then adds another complication by requiring multiple Fourier transforms (see Section 3.2 for details). Due to these difficulties, previous works relied only on simulations to examine the reconstruction. Incorporating all the theoretical segments into one framework, we derive analytic expressions for cross and power spectrum. These expressions yield the correct structure in power spectrum and quantify the expectations in reconstruction efficiency due to lost modes and peculiar velocities. This framework also promises an estimate for Wiener filter when higher order terms are considered.

We start by defining the Gaussianized over-density field $\delta_{G}(\textbf{x})=f_{G}[\delta_{R}(\textbf{x});\sigma]$, where $\delta_{R}$ is smoothed with a Gaussian window function with radius $R$ and $f_{G}$ is the Gaussian mapping. Then, we rewrite equation (31) as $\delta_{G}^{w_{i}}(\textbf{k})=\delta_{G}(\textbf{k})w(\textbf{k})(ik_{i}/k)$. We transform equations (33) and (34) for $\gamma_{i}(\textbf{x})$ into $k$ integrals by substituting Fourier transforms of $\delta_{G}^{w_{i}}(\textbf{x})$.

[TABLE]

where

[TABLE]

We warn the reader that q in this context does not represent Lagrangian coordinates.

Note that $w(\textbf{k})=w(-\textbf{k})$ and therefore $T_{\gamma_{i}}$ is symmetric and has the following properties:

[TABLE]

Now, we find the Fourier transforms of $\gamma$’s.

[TABLE]

The integral over $x$ yields a Dirac delta function which integrates out $p$. Consequently, the local estimates for each $\gamma$ becomes

[TABLE]

in Fourier space. Equation (45) easily extends to an expression for the three dimensional convergence field $\kappa_{3D}(\textbf{k})$. We prefer defining $\textbf{p}\equiv\textbf{k}-\textbf{q}$ to further simplify the notation. Using equation (35), we arrive at

[TABLE]

where $\tilde{T}_{\kappa_{3D}}(\textbf{q},\textbf{p})$ is given by equation (35) with $\tilde{T}_{\gamma_{i}}$ instead of $\gamma_{i}$. The exact form of $\tilde{T}_{\kappa_{3D}}$ can be found in Appendix B.

Our conventions for the cross spectrum and the power spectrum are the same: $\langle\kappa_{3D}(\textbf{k})\delta(\textbf{k}^{\prime})\rangle=(2\pi)^{3}\delta_{D}(\textbf{k}+\textbf{k}^{\prime})P_{\kappa_{3D}\delta}(\textbf{k})$. Having a direct link between $\delta$ and $\kappa_{3D}$, we are now able to find these spectra. The cross spectrum is

[TABLE]

where we have defined a modified bispectrum $B_{T}$ as

[TABLE]

Although assuming $\delta_{G}$ as the linear field $\delta_{l}$ is tempting, this assumption results in a very small bispectrum in low $k$, since $\delta$ also approximates the linear field at these scales. Given $\delta_{l}$ is Gaussian and has zero bispectrum, $B_{T}\propto\langle\delta_{l}^{3}\rangle=0$. From this short discussion, we conclude that second order corrections in $\delta_{G}$ have to be taken into account. This makes sense considering tidal reconstruction makes use of mode coupling terms which do not exist in linear theory. To find these necessary higher order contributions, we need to take a closer look at the Gaussianization procedure. Since the logarithmic transform is easier to handle (Neyrinck et al., 2011; McCullagh et al., 2016), whereas the exact Gaussian mapping is analytically complicated, we assume

[TABLE]

The constant average term is only in $\textbf{k}=0$ mode and can be ignored. We ignore smoothing for simplicity, and refer the reader to Appendix C for expressions with smoothing. Taking the Fourier transform and using second-order perturbation results (Jain & Bertschinger, 1994), we arrive at

[TABLE]

where $F_{2,G}=F_{2}-\frac{1}{2}$ and

[TABLE]

Having these expressions, we can now calculate the modified bispectrum $B_{T}$. Its only difference from normal bispectrum $B$ calculation is $-1/2$ term in $F_{2,G}$, so this step is not too cumbersome. Then,

[TABLE]

where $P_{l}$ is the linear power spectrum and

[TABLE]

A similar calculation shows that the four-point function gives the power spectrum of $\kappa_{3D}$. We remind the reader that the four-point function $\langle\delta^{4}\rangle$ can be decomposed into three cyclic $\langle\delta^{2}\rangle^{2}$ terms and a connected term. The trispectrum is the Fourier transform of the connected four-point correlation function. In this estimation, assuming $\delta_{G}$ to be the linear field yields a zero trispectrum, but a non-zero four-point function with one vanishing cyclic term.

[TABLE]

As we have discussed in the cross spectrum case, $\delta_{G}\approx\delta_{l}$ produces uncorrelated $\kappa_{3D}$. Hence, we would expect this equation to underestimate the true power.

With all limitations in mind, we can build a theoretical Wiener filter using equations (47) and (55), and find the reconstructed power spectrum $P_{\kappa}$.

[TABLE]

5.1 Results

Having performed extensive numerical analyses in Section 4, we now test our theoretical framework to understand the emergent features and limited efficiency of cosmic tidal reconstruction. Our focus is on equations (47) and (55). We include smoothing in our calculations (see Appendix C) and take $R=1.5$ Mpc$/h$ as before.

The volume integration is simpler in spherical coordinates.

[TABLE]

where $\theta$ is the azimuthal angle with respect to the line-of-sight and $v=\cos\theta$. We apply two more transformations. First, we write the $\phi$ integral in terms of $t=\cos\phi$, since this integrand depends only on $\cos\phi$ and it is symmetric under $\phi\rightarrow-\phi$. Second, we use logarithmic spacing in $q$ integration. Putting these together, the volume integration becomes

[TABLE]

We use GSL’s Monte-Carlo integration library to compute these integrals. Assuming the same cosmology and CAMB power spectra as before, we integrate from $q_{i}=10^{-4}$ $h/$Mpc to $q_{f}=2$ $h/$Mpc at $z=0$. We form an equally spaced, 50x50 grid for $(k_{\bot},k_{z})$ with values between $0.001$ $h/$Mpc and $0.5$ $h/$Mpc including both ends.

Fig.  5 shows the theoretical $P_{\kappa}(k_{\bot},k_{z})$ on the bottom panel and a representative simulation on the top panel. While their structures and shapes agree, the theoretical $\kappa$ power spectrum is off by one order of magnitude. As we have discussed in the previous section, our theoretical $P_{\kappa_{3D}}$ estimate lacks significant power from reconstructed field because of the lowest order expansion. Since $P_{\kappa}=P^{2}_{\kappa_{3D}\delta}/P_{\kappa_{3D}}$, the lack of power in $\kappa_{3D}$ results in high theoretical $P_{\kappa}$.

After confirming our theoretical framework in anisotropic power spectrum, we assess which modes contribute the most to the reconstructed field $\kappa$. We start by rewriting the cross spectrum.

[TABLE]

Because $I(q,v)$ is not symmetric with respect to $v$ (i.e. $I(q,v)\neq I(q,-v)$), we define $\bar{I}(q,v)=I(q,v)+I(q,-v)$ such that $v\in[0,1]$. To quantify the importance of modes, we integrate the cross spectrum from $(q_{i},0^{\circ})$ to $(q,\theta)$, then divide this integral by the true value.

[TABLE]

We evaluate the ratio $\alpha(q,\theta)$ at a fixed $(k_{\bot}=0.05h/$Mpc$,k_{z}=0.02h/$Mpc$)$ value at $z=0$ on a grid equally sampled with 100 $q$ points in $(10^{-4},1.5]$ $h/$Mpc and 90 equally spaced $\theta$ points in $(0,90]^{\circ}$. $P_{\kappa_{3D}\delta}$ is evaluated with $q_{f}=2.5$ $h/$Mpc.

The value of a point $(q,\theta)$ in Fig. 6 represents the fraction of cross correlation signal that can be achieved using modes up to $q$ and angle $\theta$ with respect to the line-of-sight. The reconstruction prefers $q_{\bot}$ modes and reaches its peak efficiency at $q\approx 1.2$ $h/$Mpc when modes up to $\theta=\pi/2$ are accessible as the smallest fluctuations are suppressed by Gaussian smoothing. We added the wedge line for $z=1$ for reference. Foreground wedge impedes access to higher $\theta$ modes and restricts the efficiency of reconstruction to approximately $0.2$. Redshift space distortions affect the reconstruction in reverse direction, contaminating low $\theta$ modes. These are not significant in the reconstruction, explaining the weak dependence in Fig. 1.

6 Summary

21-cm intensity mapping surveys lose large-scale radial signal to astrophysical smooth foregrounds. These foregrounds then contaminate smaller-scale radial and angular modes. Recovering the lost data (low $k_{\|}$ modes) will improve BAO measurements and cross correlations of 21-cm signal with CMB measurements and photo-z galaxy surveys.

The large-scale over-density field deforms the local universe around an observer through a non-zero second derivative of the long-wavelength gravitational potential. This deformation imprints anisotropic features in the local small-scale power spectrum and allows us to reconstruct the lost data.

The second derivative of $\Phi_{L}$ gives the tidal field $t_{ij}(\tau)\equiv\Phi_{L,ij}(\textbf{0})-\delta^{K}_{ij}\Phi_{L,kk}(\textbf{0})/3$. We choose two components ($\gamma_{1}$ and $\gamma_{2}$) to minimize redshift space distortions on the reconstruction. The estimators for each $\gamma$ can be constructed assuming the observed $\delta$ is a Gaussian field. In return, this assumption requires us to Gaussianize the observed $\delta$. The three dimensional convergence field $\kappa_{3D}$ is a linear combination of $\gamma_{1}$ and $\gamma_{2}$; and it is a better estimate for the underlying large-scale over-density field. We correct $\kappa_{3D}$ for bias and noise with a Wiener filter constructed from ten N-body simulations.

In this work, we reviewed non-linear tidal interactions using Lagrangian perturbation theory and explored the efficiency of reconstructing long-wavelength modes from local small-scale fluctuations. We performed new tests on cosmic tidal reconstruction by adding redshift space distortions, investigating a range of $k_{\|}$ and $k_{\bot}$ data loss and foreground wedge data loss. We also presented a novel theoretical framework to study the reconstruction, which can predict its efficiency and produce theoretical estimates for Wiener filter if improved.

We found the cross-correlation coefficient $r$ between the reconstructed field and true over-density field is above $0.7$ until $k=0.1$ $h/$Mpc in both spaces at two redshifts. Our numerical tests also confirmed the reconstruction is robust against peculiar velocities with minor degradation (less than 5%) in $r$. This validates choosing quadrupolar distortions in the plane perpendicular to the line-of-sight, $\gamma_{1}$ and $\gamma_{2}$.

To assess 21-cm intensity mapping challenges, we tested the reconstruction against astrophysical foregrounds. We modelled spectrally smooth foregrounds as a sharp cutoff below $k^{c}_{\|}$ and found that the cross-correlation coefficient decreases by 20%, but it is not sensitive to the cutoff value when $k^{c}_{\|}<0.1$ $h/$Mpc. However, tidal reconstruction does not recover modes isotropically and specifically does not recover low $k_{\bot}$-high $k_{\|}$ modes. Given this feature, the reconstruction recovers modes $k_{\|}<0.1$ $h/$Mpc up to $k_{\bot}=0.2$ $h/$Mpc. Moreover, smooth foregrounds will leak into higher $k_{\bot}-k_{\|}$ modes due to imperfect instrumentation. We modelled this foreground wedge by removing modes obeying $|k_{\|}|\leq k_{\bot}D_{M}(z)H(z)/(1+z)$. After these modes are subtracted, $k^{c}_{\|}$ becomes irrelevant and the reconstruction deteriorates to $r(k<0.1h/$Mpc$)\la 0.5$. In the anisotropic picture, modes around the wedge slope can be recovered with $r\sim 0.5$. On $k_{\|}=0$ plane, the wedge case performance shows mild improvement, whereas other cases do significantly better.

Finally, we incorporated all the theoretical segments into one framework, which showed that the cross correlations arise from a modified bispectrum and the power spectrum of $\kappa_{3D}$ originates from the four-point function. Our framework revealed a similar structure in $P_{\kappa}$, but missed the reconstructed signal due to lowest order expansion in the four-point function. Using the same framework, we estimated what modes are most important for efficient reconstruction. Our theoretical investigation predicted performance variations: steep decline when wedge modes are lost and minor degradation in redshift space. We expect one-loop corrections will decrease the discrepancy between simulations and theoretical predictions; it could construct a capable Wiener filter as well.

We have investigated the angular resolution and foregrounds individually. The available modes will depend on the details of the instrument. An exact treatment should be calculated into $P_{N}$ for real data applications. Since the modes above the wedge are less valuable, this should still produce correlations $r\sim 0.4$.

Even though combining low angular resolution and the foreground wedge will not erase the correlations completely, we find the overall performance not as high as we would like. These challenges hint at utilizing the distortions in the radial direction. The trade-off for using all $t_{ij}$ components is to extend the study to redshift space. Since small scales matter, the instrument design should take into account of resolving smaller scales and mitigating foreground wedge as much as possible as well. We hope that our work stresses the necessity and importance of such further investigations.

Throughout our work we assumed neutral hydrogen perfectly traces dark matter. The neutral hydrogen bias should be modelled for further investigation. A simple proposal would be to take haloes as tracers.

Acknowledgements

N. G. K. and N. P. are partially supported by the DOE DE-SC0017660.

We thank Dongzi Li, Hongming Zhu and Ue-Li Pen for helpful discussion.

Appendix A Time Dependent Functions

The redshift dependent function $f$ is given by

[TABLE]

where $D(z)$ is the linear growth function and

[TABLE]

We defined the following intermediate functions to simplify expressions,

[TABLE]

Appendix B Details on Reconstruction Theory

Three dimensional convergence field is given by the integration in equation (46). We start by writing $\tilde{T}_{\kappa_{3D}}$.

[TABLE]

Since $T_{\gamma_{i}}$ is symmetric, $\tilde{T}_{\kappa_{3D}}(\textbf{k}^{\prime},\textbf{k}-\textbf{k}^{\prime})$ is symmetric as well.

[TABLE]

We can reshape this expression into the following:

[TABLE]

We want to evaluate $\tilde{T}_{\kappa_{3D}}(\textbf{q},\textbf{k}-\textbf{q})$ in spherical coordinates. We start by fixing $\textbf{k}_{\bot}$ and assigning an angle $\phi$ between $\textbf{k}_{\bot}$ and $\textbf{q}_{\bot}$. We continue using $\textbf{p}\equiv\textbf{k}-\textbf{q}$ and define $t\equiv\cos\phi$. Then, the following relations hold:

[TABLE]

We further define $q_{z}\equiv q\cos\theta=qv$. Then $q_{\bot}=q\sqrt{1-v^{2}}$. We can rewrite the relations above in terms of these new spherical coordinates.

[TABLE]

Finally, $\tilde{T}_{\kappa_{3D}}$ can be expressed in terms of the same coordinates.

[TABLE]

Appendix C Smoothing

Let us consider the smoothing: $\delta_{R}(\textbf{k})=S_{R}(k)\delta(\textbf{k})$. Then,

[TABLE]

where we have modified the kernel to include smoothing.

[TABLE]

There are two modifications to previous modified bispectrum $B_{T}$ expression. First, we multiply $B_{T}$ by $S_{R}(q)S_{R}(p)$. Second, we substitute $F_{2,G}\rightarrow F_{2,G,R}$.

[TABLE]