Iterative initial condition reconstruction

TL;DR

This paper introduces an iterative algorithm to reconstruct initial cosmic density conditions from dark matter distributions, significantly improving the accuracy of large-scale structure analysis and BAO measurements.

Contribution

The paper presents a novel iterative reconstruction method that effectively reverses nonlinear gravitational effects up to one-loop order, enhancing the recovery of initial conditions from dark matter data.

Findings

Reconstruction achieves over 95% correlation with initial density at certain scales.

Reduces differences between reconstructed and initial fields by over 100 times at specific scales.

Improves BAO signal-to-noise ratio by approximately 2.5 to 2.7 times, approaching the linear density limit.

Abstract

Motivated by recent developments in perturbative calculations of the nonlinear evolution of large-scale structure, we present an iterative algorithm to reconstruct the initial conditions in a given volume starting from the dark matter distribution in real space. In our algorithm, objects are first moved back iteratively along estimated potential gradients, with a progressively reduced smoothing scale, until a nearly uniform catalog is obtained. The linear initial density is then estimated as the divergence of the cumulative displacement, with an optional second-order correction. This algorithm should undo nonlinear effects up to one-loop order, including the higher-order infrared resummation piece. We test the method using dark matter simulations in real space. At redshift $z=0$, we find that after eight iterations the reconstructed density is more than $95\%$ correlated with the initial density at $k\le 0.35\; h\mathrm{Mpc}^{-1}$. The reconstruction also reduces the power in the difference between reconstructed and initial fields by more than 2 orders of magnitude at $k\le 0.2\; h\mathrm{Mpc}^{-1}$, and it extends the range of scales where the full broadband shape of the power spectrum matches linear theory by a factor of 2-3. As a specific application, we consider measurements of the baryonic acoustic oscillation (BAO) scale that can be improved by reducing the degradation effects of large-scale flows. In our idealized dark matter simulations, the method improves the BAO signal-to-noise ratio by a factor of 2.7 at $z=0$ and by a factor of 2.5 at $z=0.6$, improving standard BAO reconstruction by $70\%$ at $z=0$ and $30\%$ at $z=0.6$, and matching the optimal BAO signal and signal-to-noise ratio of the linear density in the same volume. For BAO, the iterative nature of the reconstruction is the most important aspect.

Tables

Table 1. Table 1: Overview of some recently proposed BAO reconstruction techniques that estimate first a displacement field from initial to final conditions, and then use that to infer the linear density. The methods in the bottom three rows are new to this paper, and a crucial ingredient separating them from other methods is that the smoothing scale is progressively reduced with the number of iteration steps to access smaller scales. The list is by no means complete; see main text for more methods.

Reconstruction algorithm		Stage I: Estimate displacement		Stage II: Estimate linear density
Standard reconstruction Eisenstein et al. (2007a)		Zeldovich ${\widehat{𝝌}}_{\mathrm{ZA}}=\frac{i𝒌}{k^2}W{\delta }_{\mathrm{NL}}$		${\hat{\delta }}_0={\delta }_d\left[{\widehat{𝝌}}_{\mathrm{ZA}}\right]-{\delta }_s\left[{\widehat{𝝌}}_{\mathrm{ZA}}\right]$
Iterative standard rec Seo et al. (2010)		Iterative Zeldovich with fixed $W$		${\hat{\delta }}_0={\delta }_d\left[{\widehat{𝝌}}_{\mathrm{ZA}}\right]-{\delta }_s\left[{\widehat{𝝌}}_{\mathrm{ZA}}\right]$
Improved standard rec Tassev and Zaldarriaga (2012)		Zeldovich and iterative Newton-Raphson		${\hat{\delta }}_0={\hat{\delta }}_{\chi }=\nabla \cdot \widehat{𝝌}$
Standard rec with pixels Seo and Hirata (2016); Obuljen et al. (2016)		Zeldovich ${\widehat{𝝌}}_{\mathrm{ZA}}=\frac{i𝒌}{k^2}W{\delta }_{\mathrm{NL}}$		Move pixels instead of galaxies
Eulerian growth-shift rec Schmittfull et al. (2015)		Zeldovich ${\widehat{𝝌}}_{\mathrm{ZA}}=\frac{i𝒌}{k^2}W{\delta }_{\mathrm{NL}}$		${\hat{\delta }}_0={\delta }_{\mathrm{NL}}-{\widehat{𝝌}}_{\mathrm{ZA}}\cdot \nabla \delta -{\delta }^2$
Nonlinear isobaric rec Zhu et al. (2016a, b)		Solve diff. eqn. with multigrid algorithm		${\hat{\delta }}_0={\hat{\delta }}_{\chi }=\nabla \cdot \widehat{𝝌}$
New $𝒪\left(1\right)$ rec		Iteratively solve $T\left(k\right)F_Z\left[\nabla \cdot 𝝌\right]=W{\delta }_{\mathrm{NL}}$ for $𝝌$		${\hat{\delta }}_0={\hat{\delta }}_{\chi }=\nabla \cdot \widehat{𝝌}$
New $𝒪\left(2\right)$ rec		Iteratively solve $T\left(k\right)F_Z\left[\nabla \cdot 𝝌\right]=W{\delta }_{\mathrm{NL}}$ for $𝝌$		${\hat{\delta }}_0=t_1\left(k\right){\hat{\delta }}_{\chi }+t_2\left(k\right){\int }_{𝒑}{\kappa }_2{\hat{\delta }}_{\chi }\left(𝒑\right){\hat{\delta }}_{\chi }\left(𝒌-𝒑\right)$
Extended std rec [appdx. B]		Iteratively solve $T\left(k\right)F_Z\left[\nabla \cdot 𝝌\right]=W{\delta }_{\mathrm{NL}}$ for $𝝌$		${\hat{\delta }}_0={\delta }_d\left[\widehat{𝝌}\right]-{\delta }_s\left[\widehat{𝝌}\right]$

Table 2. Table 2: Systematic bias of the BAO scale estimated from the best-fit BAO scale from the power spectrum of ten large-volume simulations at z = 0 𝑧 0 z=0 . The BAO scale from the nonlinear density is biased high by 0.3 % percent 0.3 0.3\% . Reconstruction eliminates that bias Padmanabhan  et al. ( 2009 ); Noh  et al. ( 2009 ); Sherwin and Zaldarriaga ( 2012 ); Tassev and Zaldarriaga ( 2012 ); Seo  et al. ( 2010 ); White ( 2015 ); Mehta  et al. ( 2011 ) . The residual biases after reconstruction are small and likely consistent with zero because the estimates are derived from only ten simulations. The left column shows the sample mean of the best-fit BAO scale relative to the fiducial theoretical value, ⟨ r ^ BAO ⟩ − r BAO fid delimited-⟨⟩ subscript ^ 𝑟 BAO superscript subscript 𝑟 BAO fid \langle\hat{r}_{\mathrm{BAO}}\rangle-r_{\mathrm{BAO}}^{\mathrm{fid}} ; the right column is relative to the initial condition of each simulation, ⟨ r ^ BAO − r ^ BAO lin ⟩ delimited-⟨⟩ subscript ^ 𝑟 BAO superscript subscript ^ 𝑟 BAO lin \langle\hat{r}_{\mathrm{BAO}}-\hat{r}_{\mathrm{BAO}}^{\mathrm{lin}}\rangle , canceling cosmic variance.

		Mean BAO scale
Field		vs lin. theory		vs lin. realization
Initial conds.		$+0.05\mathrm{Mpc}$ [$+0.03\%$]		$+0.00\mathrm{Mpc}$ [$\pm 0.00\%$]
Final conds.		$+0.49\mathrm{Mpc}$ [$+0.33\%$]		$+0.44\mathrm{Mpc}$ [$+0.30\%$]
Standard rec		$+0.02\mathrm{Mpc}$ [$+0.01\%$]		$-0.03\mathrm{Mpc}$ [$-0.02\%$]
New $𝒪\left(1\right)$ rec		$+0.05\mathrm{Mpc}$ [$+0.03\%$]		$\pm 0.00\mathrm{Mpc}$ [$\pm 0.00\%$]
New $𝒪\left(2\right)$ rec		$+0.06\mathrm{Mpc}$ [$+0.04\%$]		$+0.02\mathrm{Mpc}$ [$+0.01\%$]

Table 3. Table 3: Left column: Root-mean-square scatter of the best-fit BAO scale between ten 2.6 ​ h − 3 ​ Gpc 3 2.6 superscript ℎ 3 superscript Gpc 3 2.6\,h^{-3}\mathrm{Gpc}^{3} simulations at z = 0 𝑧 0 z=0 . This is a Monte Carlo estimate for the expected statistical 1 ​ σ 1 𝜎 1\sigma uncertainty when measuring the BAO scale from the power spectrum in a single 2.6 ​ h − 3 ​ Gpc 3 2.6 superscript ℎ 3 superscript Gpc 3 2.6\,h^{-3}\mathrm{Gpc}^{3} volume. Right column: Rms scatter of the BAO scale relative to that in the initial conditions of each simulation, r ^ BAO − r ^ BAO lin subscript ^ 𝑟 BAO superscript subscript ^ 𝑟 BAO lin \hat{r}_{\mathrm{BAO}}-\hat{r}_{\mathrm{BAO}}^{\mathrm{lin}} , which is sourced by nonlinear shift terms as discussed in Section  IV.4 . All numbers are somewhat uncertain because they were estimated from the scatter of only ten simulations.

		Rms scatter of BAO scale
Field		vs lin. theory		vs lin. realization
Initial conds.		$0.35\mathrm{Mpc}$ [$0.24\%$]		$0\mathrm{Mpc}$ [$0\%$]
Final conds.		$0.99\mathrm{Mpc}$ [$0.66\%$]		$1.20\mathrm{Mpc}$ [$0.81\%$]
Standard rec		$0.63\mathrm{Mpc}$ [$0.42\%$]		$0.55\mathrm{Mpc}$ [$0.37\%$]
New $𝒪\left(1\right)$ rec		$0.44\mathrm{Mpc}$ [$0.29\%$]		$0.13\mathrm{Mpc}$ [$0.08\%$]
New $𝒪\left(2\right)$ rec		$0.37\mathrm{Mpc}$ [$0.25\%$]		$0.08\mathrm{Mpc}$ [$0.05\%$]

Table 4. Table 4: Systematic bias of the best-fit BAO scale at z = 0 𝑧 0 z=0 as in Table  2 but for the extended standard reconstruction (’Ext std rec’) method of Appendix  B and for the scaled 𝒪 ​ ( 1 ) 𝒪 1 \mathcal{O}(1) reconstruction given by t 1 ​ ( k ) ​ δ χ ​ ( 𝒌 ) subscript 𝑡 1 𝑘 subscript 𝛿 𝜒 𝒌 t_{1}(k)\delta_{\chi}({\bm{k}}) , where δ χ = ∇ ⋅ 𝝌 ( 8 ) subscript 𝛿 𝜒 ⋅ bold-∇ superscript 𝝌 8 \delta_{\chi}={\bm{\nabla}}\cdot\bm{\chi}^{(8)} , and t 1 subscript 𝑡 1 t_{1} is from Appendix  A . There is no evidence for a systematic BAO bias for any of the methods.

		Mean BAO scale
Field		vs lin. theory		vs lin. realization
Ext std rec		$+0.02\mathrm{Mpc}$ [$+0.01\%$]		$-0.03\mathrm{Mpc}$ [$-0.02\%$]
Scaled $𝒪\left(1\right)$ rec		$+0.06\mathrm{Mpc}$ [$+0.04\%$]		$+0.01\mathrm{Mpc}$ [$+0.01\%$]

Table 5. Table 5: Rms scatter of the best-fit BAO scale at z = 0 𝑧 0 z=0 as in Table  3 but for the extended standard reconstruction method of Appendix  B and for the scaled 𝒪 ​ ( 1 ) 𝒪 1 \mathcal{O}(1) reconstruction given by t 1 ​ ( k ) ​ δ χ ​ ( 𝒌 ) subscript 𝑡 1 𝑘 subscript 𝛿 𝜒 𝒌 t_{1}(k)\delta_{\chi}({\bm{k}}) . The left column shows that in both cases the BAO uncertainty is consistent with that of the BAO scale in the initial conditions ( 0.24 % percent 0.24 0.24\% ) within the uncertainty of the ten simulations. The right column represents the nonlinear noise contribution due to residual shift terms after reconstruction as discussed in Section  IV.4 . For BAO, the methods are thus comparable to the method described in the main text.

		Rms scatter of BAO scale
Field		vs lin. theory		vs lin. realization
Ext std rec		$0.31\mathrm{Mpc}$ [$0.21\%$]		$0.18\mathrm{Mpc}$ [$0.12\%$]
Scaled $𝒪\left(1\right)$ rec		$0.4\mathrm{Mpc}$ [$0.27\%$]		$0.11\mathrm{Mpc}$ [$0.07\%$]

Table 6. Table 6: Performance of the second-order reconstruction for some choices of the initial smoothing scale R 𝑅 R , the reduction factor ϵ R subscript italic-ϵ 𝑅 \epsilon_{R} , and the displacement amplitude ϵ s subscript italic-ϵ 𝑠 \epsilon_{s} (all other parameters are set to their default values described in Section  III.1.3 ). The table shows one minus the squared correlation coefficient between reconstructed and linear density, 1 − r 2 ​ ( δ 0 , δ rec ) 1 superscript 𝑟 2 subscript 𝛿 0 subscript 𝛿 rec 1-r^{2}(\delta_{\mathrm{0}},\delta_{\mathrm{rec}}) , after 2, 8 and 16 iteration steps, at k = 0.1 ​ h ​ Mpc − 1 𝑘 0.1 ℎ superscript Mpc 1 k=0.1\;h\text{Mpc}^{-1} and k = 0.6 ​ h ​ Mpc − 1 𝑘 0.6 ℎ superscript Mpc 1 k=0.6\;h\text{Mpc}^{-1} , in our L = 500 ​ h ​ Mpc − 1 𝐿 500 ℎ superscript Mpc 1 L=500\;h\text{Mpc}^{-1} simulation at z = 0.6 𝑧 0.6 z=0.6 . Lower numbers correspond to better reconstruction. Going from 2 to 8 steps improves performance substantially, while more than 8 steps yields diminishing returns. After 8 or 16 steps steps, all reconstruction algorithms with smoothing reduction factor ϵ R = 0.5 subscript italic-ϵ 𝑅 0.5 \epsilon_{R}=0.5 perform similarly well. (Note that the results of the table were obtained using suboptimal transfer functions, which is why results for ( R , ϵ R , ϵ s ) = ( 10 , 0.5 , 1 ) 𝑅 subscript italic-ϵ 𝑅 subscript italic-ϵ 𝑠 10 0.5 1 (R,\epsilon_{R},\epsilon_{s})=(10,0.5,1) differ slightly from Fig.  13 ; the overall trends should not depend on this though.)

				2 steps				8 steps				16 steps
$R$	${ϵ}_R$	${ϵ}_s$		$k=0.1$		$k=0.6$		$k=0.1$		$k=0.6$		$k=0.1$		$k=0.6$
5	0.5	0.5		$5.9\times {10}^{-3}$		0.97		$4.4\times {10}^{-4}$		0.19		$3.0\times {10}^{-4}$		0.17
10	0.5	0.5		$6.5\times {10}^{-3}$		0.9		$4.1\times {10}^{-4}$		0.19		$2.7\times {10}^{-4}$		0.17
10	0.5	1		$6.1\times {10}^{-3}$		0.94		$4.0\times {10}^{-4}$		0.19		$2.9\times {10}^{-4}$		0.18
20	0.5	1		$1.9\times {10}^{-2}$		0.94		$4.1\times {10}^{-4}$		0.2		$2.9\times {10}^{-4}$		0.18
1	1	0.33		$2.2\times {10}^{-2}$		0.86		$3.0\times {10}^{-3}$		0.33		$3.4\times {10}^{-3}$		0.31
10	1	0.5		$8.4\times {10}^{-3}$		0.92		$4.2\times {10}^{-3}$		0.92		$4.0\times {10}^{-3}$		0.92
10	1	1		$8.6\times {10}^{-3}$		0.93		$4.2\times {10}^{-3}$		0.93		$4.1\times {10}^{-3}$		0.93

Equations

$${\delta }_{\mathrm{NL}}=T(k){\delta }_Z[{\delta }_{\chi }]+{\delta }_{error}\equiv {\delta }_{PT}+{\delta }_{error},$$

$${\delta }_{\chi }(\mathbf{k})$$

$$\times \left(1-\frac{({\mathbf{p}}_1\cdot {\mathbf{p}}_2{)}^2}{{\mathbf{p}}_1\cdot {\mathbf{p}}_1{\mathbf{p}}_2\cdot {\mathbf{p}}_2}\right){\delta }_0({\mathbf{p}}_1){\delta }_0({\mathbf{p}}_2)+\cdots .$$

$${\hat{\delta }}_{PT}(\mathbf{k})=W(k){\delta }_{\mathrm{NL}}(\mathbf{k}),$$

$$T(k){\delta }_Z[{\hat{\delta }}_{\chi }]={\hat{\delta }}_{PT}=W{\delta }_{\mathrm{NL}}.$$

$${\hat{\delta }}_0(\mathbf{k})$$

$$\times \left(1-\frac{({\mathbf{p}}_1\cdot {\mathbf{p}}_2{)}^2}{{\mathbf{p}}_1\cdot {\mathbf{p}}_1{\mathbf{p}}_2\cdot {\mathbf{p}}_2}\right){\hat{\delta }}_{\chi }({\mathbf{p}}_1){\hat{\delta }}_{\chi }({\mathbf{p}}_2),$$

$${\delta }^{\prime }+\nabla \cdot [(1+\delta )\mathbf{v}]=0.$$

$$\frac{\mathrm{D}}{\mathrm{D}\tau }\ln (1+\delta )=-\nabla \cdot \mathbf{v}=-\nabla \cdot \frac{\mathrm{d}\mathbf{\chi }}{\mathrm{d}\tau }=-\frac{\mathrm{d}}{\mathrm{d}\tau }\nabla \cdot \mathbf{\chi }.$$

$$\nabla \cdot \mathbf{\chi }=-\ln (1+{\delta }_R)\approx -{\delta }_R.$$

$$R=\max ({ϵ}_R^{n-1}{R}_{\mathrm{init}},R_{\min }).$$

$${\hat{\delta }}_0^{[2]}(\mathbf{k})=t_1(k){\delta }_{\chi }(\mathbf{k})+t_2(k){\delta }_{\chi }^{[2]}(\mathbf{k}).$$

$${\sigma }_{\mathrm{survey}}^2=\frac{V_{\mathrm{sim}}}{V_{\mathrm{survey}}}\left({\sigma }_{\mathrm{lin}}^2+{\sigma }_{\mathrm{NL}}^2\right).$$

$${\hat{\delta }}_0=t_1(k){\delta }_{\chi }(\mathbf{k})+t_2(k){\delta }_{\chi }^{[2]}(\mathbf{k}),$$

$${\delta }_{\chi }^{[2]}(\mathbf{k})={\int }_{{\mathbf{p}}_1}{\kappa }_2({\mathbf{p}}_1,{\mathbf{p}}_2){\bar{t}}_1(p_1){\delta }_{\chi }({\mathbf{p}}_1){\bar{t}}_1(p_2){\delta }_{\chi }({\mathbf{p}}_2),$$

$$t_1(k)$$

$$t_2(k)$$

$$r_{12}^2(k)=\frac{⟨{\delta }_{\chi }{\delta }_{\chi }^{[2]}{⟩}^2}{⟨{\delta }_{\chi }{\delta }_{\chi }⟩⟨{\delta }_{\chi }^{[2]}{\delta }_{\chi }^{[2]}⟩}.$$

$${\bar{t}}_1(k)=\frac{⟨{\delta }_0{\delta }_{\chi }⟩}{⟨{\delta }_{\chi }{\delta }_{\chi }⟩},$$

$$t_2(k)=\frac{1}{1-r_{12}^2}\frac{⟨({\delta }_0-{\bar{t}}_1{\delta }_{\chi }),{\delta }_{\chi }^{[2]}⟩}{⟨{\delta }_{\chi }^{[2]}{\delta }_{\chi }^{[2]}⟩}.$$

$$\underset{k\to 0}{\lim }{t}_1(k)=\underset{k\to 0}{\lim }{\bar{t}}_1(k)=1$$

$$\underset{k\to 0}{\lim }{t}_2(k)=-\frac{3}{14}.$$

$$m(k)=e^{-(k\Sigma {)}^2/2}\left[O(k/\alpha )-1\right].$$

$${\sigma }^2(k)=\frac{2}{N_{\mathrm{modes}}(k)}⟨\frac{[{\hat{P}}_{\mathrm{wiggle}}(k){]}^2}{⟨{\hat{P}}_{\mathrm{nowiggle}}(k){⟩}^2}⟩$$

$$\bar{\rho }{d}^3\mathbf{q}=\bar{\rho }(1+\delta )d^3\mathbf{x},$$

$$A_{ij}\equiv \frac{\partial x_i}{\partial q^j}={\delta }_{ij}^K+{\psi }_{i,j}$$

$$\delta(\bm{x})=\Big{(}{1\over J_{q}({\bm{q}})}-1\Big{)}_{\bm{x}={\bm{q}}+\bm{\psi}}$$

$$\mathbf{x}=\hat{\mathbf{q}}+\mathbf{\chi }(\hat{\mathbf{q}}).$$

$$(1+{\delta }_{\hat{q}})d^3\hat{\mathbf{q}}=(1+\delta )d^3\mathbf{x},$$

$$(1+{\delta }_{\hat{q}})=[1+\delta (\mathbf{x}(\hat{\mathbf{q}}))]J_{\hat{q}}(\hat{\mathbf{q}}),$$

$$B_{ij}\equiv \frac{\partial x_i}{\partial {\hat{q}}^j}={\delta }_{ij}^K+{\chi }_{i,j}.$$

Full text

I Introduction

Studies of large scale structure (LSS) and the cosmic microwave background (CMB) have played a very important role in establishing the standard model of cosmology. Successful CMB experiments over the past few decades have managed to extract almost all of the available information in the primary anisotropies, and thus the importance of LSS studies is expected to increase in the future. Because of the sheer number of independent modes that can be measured in our observable Universe, in principle three-dimensional maps of the late Universe contain a vast amount of statistical information about the initial conditions for structure formation as well as cosmological parameters. This is especially true if one can extract information from small spatial scales.

The fluctuations we observe on small scales in the late Universe have been heavily processed by their dynamical evolution progressively scrambling the cosmological information to the extent that sufficiently small scales are usually excluded from cosmological analysis. The problem is made even worse by the fact that many of the remaining open questions in cosmology such as measuring the masses of neutrinos or constraining the Gaussianity of the initial seeds require measuring very small effects.

Thus there is great motivation for and a long history of trying to undo the cosmological evolution in a process usually called reconstruction. The hope is to unscramble the cosmological information that might still be present in higher order moments of the data and recover the initial conditions directly which can then be described using a power spectrum and perhaps a few low-order statistics. These techniques have been used most successfully in the context of baryonic acoustic oscillation (BAO) measurements. The BAO feature in the correlation function acts like a standard ruler whose physical size is calibrated by CMB measurements, so that measuring its angular extent on the sky as a function of redshift measures the expansion history of the Universe. However, the BAO feature is degraded (i.e. it gets broader) as a result of nonlinear evolution. The process of reconstruction sharpens back the peak, increasing the signal-to-noise ratio.

The operational procedure of reconstruction methods is usually motivated by analytical models that connect the linear initial density with the observed nonlinear density in the perturbative regime. In the past few years a lot of work has been put into trying to find improved models; for example a recent discussion of Lagrangian-space models can be found in Matsubara (2008, 2008, 2014); Porto et al. (2014); Carlson et al. (2013); Vlah et al. (2015); Seljak and Vlah (2015); Vlah et al. (2015, 2016a); McQuinn and White (2016); Baldauf et al. (2016a, b). Here we are motivated by developments connected with the effective field theory (EFT) approach to large-scale structure Baumann et al. (2012); Carrasco et al. (2012); Carroll et al. (2014); Senatore and Zaldarriaga (2015) and ask whether these theoretical modeling advances can be used to develop an improved reconstruction method. This is the goal of our paper, where we will focus on a Lagrangian EFT model developed in Baldauf et al. (2016b).

Related to the motivation for standard reconstruction Eisenstein et al. (2007a), the key motivation for our method is that the dominant nonlinear terms that reduce the correlation of the nonlinear density with initial conditions and degrade the BAO scale are nonlinear shift terms (rather than intrinsic nonlinearities like nonlinear growth or tidal terms). These shift terms can be reduced if we have a good estimate of the nonlinear displacement field from initial Lagrangian to final Eulerian coordinates. Motivated by Baldauf et al. (2016b), we estimate the nonlinear displacement field by iteratively moving particles back along estimated Zeldovich displacements. The resulting cumulative displacement field significantly improves the single-step Zeldovich displacement used in standard reconstruction. This is the most important aspect of our method. We refine the method further by adding second order corrections when estimating the linear density from the nonlinear displacement field, but this is less important for most applications. For simplicity we will only consider the idealized toy model of dark matter in real space in this paper. In future work we plan to extend the method to halos in redshift space, accounting also for realistic shot noise and survey selection functions.

The subject of reconstruction was pioneered in Refs. Peebles (1989, 1990), where the initial positions of eight Local Group galaxies were estimated by minimizing the action of trial orbits assuming negligible initial peculiar velocities, and in Ref. Bertschinger and Dekel (1989), where the assumption of irrotational Lagrangian velocities and the Zeldovich approximation Zel’dovich (1970) were used to estimate the linear velocity potential and linear mass density from observed radial peculiar velocities. In a different approach called Gaussianization, the nonlinear density was transformed in a rank-preserving way such that its one-point probability distribution function becomes Gaussian Weinberg (1992). In a method more similar to most current methods, including ours, the Zeldovich approximation was used to describe the time evolution of the velocity potential that can be integrated back in time without generating vorticity Nusser and Dekel (1992). Mass conservation can be imposed using the continuity equation Gramann (1993), and an iterative higher-order scheme related to this was proposed in Ref. Monaco and Efstathiou (1999). In the path interchange Zeldovich approximation (PIZA) method Croft and Gaztanaga (1997) the least action principle was applied assuming trajectories are straight lines, which can be combined with the Gaussianization method Narayanan and Croft (1999), or generalized beyond straight trajectories Nusser and Branchini (2000); Branchini et al. (2002). More recently, this was generalized to more realistic dynamics, finding that different methods perform similarly for velocity reconstruction if applied to realistic mock catalogs Keselman and Nusser (2016). In another generalization known as Monge-Ampere-Kantorovich reconstruction, an optimal mass assignment problem is solved Frisch et al. (2002); Brenier et al. (2003); Mohayaee et al. (2006); Sirko (2007); Mohayaee and Sobolevskiĭ (2008).

A related method now referred to as standard reconstruction Eisenstein et al. (2007a) first computes the Zeldovich displacement from the filtered nonlinear density (corresponding to the first step of the iterative procedure of our method), and then moves the clustered and a random catalog by that displacement to estimate the linear density from the density difference of the two displaced catalogs (this procedure removes the nonlinear shift terms that otherwise degrade the BAO scale, also see Padmanabhan et al. (2009); Tassev and Zaldarriaga (2012)). That method was successfully applied to real galaxy survey data from SDSS BOSS Padmanabhan et al. (2012); Anderson et al. (2012); Ross et al. (2015); Anderson et al. (2014); Beutler et al. (2016); Alam et al. (2016) and WiggleZ Kazin et al. (2014), improving BAO measurements typically by a factor of $\sim 2$. Similar improvements are expected for future galaxy redshift surveys including, for example, DESI Levi et al. (2013), Euclid Laureijs et al. (2011), and LSST LSST Science Collaboration et al. (2009).

A conceptually different approach to reconstruction is to sample the initial condition density, propagate it forward using an $N$-body or 2LPT code, compare against observations, resample modes that did not fare well, etc, which becomes computationally feasible using Hamiltonian sampling techniques Jasche and Wandelt (2013); Kitaura (2013). Another method to reconstruct the initial density from the observed density using the Hamiltonian MCMC technique was presented in Wang et al. (2013, 2009). While these methods fix the initial linear power spectrum and cosmology to constrain the initial modes, a converse approach Seljak et al. is to marginalize over the modes and constrain the initial power spectrum and cosmology using optimization techniques.

Following the success of the standard reconstruction method Eisenstein et al. (2007a) for improving measurements of the acoustic peak, a series of related reconstruction procedures were proposed recently to achieve further improvements or enable broader applications beyond BAO measurements Seo et al. (2010); Tassev and Zaldarriaga (2012); Seo and Hirata (2016); Obuljen et al. (2016); Schmittfull et al. (2015); Zhu et al. (2016a, b). These methods have a common structure: First, they estimate the displacement field from initial to final positions, and second, they estimate the linear density from that displacement. However, the concrete procedure for these two stages differs for the various proposed methods as summarized in Table 1. The practical application of standard reconstruction to real galaxy survey data and improved implementations are discussed, for example, in Padmanabhan et al. (2012); Burden et al. (2014); Achitouv and Blake (2015); Burden et al. (2015). In our method described below, we will modify both the first and second stage based on a concrete analytical model for structure formation.

An earlier similar attempt to improve the displacement field using iterative Zeldovich displacements and higher-order corrections Seo et al. (2010) found only negligible improvements over standard reconstruction. We attribute the main reason for this to the fact that the filtering smoothing scale was kept fixed in the iterative procedure, so that only shifts on very large scales were mitigated. Among other changes and differences in theoretical motivations, we improve that method by progressively reducing the smoothing scale with each iteration step (also see Appendix B).

Recently, a related reconstruction method was proposed in a series of papers by Zhu et al. Zhu et al. (2016a, b). This method is also based on improved estimates of the nonlinear displacement between initial Lagrangian and final Eulerian coordinates, finding substantial improvements over standard reconstruction in terms of the cross-correlation with the linear density Zhu et al. (2016a, b), Fisher information Pan et al. (2016), BAO Wang et al. (2017), as well as when applying it to halos Yu et al. (2017). In that approach, the displacement is estimated by writing the deformation of the observed density grid to an initially uniform grid under the continuity equation as a differential equation that is then solved numerically using a multigrid moving mesh algorithm. Where comparisons are possible, this method seems to perform very similarly to ours. Despite differences in the theoretical motivation and operational procedure for constructing the displacement field and for inferring the linear density from that displacement, the similarity in final performance is likely related to the fact that both methods rely on significantly improved displacement fields compared to standard reconstruction. We believe both approaches shall be useful in the future to improve BAO measurements and realize other applications of reconstruction.

Our paper is structured as follows. We start in Section II with the theoretical motivation for the particular reconstruction algorithm that we propose. Section III describes the implementation of the algorithm and the numerical setup and simulations. In Section IV we assess the performance of the reconstruction based on the one-point probability distribution function (pdf) of the reconstructed density, its cross-correlation with the linear initial conditions, the BAO signature in the power spectrum, and the full shape of the power spectrum. We conclude in Section V. Appendix A discusses transfer functions that are needed when including second-order corrections to the reconstruction. In Appendix B we describe a simple extension of standard reconstruction that performs worse than the method described in the main text, but still much better than standard reconstruction. Appendix C describes our procedure to fit the BAO signature in the power spectrum. Appendix D and Appendix E discuss results at higher redshift, the choice of reconstruction parameters, and simulation convergence tests. Appendix F provides a perturbative analysis motivating the reconstruction algorithm further.

II Theoretical motivation

Perturbation theory has been successful in modeling the large-scale statistical properties of the evolved LSS given the initial conditions. In the past few years the application of effective field theory (EFT) ideas has allowed additional improvements on scales larger than the nonlinear scale. Regardless of these improvements, perturbation theory breaks down on small scales once shell crossing has occurred. In addition, we also expect reconstruction algorithms to break in the same regime as it is not possible to uniquely infer the initial density field from the final density once shells have crossed. One could think of a shell of matter of some radius; given just the density one cannot know if the shell is collapsing for the first time or moving back out. Equivalently one could imagine finding a solution for the displacement of a uniform set of particles that after being displaced lead to an observed density field. Given any such solution, one can easily find another one by swapping particles around. Thus if perturbation theory and any reconstruction algorithm are expected to fail on the same scales, it seems relevant to try to use our better understanding of perturbation theory to improve reconstruction, as an algorithm based on perturbation theory might do as well as one possibly can do. We will take the first steps along this direction in this paper.

In particular, although BAO reconstruction has been very successful, we know it does not recover all the information that is there in the linear field. Thus, if we manage to improve the reconstruction algorithm, we might expect (or hope) that for sufficiently good data there might be some room for improvement in BAO related measurements. There is also the hope that reconstruction might be useful for studying other questions.

Let us start by being more precise about our theoretical expectations:

• •

Perturbation theory only works on large scales, so we should only try to use the large scale nonlinear density as an input to the reconstruction algorithm.

• •

In perturbation theory there are two sources of nonlinearities: shift terms and true changes to the small scale dynamics, for example corrections to the growth rates. These terms are of different sizes. We see this directly when we compare one-loop corrections to the displacement and density fields (the displacement is a Lagrangian quantity that is not corrected by the shift terms). The scale at which the one-loop correction to the displacement power is equal to linear displacement power is much smaller (at higher $k$) than that for which the one-loop corrections to the density power are equal to the linear density power (e.g. Baldauf et al. (2016a, b)). This shows that the biggest effect in the broadband error of linear theory is from shift terms, which affect the density but not the displacement. We discuss this more quantitatively in Appendix F.

• •

If we are interested in the BAO feature, we also know that it is the shift terms that damp the oscillations (broaden the peaks) Eisenstein et al. (2007b, a).

• •

We also know that the correlation function around the BAO scale computed in the Zeldovich approximation is almost the same as the fully nonlinear one. It is very difficult to spot the additional corrections from the other one-loop terms (i.e. the one-loop IR-resummed correlation function is difficult to tell apart from Zeldovich). Thus we expect that undoing Zeldovich is almost all that is necessary.

• •

If one assumes just one stream, there is a unique way to invert Zeldovich. But once shell crossing has occurred, no unique solution exists, and the problem is no longer invertible. So on small scales the displacement inferred by making the density uniform (inverting Zeldovich) cannot be related to the initial conditions in a simple form. Finding the correct displacement might not even be possible with simulations rather than perturbation theory. At best one can rely on a prior to rank the probability of different solutions.

These considerations suggest estimating the perturbative part of the density field by filtering the nonlinear density field, then inverting the relation between displacement and density to obtain the displacement that would transform a uniform density to the filtered nonlinear density, and then estimate the linear density from that displacement.

To see this more quantitatively, we can start from the Lagrangian model built in Baldauf et al. (2016b)

[TABLE]

where

[TABLE]

Here $\delta_{\mathrm{NL}}$ is the nonlinear final density and $\delta_{0}$ is the linear initial density. $a_{1}(k)$, $a_{2}(k)$ and $T(k)$ are transfer functions that only depend on the modulus of the Fourier wave vector. $\delta_{Z}[\delta_{\chi}]$ is the density obtained by moving a uniform set of particles by a displacement $\bm{\chi}$ whose divergence is $\delta_{\chi}\equiv{\bm{\nabla}}\cdot\bm{\chi}$. This model contains all the terms relevant for a one-loop calculation. In principle one should add a cubic term but up to one loop only the part that correlates with the linear density enters. Thus the cubic term is included in the $a_{1}$ transfer function. In the EFT approach one has to include counterterms, which at this order are described by the so-called speed of sound and related terms. These contributions are absorbed in $a_{1}$ and $T$. Finally, since we are working in Lagrangian space, all shift terms are included to all orders. Thus this scheme already does the so-called infrared resummation.

One should think of $\bm{\chi}$ as the nonlinear displacement from initial Lagrangian coordinates to late-time Eulerian coordinates, which can be modeled perturbatively in the initial (linear) density field $\delta_{0}$ as in Eq. (II), with ${\bm{k}}={\bm{p}}_{1}+{\bm{p}}_{2}$. The error density $\delta_{\mathrm{error}}$ accounts for the difference between the perturbatively based prediction and the simulation results. Its power spectrum was presented in Baldauf et al. (2016b). That paper showed that not much is gained by going to higher order in perturbation theory for $\delta_{\chi}$. Even with the quadratic based prediction for $\delta_{\chi}$, the cross-correlation coefficient between the perturbative and nonlinear fields drops to 0.5 only at $k\approx 0.5\;h\text{Mpc}^{-1}$ at redshift $z=0$.

Assuming the model of Eqs. (1) and (II) is the correct connection between linear density $\delta_{0}$ and nonlinear density $\delta_{\mathrm{NL}}$, how can we estimate $\delta_{0}$ from the observed $\delta_{\mathrm{NL}}$? Eq. (1) suggests that we first want to estimate the perturbative part of the density, $\delta_{PT}$, from the nonlinear density field. If all fields were Gaussian, that would involve a simple Wiener filter. We will follow the same procedure here and estimate $\delta_{PT}$ by filtering the nonlinear density field,

[TABLE]

where $W$ is a filter which we choose to be a Gaussian for simplicity. Then we estimate the divergence $\delta_{\chi}$ of the Lagrangian-to-Eulerian displacement by demanding

[TABLE]

This is a nonlinear equation for $\hat{\delta}_{\chi}$, which we will solve iteratively. The solution is an estimate $\hat{\bm{\chi}}$ of the Lagrangian-to-Eulerian displacement, and represents the first stage of the reconstruction. In a second stage, we estimate the linear density $\hat{\delta}_{0}$ corresponding to that displacement $\hat{\bm{\chi}}$. If $\hat{\bm{\chi}}$ contained only linear displacements, the linear density would simply be given by its divergence. If $\hat{\bm{\chi}}$ also contains nonlinear displacements, we can model the relation between nonlinear displacement and linear density using Eq. (II), and invert that to estimate the linear density from the nonlinear displacement using

[TABLE]

with transfer functions $t_{1}$ and $t_{2}$ calibrated using cross-correlations in simulations as described in Appendix A.

We still need to specify how we solve Eq. (4) for the Lagrangian-to-Eulerian displacement $\bm{\chi}$, or its divergence $\delta_{\chi}$. One could solve this equation perturbatively; we will present the relevant formulas in Appendix F. Here we will think of $\bm{\chi}$ as a ‘migration’ map that shows where each final dark matter particle originated. We can hope to estimate this because (a) dark matter particles do not migrate very far over the history of the Universe (typically a few comoving Mpc), (b) they preferentially travel along potential gradients that can be estimated, and (c) they originate from a nearly uniform initial distribution because the clustering power spectrum is small at early times. This suggests to compile the migration map $\bm{\chi}$ iteratively by estimating potential gradients from the observed late-time catalog of galaxies or matter overdensities, moving them back along those gradients, computing potential gradients again from the new catalog, moving them further back, and so on. The density of the displaced particles decreases from step to step, until it converges to a set of nearly uniformly distributed particles that will not be moved further because potential gradients nearly vanish. The total displacement from the late-time nonlinear catalog to the nearly uniform catalog is our estimate of the migration map or Lagrangian-to-Eulerian displacement $\bm{\chi}$.

An alternative motivation for this iterative procedure to solve Eq. (4) follows from the Eulerian continuity equation

[TABLE]

This can be written as a convective derivative

[TABLE]

We can solve this equation iteratively in small time steps, where in each step

[TABLE]

Here $\delta_{R}$ is the nonlinear density filtered with a progressively smaller smoothing scale $R$ that ensures $\delta_{R}\lesssim 1$ for each step.

Note that the solution of Eq. (4) is not unique after shell crossing. So we should not be too concerned with the specific scheme we use to solve the equation. Different choices will result in solutions that differ on small scales where they cannot be trusted.

III Numerical setup

Having motivated and described our reconstruction algorithm in the last section, we next describe a concrete numerical implementation of that algorithm and the numerical setup that we will use in Section IV to assess the performance of the reconstruction.

III.1 New reconstruction algorithm

III.1.1 Procedure

As described above, our reconstruction algorithm consists of two stages. In the first stage, we compute a displacement field $\bm{\chi}$ using an iterative scheme. This is an approximation of the true displacement between initial and final conditions. In the second stage, that displacement is used to estimate the linear density. In detail, we implement the new reconstruction using the following recipe.

Stage I: Iteratively solve for displacement field $\bm{\chi}$

Iteratively displace objects as follows:

1. (a)

From the catalog, compute the fractional overdensity, $\delta(\bm{x})=\rho(\bm{x})/\bar{\rho}-1$, on a regular grid using cloud-in-cell. 2. (b)

Apply Gaussian smoothing with kernel $W_{R}=e^{-(kR)^{2}/2}$, with smoothing scale

[TABLE]

in the $n$th iteration step, and truncate small-scale modes by setting $\delta({\bm{k}})=0$ if $k>k_{\mathrm{max}}$. 3. (c)

Compute the Zeldovich displacement $\bm{s}({\bm{k}})=-\epsilon_{s}\frac{i{\bm{k}}}{k^{2}}\delta({\bm{k}})$ on the grid. 4. (d)

For every object in the catalog, interpolate the displacement $\bm{s}(\bm{x})$ from the regular grid to the object’s position, move the object by that displacement, and update the object’s position in the catalog. 2. 2.

For every object in the catalog, compute its total displacement $\bm{\chi}(\hat{\bm{q}}_{\mathrm{end}})=\hat{\bm{q}}_{\mathrm{end}}-\bm{x}_{\mathrm{start}}$, where $\bm{x}_{\mathrm{start}}$ is the observed starting position before the first iteration step, and $\hat{\bm{q}}_{\mathrm{end}}$ is the end position of the object after iteratively displacing it ($\bm{x}_{\mathrm{start}}$ is the observed Eulerian position, and $\hat{\bm{q}}_{\mathrm{end}}$ is the estimated associated Lagrangian position). 3. 3.

Paint111At every grid point, we average the displacement $\bm{\chi}(\hat{\bm{q}}_{\mathrm{end}})$ over all nearby objects using the cloud-in-cell distance weight. If a grid cell is empty, which can happen if the grid is finer than the typical particle separation, we fill $\bm{\chi}$ with a random neighbor cell until all cells are filled. This avoids a broadband power deficit that would result from wrongly setting $\bm{\chi}=0$ in empty cells, but it does not affect the correlation coefficient or BAO much. the displacement $\bm{\chi}(\hat{\bm{q}}_{\mathrm{end}})$, defined at the estimated Lagrangian positions $\hat{\bm{q}}_{\mathrm{end}}$, to a regular grid, and truncate small-scale modes by setting $\bm{\chi}({\bm{k}})=0$ if $k>k_{\mathrm{max}}$.

Stage II: Given $\bm{\chi}$, estimate linear initial density

• •

$\mathcal{O}(1)$ reconstruction: Compute $\hat{\delta}_{0}^{[1]}=\delta_{\chi}={\bm{\nabla}}\cdot\bm{\chi}$.

• •

$\mathcal{O}(2)$ reconstruction: Using the transfer functions and definitions in Appendix A, compute

[TABLE]

We apply this procedure to dark matter (DM) particles in simulations in real space, leaving important extensions like galaxy bias and redshift space distortions for future work.

III.1.2 Implementation

The algorithm involves only simple operations such as moving objects by a displacement field or computing the divergence of the displacement, which are the same basic operations also used by standard reconstruction. Implementing our method is therefore not more difficult than standard reconstruction Eisenstein et al. (2007a). The computational cost is rather modest, only a few times higher than that of standard reconstruction if running with our default of eight iteration steps. We have implemented the reconstruction algorithm in a simple Python code. With this, reconstructing initial conditions from a catalog of $85\times 10^{6}$ particles using a $512^{3}$ grid takes only $3$ CPU hours, which is negligible compared to other computations that are required to analyze LSS data. The low computational cost and simple operational structure allow the algorithm to be applied to large catalogs, requiring only minimal modifications of codes that already implement the standard reconstruction of Eisenstein et al. (2007a).

For application to real data from galaxy redshift surveys, the algorithm would have to be adopted to deal, for example, with the survey selection function including possible gaps in the data, finite shot noise, and redshift space distortions. While we do not address these issues here, we note that Padmanabhan et al. (2012); Kazin et al. (2014) developed approaches to apply standard reconstruction to real data, and those should in principle also be applicable to our method. For example, the Zeldovich displacement in each iteration step could be computed by numerically solving the differential equation relating the displacement field and the density in configuration space rather than working on a grid and solving the equation in Fourier space.

III.1.3 Parameters

Our reconstruction algorithm has several parameters. Throughout the paper we use the following default parameters unless explicitly stated otherwise. We start the iteration with the smoothing scale $R_{\mathrm{init}}=10\ h^{-1}\text{Mpc}$. We then run $N_{\mathrm{steps}}=8$ iteration steps and halve the smoothing scale from one step to the next, $\epsilon_{R}=0.5$, until we reach the minimum smoothing scale $R_{\mathrm{min}}=1.01\,L/N_{\mathrm{grid}}$. In each iteration step we displace particles by the full Zeldovich displacement, $\epsilon_{s}=1$. To represent fields on a regular grid, we choose $N_{\mathrm{grid}}=512$ and truncate modes with $k>k_{\mathrm{max}}=2\pi/L\times N_{\mathrm{grid}}/2$.

We discuss these choices in Appendix E.1, noting that results are quite robust as long as the smoothing scale is progressively reduced from step to step, $\epsilon_{R}\sim 0.5$, and the initial smoothing scale is chosen relatively large, $R_{\mathrm{init}}\gtrsim 5\ h^{-1}\text{Mpc}$, to ensure $W_{R}\delta\lesssim 1$.

The $\mathcal{O}(2)$ reconstruction requires two transfer functions $t_{1}(k)$ and $t_{2}(k)$ that are smooth functions of the wave vector modulus $k$. We calibrate them using simulations as described in Appendix A. The $\mathcal{O}(1)$ reconstruction, in contrast, is just given by $\delta_{\chi}={\bm{\nabla}}\cdot\bm{\chi}$ and does not involve any transfer function.

For the case of dark matter in real space considered in this paper the reconstruction does not assume any cosmological model. This may change if the method is extended to biased tracers in redshift space.

III.2 Standard reconstruction

To compare methods, we also run the standard reconstruction method of Eisenstein et al. (2007a) on all our simulations. We implement this by making the following changes to the above recipe: Modify the first stage to run only one iteration step with $\epsilon_{s}=1$ and $R=10\ h^{-1}\text{Mpc}$; define $\tilde{\bm{\chi}}(\bm{x}_{\mathrm{start}})=\hat{\bm{q}}_{\mathrm{end}}-\bm{x}_{\mathrm{start}}$ as a function of the Eulerian starting position before displacing any objects; paint this $\tilde{\bm{\chi}}$ to a regular grid; generate a uniform catalog with particles on that regular grid and displace them by $\tilde{\bm{\chi}}$; call the density of that shifted uniform catalog $\delta_{s}$; also compute the density of the particles at $\hat{\bm{q}}_{\mathrm{end}}$ after they were displaced, and call the corresponding density $\delta_{d}$; and then $\hat{\delta}_{0}^{\mathrm{std}}=\delta_{d}-\delta_{s}$ is the reconstructed density.

A simple extension of this standard method is to use more than one iteration step for the displacement $\bm{\chi}$, similar to Seo et al. (2010). We discuss this in Appendix B.

III.3 Simulations

We use simulations produced with the FastPM particle-mesh code Feng et al. (2016). They closely resemble full $N$-body simulations at much lower computational cost. To study how far into the nonlinear regime reconstruction works, we run a small-volume simulation with box size $L=500\ h^{-1}\text{Mpc}$ and $120$ time steps linearly spaced between $a=0.01$ and $a=1$. To study the BAO signature that is only measurable on relatively large scales, $k\lesssim 0.5\;h\text{Mpc}^{-1}$, we use ten large-volume simulations from Ding *et al. *Ding et al. with $L=1380\ h^{-1}\text{Mpc}$ and $40$ time steps linearly spaced between $a=0.1$ and $a=1$. Each of the ten realizations was evolved from initial conditions with and without BAO signature in the initial power spectrum to allow for cosmic variance cancellation when studying the BAO Prada et al. (2016); Schmittfull et al. (2015); Ding et al. ; Wang et al. (2017). The ’nowiggle’ power spectrum without the BAO signature was matched to the broadband shape of the ’wiggle’ power spectrum with the BAO signature, using the method described in the appendix of Ref. Vlah et al. (2016b). The simulations evolve $2048^{3}$ particles by computing forces on a $4096^{3}$ particle-mesh grid. We take snapshots at redshift $z=0.6$, which is slightly higher than the peak redshift of the SDSS BOSS CMASS galaxy sample (e.g. Alam et al. (2016)), and today at redshift $z=0$. To speed up computations, we select a $1\%$ subsample that contains roughly $85\times 10^{6}$ DM particles. To study the BAO, we subsample the same particles in wiggle and nowiggle simulations by keeping every particle whose particle ID is divisible by 101 (a prime number that avoids undesired regularities in the selection).

The simulations assume a flat $\Lambda$CDM cosmology with $\Omega_{m}=0.3075,\Omega_{c}h^{2}=0.1188$, $\Omega_{b}h^{2}=0.0223$, $h=0.6774$ and $\sigma_{8}=0.8159$ based on Planck Collaboration et al. (2016). The fiducial BAO scale, given by the radius of the sound horizon at the drag epoch when photons decoupled from baryons, evaluates to $r_{\mathrm{BAO}}^{\mathrm{fid}}\approx 147.5\,\mathrm{Mpc}=99.9\ h^{-1}\text{Mpc}$, assuming the approximation of Eq. (16) in Aubourg et al. (2015) for simplicity.

IV Results

We first check the basic functionality of the iterative procedure in the first reconstruction stage that progressively displaces the DM particles along their potential gradients until a nearly uniform set of particles is obtained. The upper panel of Fig. 1 shows the power spectrum of the displaced catalog after $1,2,4,8$, and $16$ steps. As expected, the power decreases from step to step. After eight steps the large-scale power is reduced by a factor of $10^{8}$ at $z=0$. Using more than eight steps does not reduce the large-scale power much further, so that our default choice is to stop after eight iteration steps. On smaller scales, the power of the displaced catalog is still reduced, but by a much smaller factor, for example a factor of $\sim 100$ at $k=1\;h\text{Mpc}^{-1}$. In this regime, the displacement obtained from the iterative procedure is not unique because of shell crossing as discussed in Section II above. The lower panel of Fig. 1 shows how the power of the cumulative displacement $\bm{\chi}$ progressively increases with more iterative steps as the power of the displaced catalog decreases. On large scales, the power of the displacement divergence approaches that of the linear field, but it misses power on small scales where shell crossing happens. This shows that the basic concepts motivating the reconstruction algorithm in Section II work in the simulations.

To evaluate the numerical performance of our method in more detail we consider several metrics in the next sections, including density slices and histograms, the cross-correlation coefficient with the initial conditions, the BAO signature in the power spectrum, and the full broadband shape of the power spectrum.

IV.1 Qualitative behavior of the reconstructed density

Fig. 2 compares the reconstructed density with the initial conditions in two-dimensional slices through our $L=500\ h^{-1}\text{Mpc}$ simulation. The reconstruction successfully recovers the spatial structure of the initial conditions on large scales and in regions with low to moderate overdensities. On small scales and in regions with high late-time overdensity, however, the reconstruction is less successful at recovering the correct initial conditions because of shell crossing as discussed above.

Besides spatial structure, we can also consider the one-point pdf or histogram of the smoothed mass density as shown in Fig. 3. By construction, the linear initial conditions have a Gaussian pdf centered at $\delta=0$. Under nonlinear gravitational evolution the density turns into a rather non-Gaussian field with a skewed one-point pdf (see Bernardeau et al. (2002) for a review). Fig. 3 demonstrates that reconstruction transforms that nonlinear density back to a field with a one-point pdf that is much closer to a Gaussian again. In other words, our reconstruction Gaussianizes the mass density and moves cosmological information from the skewness and other higher-order moments back to the variance of the density, similar to findings in Schmittfull et al. (2015); Slepian et al. (2015); Hikage et al. (2017).

IV.2 Correlation with the linear initial conditions

The two-dimensional (2D) slices and one-point histograms of the mass density are useful qualitative measures, indicating that the reconstruction algorithm can at least partially recover the initial conditions. To check this more quantitatively and rigorously we turn to Fourier space.

Fig. 4 shows the cross-correlation coefficient $r(k)$ between the reconstructed density and the linear initial conditions as a function of Fourier wavenumber $k$ measured from our $L=500\ h^{-1}\text{Mpc}$ simulation at $z=0$. Using our first-order method, the reconstructed density is more than $95\%$ correlated with the initial conditions on scales $k\leq 0.31\;h\text{Mpc}^{-1}$. Using our second-order method, this is the case on scales $k\leq 0.35\;h\text{Mpc}^{-1}$. For comparison, the wavenumber where the correlation with the initial conditions drops below $95\%$ is $k=0.18\;h\text{Mpc}^{-1}$ for standard reconstruction, and $k=0.07\;h\text{Mpc}^{-1}$ for the nonlinear density without reconstruction, in the same simulation. Based on this correlation coefficient with the linear density, our reconstruction thus improves the $k$ range by a factor of $2$ over standard reconstruction, and by a factor of $5$ compared to performing no reconstruction. At redshift $z=0.6$ the improvement factors are similar; see Fig. 13 in the appendix.

Related to this, the lower panel of Fig. 4 shows the fractional error of the reconstructed density phases. This is represented by the power spectrum of the difference between reconstructed density and true linear density, in units of the linear power spectrum. This fractional nonlinear error power can be shown to reduce to one minus the squared correlation with the initial conditions Baldauf et al. (2016b).222Focusing on the density phase correlation, we have implicitly assumed here that the reconstructed density is rescaled by $\mathrm{TF}(k)\equiv\langle\hat{\delta}_{0}\delta_{0}\rangle/\langle\delta_{0}\delta_{0}\rangle$ as in Eq. 4.6 in Baldauf et al. (2016b). We will discuss the density amplitude later in Section IV.5 and Appendix A. Fig. 4 demonstrates that our reconstruction significantly reduces this nonlinear error and improves the correlation with the initial conditions on all scales, outperforming the standard method. For our second-order method, the nonlinear error power relative to the linear power is $1-r^{2}\simeq(10^{-6},10^{-4},10^{-3},10^{-2})$ at $k=(0.02,0.06,0.1,0.2)\;h\text{Mpc}^{-1}$ at $z=0$. At higher redshift, $z=0.6$, the nonlinear error at the same scales is slightly smaller, so that $1-r^{2}\simeq(10^{-6},10^{-4},10^{-3},10^{-2})$ at $k=(0.02,0.1,0.15,0.3)\;h\text{Mpc}^{-1}$ as shown in Fig. 13 in the appendix. At both redshifts the nonlinear error is very small, so that for most practical purposes we can regard the reconstructed and initial densities as identical on large scales.

At some point, reconstruction should be limited by the stochastic displacement term identified in Baldauf et al. (2016a), because it prohibits a deterministic mapping between initial and final conditions, at least using perturbation theory. This would imply that one cannot improve over $1-r^{2}\simeq(5\times 10^{-5},5\times 10^{-4},10^{-2})$ at $k=(0.06,0.1,0.2)\;h\text{Mpc}^{-1}$ at $z=0$, as shown in Fig. 22 of Baldauf et al. (2016a). Our second-order reconstruction reaches that limit within a factor of about 2. We therefore expect that other reconstruction methods could improve over our method by at most a factor of 2 on large scales.

As motivated in the Introduction, using our method to recover initial conditions over a wide range of scales can substantially improve many of the science goals of galaxy surveys by increasing the number of linear Fourier modes amenable to cosmological analysis. An important caveat is, however, that our numerical setup is rather idealized because we work with DM particles and ignore galaxy biasing and redshift space distortions. Both effects will certainly degrade the performance of reconstruction in practical applications. We plan to study this in future work.

As described in Section III, we use the eight-step displacement $\bm{\chi}^{(8)}$ for our new reconstruction, but the one-step displacement $\bm{\chi}^{(1)}$ for the standard reconstruction, because this is what has been used in the literature so far. Fig. 12 in Appendix D explores how the performance depends on the number of steps used to construct the displacement field $\bm{\chi}$. We see only little benefit in using more than eight iteration steps, indicating that the algorithm has converged after eight steps. Extending the standard reconstruction by applying it to the eight-step displacement $\bm{\chi}^{(8)}$ improves over using $\bm{\chi}^{(1)}$ on most scales, but it still performs worse than our second-order reconstruction (see Appendix B for discussion).

For the correlation coefficient shown in Fig. 4, the first order reconstruction does not depend on transfer functions because any rescaling by a function of $k$ would not affect the correlation coefficient. In contrast, the second order correction does require transfer functions that were calibrated to simulations as described in Appendix A. Just as in the case of forward modeling in Baldauf et al. (2016a), the shape of the transfer functions can probably be understood using the EFT approach, but we leave this for future work.

IV.3 Baryonic acoustic oscillations

Measurements of the BAO scale from the galaxy power spectrum are a prime example for the application of reconstruction, because it reverses or avoids large-scale shifts that would otherwise wash out the BAO wiggles in the observed galaxy power spectrum, degrading the measurement Eisenstein et al. (2007a, b); Padmanabhan et al. (2009). As mentioned above, the standard reconstruction technique has been successfully applied to several redshift surveys, improving the precision of the measured BAO scale typically by a factor of $\sim 2$ Padmanabhan et al. (2012); Anderson et al. (2012); Ross et al. (2015); Anderson et al. (2014); Kazin et al. (2014); Beutler et al. (2016); Alam et al. (2016), with similar improvements expected for future surveys. It is therefore exciting to see if our method can improve BAO measurements further. To answer this, we use ten large-volume simulations with $L=1380\ h^{-1}\text{Mpc}$ that were produced by Ding et al. Ding et al. as described in Section III.3 above.

Fig. 5 shows the fractional BAO signal in the simulations. Our method restores the BAO signal of the linear density perfectly, reversing the nonlinear damping. This is not surprising given that the BAO signal is only visible at $k<0.5\;h\text{Mpc}^{-1}$, where the reconstructed density is more than $90\%$ correlated with the linear density as we already found in Fig. 4. Standard reconstruction (green line in Fig. 5) also reduces the nonlinear damping, but it does not recover the full linear BAO wiggles at $k\gtrsim 0.2\;h\text{Mpc}^{-1}$.

To see if the signal-to-noise ratio of the BAO scale estimated from the power spectrum is also improved by reconstruction, we need to characterize the noise of the estimated BAO scale. This would be straightforward if we knew the covariance between power spectrum bins after reconstruction, but that is difficult to compute reliably. We therefore choose a simpler Monte Carlo approach and estimate the BAO uncertainty from the scatter of the best-fit BAO scale among the ten simulations. This provides a conservative estimate for the uncertainty of the best-fit BAO scale (see Appendix C, where we also describe our fitting procedure).

Fig. 6 compares the best-fit BAO scales estimated from linear initial conditions, nonlinear late-time density, and reconstructed density in each of the ten simulations, by fitting the BAO scale to the power spectrum at $k\leq k_{\mathrm{max}}=0.6\;h\text{Mpc}^{-1}$. This shows that our reconstruction recovers the linear BAO scale with high precision and on a realization-by-realization basis.

To estimate if the estimated BAO scale is systematically biased relative to the true BAO scale, we compute the expectation value of the best-fit BAO scale; see Table 2. Within the uncertainty of our ten simulations, we do not find evidence for any systematic BAO bias after any of the reconstruction methods that we tested. The reconstructions thus eliminate the systematic nonlinear BAO bias of $\sim 0.3\%$ at $z=0$ that is generated by shifts of particles that were separated by the pristine BAO scale in the initial conditions Eisenstein et al. (2007b, a), and that would be present when measuring the BAO scale from the nonlinear power spectrum without reconstruction. This is consistent with previous findings Seo et al. (2010); Padmanabhan et al. (2009); Noh et al. (2009); Sherwin and Zaldarriaga (2012); Tassev and Zaldarriaga (2012); White (2015); Mehta et al. (2011).

To estimate the statistical $1\sigma$ uncertainty corresponding to measuring the BAO scale from the power spectrum in a $2.6\,h^{-3}\mathrm{Gpc}^{3}$ volume, we compute the root-mean-square (rms) scatter of the best-fit BAO scale between the ten simulations; see Table 3 and Fig. 7.

The uncertainty of the BAO scale from the nonlinear power spectrum is increased by a factor of 2.8 at $z=0$ and by a factor of 2.6 at $z=0.6$ relative to the uncertainty from the linear power spectrum. This is again caused by shifts of particles that were separated by the BAO scale in the early universe. By reducing those shifts, standard reconstruction Eisenstein et al. (2007a) reduces the statistical BAO uncertainty by a factor of 1.6 at $z=0$ and by a factor of 1.9 at $z=0.6$ relative to performing no reconstruction. Standard reconstruction thus reduces the nonlinear degradation of the BAO signature, but it does not reverse it entirely. In contrast, our iterative second-order reconstruction reduces the BAO uncertainty by a factor of 2.7 at $z=0$ and by a factor of 2.5 at $z=0.6$ relative to performing no reconstruction, i.e. the nonlinear degradation is fully reversed within the uncertainty of our estimates. Relative to the standard reconstruction, the improvement from our method ranges from $70\%$ at $z=0$ to $30\%$ at $z=0.6$.

At both redshifts, $z=0.6$ and $z=0$, the statistical BAO uncertainty after the iterative second order reconstruction matches that of the linear initial conditions ($\simeq 0.24\%$), so that our reconstruction cannot be further improved for measuring BAO in an observed volume of $2.6\,h^{-3}\mathrm{Gpc}^{3}$. We will return to this point in the next section. Again, an important caveat is that practical issues like galaxy bias, shot noise, and redshift space distortions are expected to degrade the reconstruction performance in practice.

IV.4 Impact of residual nonlinear shifts on BAO

The measured properties of the best-fit BAO scale in the last section have a non-negligible uncertainty because they were obtained from a somewhat small number of only ten realizations. To provide a more stringent test, we cancel the linear cosmic variance contribution to the BAO uncertainty by measuring the difference of the best-fit BAO scale between initial and reconstructed density in each individual realization.

To motivate this procedure further, it is useful to split the statistical uncertainty of the measured BAO scale into a linear cosmic variance contribution $\sigma_{\mathrm{lin}}$ that would be present even if we were given the linear density (due to the finite observed volume and finite width of the linear acoustic peak) and an additional nonlinear noise term $\sigma_{\mathrm{NL}}$ determined by the rms of nonlinear shifts broadening the acoustic peak,

[TABLE]

We included a prefactor to indicate the expected volume scaling. Reconstruction reduces $\sigma_{\mathrm{NL}}$ by reducing nonlinear shifts. At some point the residual shift terms are small enough that their broadening does not matter compared to the intrinsic width of the acoustic peak in the linear density. Once we reach that threshold, reconstruction cannot improve BAO measurements any further Eisenstein et al. (2007a). In other words, the BAO signal-to-noise ratio can never be better than if we observed the linear density in the same volume.

After our reconstruction, we cannot distinguish the total BAO uncertainty from the linear uncertainty, i.e. the nonlinear noise term $\sigma_{\mathrm{NL}}$ is reduced so much that it is not detectable in our simulations. Our reconstruction thus reduces the rms of nonlinear shifts well below the threshold corresponding to the intrinsic width of the linear acoustic peak.

To see how well different reconstruction methods suppress the nonlinear shift noise $\sigma_{\mathrm{NL}}$ we cancel the linear cosmic variance noise term by subtracting the linear BAO scale of the initial conditions of each simulation, as illustrated in Fig. 8. This again shows that our reconstruction recovers the correct BAO scale of the initial conditions realization by realization, in a much better way than standard reconstruction. An estimate for the ratio of nonlinear noise relative to the BAO signal, $\sigma_{\mathrm{NL}}/r_{\mathrm{BAO}}$, is shown in Fig. 9 (also see right column of Table 3 for reference). This demonstrates that the nonlinear noise caused by nonlinear shift terms is substantially reduced by reconstruction. Relative to the nonlinear density, standard reconstruction reduces the nonlinear noise by a factor of 2.2, our first-order method reduces it by a factor of 9, and our second-order method reduces it by a factor of 15. We obtain similar factors at redshift $z=0.6$.

Our reconstruction method thus reduces the rms dispersion of nonlinear shifts much more effectively than standard reconstruction. However, as alluded to above, most of that improvement is not relevant in practice because the acoustic peak of the linear density has some intrinsic width so that reducing the rms of the nonlinear shifts below that width does not help to improve the BAO signal-to-noise ratio further. Our reconstruction reduces the nonlinear shifts well below that threshold, easily obtaining the BAO signal-to-noise ratio of the linear density.

The fact that our reconstruction works so much better than it needs to makes us optimistic that it can also restore the linear acoustic peak in more realistic scenarios that account for galaxy bias, redshift space distortions, or higher shot noise. The efficient suppression of nonlinear shift terms also explains the excellent correlation of the reconstructed density with the linear density. This is useful for applications beyond BAO measurements, for example to extract cosmological information from the broadband power spectrum as we will discuss in the next section, or information about the statistics of primordial fluctuations from higher-order statistics of the reconstructed density.

In summary, the above results show that our iterative reconstruction method is very powerful at reconstructing the BAO scale, and substantially improves over the standard reconstruction method. The reconstructed BAO scale is so close to the true BAO scale of the initial conditions that the two cannot be told apart within cosmic variance of current and near-future survey volumes, at least according to our idealized simulations.

IV.5 Full shape of the reconstructed power spectrum

Having discussed the correlation of the reconstructed density with the initial conditions and the BAO signature in the power spectrum, we can also ask how well our reconstruction method recovers the full broadband power spectrum of the initial conditions. This is shown in Fig. 10, where the measured power spectrum after reconstruction is divided by the linear initial power spectrum linearly scaled to the redshift $z=0$ of the simulation.

The original nonlinear power spectrum agrees with the linear power spectrum within $5\%$ at $k\leq 0.11\;h\text{Mpc}^{-1}$ in our simulations at redshift $z=0$. Our first-order reconstruction without transfer functions improves this slightly, so that reconstructed and linear power agree within $5\%$ at $k\leq 0.16\;h\text{Mpc}^{-1}$. The second-order method gives a larger improvement, agreeing with the linear power within $5\%$ at $k\leq 0.28\;h\text{Mpc}^{-1}$. At higher redshift, $z=0.6$, the power spectra agree with the linear one down to smaller scales, with $5\%$ agreement at $k\leq 0.15\;h\text{Mpc}^{-1}$ for the nonlinear density, $k\leq 0.33\;h\text{Mpc}^{-1}$ for first-order reconstruction, and $k\leq 0.43\;h\text{Mpc}^{-1}$ for second-order reconstruction (see Fig. 14 in the appendix).

In our idealized simulations, our reconstruction method thus extends the $k$ range where linear theory is valid by a factor of 2 to 3. In principle, the reconstruction could thus enable cosmological analyses to include the power spectrum on moderately nonlinear scales, $0.15\;h\text{Mpc}^{-1}\lesssim k\lesssim 0.4\;h\text{Mpc}^{-1}$ at $z\simeq 0-0.6$, where linear theory would not be valid without reconstruction. Since the precision of many cosmological constraints scales steeply with the number of Fourier modes included in the analysis, this could substantially improve cosmology constraints from present and future galaxy surveys. For this to be possible in practice, however, one would have to address several complications, in particular nonlinear halo bias and redshift space distortions that affect the broadband shape of the power spectrum.

On small scales, $k\gtrsim 0.3\;h\text{Mpc}^{-1}$ at redshift $z=0$, the reconstruction suffers from a deficit in power relative to the linear density. Qualitatively, this is not surprising because the reconstructed displacement field does not contain contributions from modes that underwent shell crossing, so that particles that ended up in halos are moved back farther than in the true Universe, leading to a suppression of power on small scales (see Fig. 2). This is the same reason for why the power spectrum in the Zeldovich approximation misses small-scale power: Particles free-stream through halos rather than undergoing nonlinear collapse to form halos, which washes out small-scale clustering power that is present in the true Universe.

More quantitatively, we show in Appendix F that reconstruction fully removes nonlinear shift terms, while it only partially removes nonlinear growth and tidal terms. While we attempt to mitigate the latter nonlinearities at second order by using the second-order reconstruction scheme, we have not attempted to remove nonshift nonlinearities at third order or higher. The residual nonlinearities likely suppress the power spectrum on small scales. We therefore expect the reconstructed power spectrum to agree with the linear power spectrum only on large scales, with a suppression relative to linear on smaller scales as seen in Fig. 10. On smaller scales we would have to add corrections to the linear power spectrum to match the reconstructed power spectrum, in particular 1-3 and $k^{2}P_{\mathrm{lin}}$ counterterms. These can be computed perturbatively, and one might hope that these corrections are easier to model than the fully nonlinear power spectrum, for example because the corrections contain no BAO wiggles, but we leave this for future work. Instead we account for the small-scale mismatch of linear and reconstructed power spectra using transfer functions $t_{1}(k)$ and $t_{2}(k)$ calibrated to the simulations as described in Appendix A.

The fact that the reconstruction misses small-scale broadband power relative to linear is not necessarily a problem. All that is needed to estimate cosmological parameters is a theoretical model predicting the reconstructed power spectrum as a function of cosmological parameters. Such a model could be developed using e.g. the Zeldovich approximation or higher-order LPT models. In Appendix F we provide a starting point for this, and we plan to compare this against the reconstructed density in simulations in future work (noting that 1-3 terms and a $k^{2}P_{\mathrm{lin}}$ term should also be added as mentioned above). For standard reconstruction several models have already been developed to describe the reconstructed power spectrum Padmanabhan et al. (2009); Sherwin and Zaldarriaga (2012); Schmittfull et al. (2015); White (2015); Baldauf et al. (2015); Cohn et al. (2016); Hikage et al. (2017), and it would be interesting to extend these to our method. To answer the question of how much cosmological information can be obtained from the full broadband shape of the reconstructed power spectrum, one additionally needs to determine the covariance between power spectrum bins after reconstruction Pan et al. (2016).

As a slightly technical remark, we note that it is important to fill empty $\bm{\chi}$ cells with random neighbor grid cells as described in Section III.1 to avoid a large-scale bias of the power spectrum that would be caused by setting $\bm{\chi}$ to zero in empty cells (which would effectively violate mass conservation). This is mainly important for the full-shape power spectrum and does not affect the correlation coefficient with initial conditions or BAO much.

V Conclusions

We have presented a new method to reconstruct the initial conditions in a given volume based on recent theoretical progress relating the linear initial density with the final nonlinear large-scale structure density Baldauf et al. (2016b). We first filter the density to suppress highly nonlinear scales where perturbation theory is not valid. We then estimate the displacement field that relates the uniform initial density in Lagrangian space to the filtered nonlinear density in Eulerian space by iteratively moving particles back with Zeldovich displacements. In this iterative procedure, we progressively reduce the smoothing scale to capture progressively smaller scales while keeping $\delta\lesssim 1$. We then estimate the linear initial density as the divergence of the cumulative Zeldovich displacement. This estimate can be improved further by adding a second-order correction.

In dark matter simulations, we find that after eight iteration steps the reconstructed density is more than $95\%$ correlated with the linear initial density of the simulation for scales $k\leq 0.35\;h\text{Mpc}^{-1}$ at redshift $z=0$, and for $k\leq 0.53\;h\text{Mpc}^{-1}$ at $z=0.6$. The power in the error between initial and reconstructed fields in units of the power spectrum of the initial field, given by one minus the squared correlation coefficient between initial and reconstructed fields, is very small, approximately $1-r^{2}=(10^{-6},10^{-4},10^{-2})$ at $k=(0.02,0.06,0.2)\;h\text{Mpc}^{-1}$ at redshift $z=0$. Our method therefore serves as an excellent estimator for the linear density.

As a concrete application, we have demonstrated that the method improves the standard reconstruction method Eisenstein et al. (2007a) to restore the initial linear-theory BAO signature from the final nonlinear density contrast of the Universe. Compared to performing no reconstruction, our method improves the BAO signal-to-noise ratio by a factor of 2.7 at redshift $z=0$ and by a factor of 2.5 at $z=0.6$, matching the optimal BAO signal-to-noise ratio of the linear density in our simulated volume for both redshifts. This improves over standard reconstruction Eisenstein et al. (2007a) by $70\%$ at $z=0$ and $30\%$ at $z=0.6$, using the same idealized simulations.

This demonstrates that our method reduces large-scale flows to a level where their effect on the width of the BAO peak is much smaller than the intrinsic width of the linear BAO peak imprinted at recombination Eisenstein et al. (2007b), i.e. the iterative reconstruction works much better than it needs to for BAO measurements. We confirmed this by showing that the reconstructed BAO scale agrees with that of the linear initial conditions in each of our simulated realizations even if cosmic variance is canceled, at least under the idealized assumptions of our simulations. Given this overly good performance in idealized situations, the iterative reconstruction might still recover the optimal linear BAO signal-to-noise ratio in more realistic situations. This should make our method interesting for future galaxy surveys like DESI Levi et al. (2013), Euclid Laureijs et al. (2011), and LSST LSST Science Collaboration et al. (2009) that measure the BAO scale to map the expansion history of the Universe. Similarly to other reconstruction methods, the improvement is larger toward lower redshifts where nonlinearities degrade the BAO signature more.

As a more general application, we have shown that the reconstruction also improves the agreement of the full broadband shape of the power spectrum with the linear power spectrum on moderately nonlinear scales. The reconstructed and linear power spectra agree within $5\%$ for $k\lesssim 0.16\;h\text{Mpc}^{-1}$ at $z=0$ and $k\lesssim 0.3\;h\text{Mpc}^{-1}$ at $z=0.6$ if we use our first-order reconstruction method without any transfer functions. This improves to $k\simeq 0.3\;h\text{Mpc}^{-1}$ at $z=0$ and $k\simeq 0.4\;h\text{Mpc}^{-1}$ at $z=0.6$ if we use second-order reconstruction that involves transfer functions calibrated to simulations. This improves the $k$ range where linear theory is valid by a factor of 2-3 relative to performing no reconstruction. If a similar improvement can be obtained for real data, this could have substantial impact on cosmological constraints from galaxy power spectra because the statistical uncertainty of the power spectrum drops rapidly with increasing wavenumber. Additional gains may be possible by modeling the shape of the reconstructed power spectrum perturbatively or by improving the reconstruction method further, but we leave this for future work.

We demonstrated that modeling the relation between the nonlinear displacement and linear density to second order and inverting it improves the correlation of the reconstruction and the initial conditions substantially on large scales. For the BAO application however, most of the improvement compared to previous methods is due to the improved iterative displacement field constructed in the first stage of the method. Indeed, for BAO measurements one could also use a simple extension of the standard reconstruction method by displacing clustered and random catalogs by the improved iterative displacement as discussed in Appendix B.

In this paper, we have only studied the idealized toy model of dark matter in real space, ignoring important practical issues related to halo biasing, shot noise, redshift space distortions, and survey selection function. These effects are expected to degrade the effectiveness of any reconstruction method because they add stochasticity and are difficult to model. We plan to study the quantitative impact of these effects on the reconstruction performance in future work.

Related to this, it is not entirely clear how to adopt our method to deal with real data that can have gaps and other complications. We note, however, that the method involves only simple operations such as estimating Zeldovich displacements from the density and moving objects, and these are the same operations as for standard reconstruction Eisenstein et al. (2007a). To apply our method to real data one could therefore try to follow the same approaches that were used to apply standard reconstruction to real data Padmanabhan et al. (2012); Kazin et al. (2014). Using our default of eight iteration steps, our method is only a few times slower than standard reconstruction, which should be acceptable given the potential gains.

A method similar to ours has recently been proposed by Zhu et al. Zhu et al. (2016a, b) and was studied further in Pan et al. (2016); Wang et al. (2017); Yu et al. (2017), finding similar results to ours where comparisons are possible. Their method differs from ours in mainly two aspects: (1) The displacement between initial and final conditions is estimated using a multigrid algorithm to smoothly deform coordinates (or solve the corresponding differential equation in small time steps), whereas we iteratively displace particles by their Zeldovich displacement. (2) To estimate the linear density from the displacement, both methods first compute the divergence of the displacement, but we also add a second-order correction that improves the correlation with initial conditions. The similar performance of the methods is likely due to the fact that the displacement field obtained by both methods improves the one-step Zeldovich displacement used by standard reconstruction, and in the regime where shell crossing is relevant, the details of the method to find the displacement likely do not matter much because the solution is not unique in that regime. Both reconstruction methods do not assume a cosmological model. They shall both be useful in the future because they likely enable similar gains but with different systematics because of the different operations involved.

We have demonstrated that a simple iterative reconstruction method improves significantly over the standard method in idealized simulations. If similar improvements can be achieved in more realistic scenarios including noise, halo biasing and redshift space distortions, this could improve BAO measurements and increase the number of power spectrum Fourier modes amenable for cosmological analysis for future galaxy surveys like DESI Levi et al. (2013), Euclid Laureijs et al. (2011), and LSST LSST Science Collaboration et al. (2009). This could significantly improve the scientific return of these surveys, including high-precision probes of the expansion history of the Universe, dark energy, neutrino mass, and the statistics of the primordial fluctuations.

ACKNOWLEDGEMENTS

We thank Zhejie Ding for sharing a suite of large-volume simulations optimized for BAO studies Ding et al. . We also thank Yu Feng for making his code base publicly available and for providing extensive help for how to use it. We are also grateful to Zvonimir Vlah for sharing linear power spectra with and without BAO wiggles (see appendix of Ref. Vlah et al. (2016b)). We acknowledge useful comments from Anže Slosar, Nikhil Padmanabhan, Uroš Seljak, and Marko Simonović. Parts of our reconstruction algorithm use nbodykit, a publicly available333https://github.com/bccp/nbodykit parallel Python postprocessing package developed by Yu Feng and Nick Hand. Our simulations were run with the publicly available444https://github.com/rainwoodman/fastpm FastPM code developed by Yu Feng Feng et al. (2016). This research used resources of the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We also acknowledge use of the Hyperion cluster at the Institute for Advanced Study. M.S. acknowledges support from the Bezos Fund through the Institute for Advanced Study.

Appendix A Transfer functions

To obtain the transfer functions needed for the second-order reconstruction, we write the reconstructed density as

[TABLE]

where $\delta_{\chi}={\bm{\nabla}}\cdot\bm{\chi}$ and

[TABLE]

where ${\bm{p}}_{2}={\bm{k}}-{\bm{p}}_{1}$, $\bar{t}_{1}$ is a filter to be specified later, and $\kappa_{2}$ is defined in Eq. (42) below. We choose the transfer functions $t_{1}$ and $t_{2}$ in Eq. (12) to minimize the difference $\langle(\hat{\delta}_{0}-\delta_{0})^{2}\rangle$ with the linear density in simulations. This gives

[TABLE]

and

[TABLE]

where we assumed $\bar{t}_{1}$ to be fixed, and we used the correlation coefficient of the first and second order contribution,

[TABLE]

To specify the transfer function $\bar{t}_{1}$ inside the second order contribution in Eq. (13) we minimize $\langle(\bar{t}_{1}\delta_{\chi}-\delta_{0})^{2}\rangle$, i.e. the error in the absence of a second order contribution. This gives

[TABLE]

Indeed, we have $\bar{t}_{1}=t_{1}$ if $\langle\delta_{\chi}\delta_{\chi}^{[2]}\rangle$ were to vanish. On large scales $\delta_{\chi}$ agrees with $\delta_{0}$ so that $\bar{t}_{1}=1$. On smaller scales $\delta_{\chi}$ differs from $\delta_{0}$ so that $\bar{t}_{1}$ approaches zero, which down-weights small-scale modes in Eq. (13). Using $\bar{t}_{1}$, we can rewrite $t_{2}$ in a form reminiscent of Gram-Schmidt orthogonalization,

[TABLE]

In the above equations, the shorthand notation $\left\langle\delta_{0}\delta_{\chi}\right\rangle\equiv P_{\delta_{0},\delta_{\chi}}(k)$ refers to the cross spectrum between $\delta_{0}$ and $\delta_{\chi}$ as a function of wavenumber $k$, and analogously for other fields. All transfer functions depend only on the modulus $k$ of wave vectors and can therefore be calibrated from simulations.

To calibrate the transfer functions, we measure the auto- and cross-spectra of $\delta_{0}$, $\delta_{\chi}$ and $\delta_{\chi}^{[2]}$ from the simulations, compute the transfer functions, and average over realizations where more than one realization is available. The result is shown in Fig. 11.

At low $k$, the transfer functions approach their theoretical low-$k$ limits given by

[TABLE]

and

[TABLE]

These limits follow by using the expansion (F.5) of $\delta_{\chi}$ in $\delta_{0}$, modeling all spectra at leading order in $P_{\mathrm{lin}}$, and using the large-scale limit given in Eq. (54). At higher $k$, the transfer functions deviate from that limit, but they are still rather smooth functions of $k$. It should therefore be possible to either model them analytically or continue the practice of calibrating them against simulations. Ultimately, the transfer functions should also include halo bias and redshift space distortions.

Appendix B Extended standard reconstruction

In this appendix we discuss a simple extension of the standard reconstruction method of Eisenstein et al. (2007a) that could be used as an alternative to the reconstruction method discussed in the main text. We regard the method in the main text as superior on theoretical grounds because the conversion from nonlinear displacement to linear density is better motivated, but the alternative method discussed here could in principle be useful for practical purposes (for example, it might be easier to deal with the survey selection function).

B.1 Method

We first apply the iterative procedure described in Section III.1 to find the displacement field $\bm{\chi}^{(8)}$ with eight iteration steps, which improves the Zeldovich displacement used in standard reconstruction. This is the same as the first stage of the method described in the main text. We can now ask how well we can do if we apply the second stage of standard reconstruction to the improved displacement field. Following Section III.2, we thus displace the clustered catalog and a uniform catalog by $\bm{\chi}^{(8)}$ defined in Eulerian coordinates $\bm{x}$ and then take the density difference. This second stage differs from the method in the main text, where we take the divergence of $\bm{\chi}^{(8)}$ defined in approximate Lagrangian coordinates $\hat{\bm{q}}$ with an optional second order correction; also see Table 1.

This extension of the standard method based on an iterative displacement field is very similar to the iterative standard reconstruction of Ref. Seo et al. (2010), where no improvement over standard reconstruction was found. However, an important difference is that the smoothing scale used to compute the iterative Zeldovich displacements was held fixed in Ref. Seo et al. (2010), whereas we progressively decrease the smoothing scale from step to step to access progressively smaller scales. We believe that this is the main difference compared to the implementation in Seo et al. (2010), but there may be additional differences. As shown in the lowest row in Table 6, we also find little improvement from iterations if we keep the smoothing scale fixed, consistent with the findings of Seo et al. (2010), but, as we discuss next, we find substantial improvements when reducing the smoothing scale from step to step, $\epsilon_{R}<1$.

B.2 Performance and discussion

The right panel of Fig. 12 shows that after applying the extended standard reconstruction described above, the large-scale density is very correlated with the linear initial conditions. The achieved correlation is similar to the $\mathcal{O}(1)$ reconstruction discussed in the main text and shown in the left panel of Fig. 12, but worse than the $\mathcal{O}(2)$ reconstruction shown in the middle panel of Fig. 12.

This demonstrates that the second order scheme of Eq. (II) is more successful at converting a given nonlinear displacement field to the linear density than just taking the divergence of that displacement or displacing clustered and random catalogs and taking their density difference. This is expected given that the second order scheme inverts the model of Eq. (II) which has been validated against simulations Baldauf et al. (2016b).

For measuring the BAO scale, reconstruction does not need to work perfectly as discussed in Eisenstein et al. (2007a) and Section IV.4. Indeed, Tables 4 and 5 show that the extended standard reconstruction works sufficiently well that it fully recovers the linear BAO signal-to-noise ratio, improving significantly over the standard one-step reconstruction. For BAO measurements, the $\mathcal{O}(1)$ and $\mathcal{O}(2)$ methods from the main text and the extended standard method discussed in this appendix thus all perform equally well, improving the standard method based on a one-step displacement.

This result shows that the main improvement compared to the standard one-step method of Eisenstein et al. (2007a) comes from the improved iterative displacement field. To obtain this, it is crucial to progressively reduce the smoothing scale in the iterations to access smaller and smaller scales with each iteration step (see Table 6). Once this improved displacement field is obtained in the first stage of reconstruction, it does not matter for BAO measurements how the displacement is converted to an estimate for the linear density in the second stage. For BAO measurements from real data one should therefore choose the second stage based on practical considerations and characteristics of the given data set. For applications beyond BAO measurements, the second order reconstruction in the main text should be used instead because it achieves the highest correlation with the initial conditions and best broadband power spectrum shape.

Appendix C Fitting BAO from the power spectrum

We employ the following procedure to estimate the BAO scale from the power spectra measured in our simulations before and after reconstruction. In each realization, we compute the fractional BAO signal of the power spectrum, $\hat{d}(k)=\hat{P}_{\mathrm{wiggle}}(k)/\langle\hat{P}_{\mathrm{nowiggle}}(k)\rangle-1$. Here, $\hat{P}_{\mathrm{wiggle}}$ is the power spectrum measured in a simulation realization with BAO wiggles; $\hat{P}_{\mathrm{nowiggle}}$ is the same but in a simulation initialized with no BAO wiggles; and $\langle\cdot\rangle$ represents the average over our ten realizations. We then fit that data in each realization with a simple theoretical model that multiplies the BAO scale by a factor $\alpha$ and suppresses the BAO wiggles using a damping factor $\Sigma$,

[TABLE]

Here $O(k)=P^{\mathrm{lin}}_{\mathrm{wiggle}}(k)/P^{\mathrm{lin}}_{\mathrm{nowiggle}}(k)$ is the ratio of the theoretical linear wiggle and nowiggle power spectra that were used to initialize the simulations. We fit for $\alpha$ and $\Sigma$ by minimizing the chi-squared $\sum_{k}[\hat{d}(k)-m(k)]^{2}/\sigma^{2}(k)$ at $0<k\leq k_{\mathrm{max}}^{\mathrm{fit}}=0.6\;h\text{Mpc}^{-1}$ for every realization.555We allow $\Sigma$ to vary between [math] and $20\ h^{-1}\text{Mpc}$, and $\alpha$ between $0.9$ to $1.1$. This is conservative given that $\alpha$ should deviate by less than $1\%$ from unity in any realization of the large-volume simulations. Best-fit values are consistent when fitting with scipy.optimize.minimize’s Nelder-Mead, Powell or TNC algorithms.

This gives the best-fit BAO scale in each realization, $\hat{r}_{\mathrm{BAO}}=\hat{\alpha}\,r_{\mathrm{BAO}}^{\mathrm{fid}}$, as shown in Fig. 6. We estimate the uncertainty of the measured BAO scale by computing the scatter of the best-fit BAO scale between the realizations. This Monte Carlo method to estimate the BAO uncertainty should provide a robust estimate of the true uncertainty because it quantifies how much the BAO scale estimated from our particular fitting procedure varies among different realizations of the Universe.

When fitting for the BAO scale, we assume a Gaussian diagonal covariance $\mathrm{cov}(\hat{d}(k),\hat{d}(k^{\prime}))=\delta_{kk^{\prime}}\sigma^{2}(k)$ with

[TABLE]

estimated from the simulations. This should be regarded just as a particular weight in the fitting process to up-weight the power spectrum on small scales where there are more modes. If the assumed Gaussian covariance is wrong, the estimator for the BAO scale is suboptimal because it employs a suboptimal weight. In that case, the uncertainty estimated from the scatter of the best-fit BAO scale between realizations overestimates the uncertainty compared to a more optimal estimator or fitting procedure based on the correct covariance. As a consequence, our estimated BAO uncertainty is a conservative estimate of the true uncertainty.

Appendix D Results at redshift $z=0.6$

Figures 13 and 14 show reconstruction results at redshift $z=0.6$ as opposed to the redshift $z=0$ that was used in the main text. The densities before and after reconstruction match the linear density better at this higher redshift than at lower redshift, which is as expected because nonlinearities are smaller at higher redshift for any given scale.

Appendix E Parameters and convergence tests

In this appendix we discuss some choices we made for the reconstruction parameters, and some basic convergence tests of our simulations.

E.1 Reconstruction parameters

Our reconstruction algorithm has several parameters as described in Section III.1.3. Table 6 shows results for different parameter choices. This demonstrates that the final performance of the method is relatively insensitive to the detailed parameter values. Some qualitative choices are important though: It is important to start with a relatively large smoothing scale, $R_{\mathrm{init}}\gtrsim 5\ h^{-1}\text{Mpc}$, so that the smoothed overdensity is less than unity for most grid points and the Zeldovich approximation is applicable. The smoothing scale should also decrease from one iteration to the next, $\epsilon_{R}<1$, to reconstruct progressively smaller scales in the iteration. However, since reconstruction likely becomes inefficient on very small scales where shell crossing dominates, we stop decreasing the smoothing scale at $R_{\mathrm{min}}$. A reasonable choice may be $R_{\mathrm{min}}\sim 1\ h^{-1}\text{Mpc}$. We work with $R_{\mathrm{min}}=1.01\,L/N_{\mathrm{grid}}$ throughout, which gives $R_{\mathrm{min}}=0.99\ h^{-1}\text{Mpc}$ for our small-volume simulation and $R_{\mathrm{min}}=2.7\ h^{-1}\text{Mpc}$ for our large-volume simulations.

For the other parameters, we can use simple heuristics. The number of iteration steps for the displacement $\bm{\chi}$ can be determined by monitoring the final quantity of interest and stopping the iteration once that quantity stops changing significantly. Fig. 12 demonstrates that less than ten iteration steps should be sufficient to achieve convergence for most applications. The size of the regular grid should be chosen such that the smallest length scale of interest is still resolved by the grid. We work with $N_{\mathrm{grid}}^{3}=512^{3}$ grid points throughout the postprocessing of the simulations, corresponding to a grid resolution of $\Delta x=0.98\ h^{-1}\text{Mpc}$ for the $L=500\ h^{-1}\text{Mpc}$ simulations, and $\Delta x=2.7\ h^{-1}\text{Mpc}$ for the $L=1380\ h^{-1}\text{Mpc}$ simulations that we use for studying the BAO scale. To reduce potential aliasing effects, we truncate small-scale modes with $k>k_{\mathrm{max}}=2\pi/L\times N_{\mathrm{grid}}/2$.

The optional displacement factor $\epsilon_{s}$ can be set to less than unity to avoid potential overshooting when displacing objects, for example if the overdensity is large and the Zeldovich displacement may not be appropriate. However, if the initial smoothing scale is chosen sufficiently large, for example, $R_{\mathrm{init}}\gtrsim 5\ h^{-1}\text{Mpc}$ at $z=0.6$, the overdensity tends to be less than unity and full displacements with $\epsilon_{s}=1$ seem to work well.

E.2 Convergence of simulations

As a basic check for convergence of the FastPM simulations, we ran a simulation with $40$ time steps linearly spaced between $a=0.1$ and $a=1$, and a second more accurate simulation with $120$ time steps linearly spaced between $a=0.01$ and $a=1$. Both simulations used $2048^{3}$ particles, box size $L=500\ h^{-1}\text{Mpc}$, and we apply reconstruction to a $1\%$ dark matter subsample at $z=0.6$. The correlation coefficient of the reconstructed density with the initial conditions differs by less than $1\%$ between these two simulations, for all reconstruction methods considered in this paper, and for any number of iterations steps used in the reconstruction procedure. This indicates that the simulations have converged in the sense that the final result is robust against changes of starting time and number of time steps used to run the simulations.

Appendix F Modeling the new reconstruction method

Given a prescription for reconstruction, we can try to model the statistics of the reconstructed density, similar to previous efforts Padmanabhan et al. (2009); Sherwin and Zaldarriaga (2012); Schmittfull et al. (2015); White (2015); Baldauf et al. (2015); Cohn et al. (2016); Hikage et al. (2017) modeling the density after the standard reconstruction of Eisenstein et al. (2007a).

F.1 Undoing shift terms

From the mapping between Lagrangian and Eulerian space $\bm{x}=\bm{q}+\bm{\psi}$ and mass conservation

[TABLE]

we have that $1+\delta=1/J_{q}$, where $J_{q}$ is the determinant of the Jacobian matrix

[TABLE]

where $\delta^{K}_{ij}$ is the Kronecker delta (not to be confused with density contrast). Spatial indices are always raised and lowered using $\delta^{K}_{ij}$ and its inverse. When convenient we will use a comma to denote spatial derivatives. We can obtain the final density directly from $J_{q}$ which we have solved explicitly in terms of ${\bm{q}}$. Using (23)

[TABLE]

This expression should be used perturbatively expanding up to a given order in the displacement field.

$J_{q}$ depends only on displacement gradients $\psi_{i,j}$ while there are nonlinear shift terms from the mapping from ${\bm{q}}$ to $\bm{x}={\bm{q}}+\bm{\psi}$ that depend on the displacement $\psi_{i}$ itself. These two sets of terms can have different sizes. The shift terms can be large and should be resummed at the BAO scale if one wants to get a good estimate of the correlation function Senatore and Zaldarriaga (2015); Baldauf et al. (2015); even a linear long-wavelength component of the displacement can move particles over a large enough distance that truncating a perturbative expansion of the form $f({\bm{q}}+\bm{\psi})=f({\bm{q}})+\bm{\psi}\cdot{\bm{\nabla}}f+\cdots$ would lead to large errors. The goal of reconstruction is to undo the shift terms so that the nonlinearities arising from these shift terms are minimized.

F.2 Density of particles shifted by $\bm{\chi}$

Given a displacement field $\bm{\chi}$, let us make a coordinate transformation from the final Eulerian coordinates $\bm{x}$ to new coordinates $\hat{\bm{q}}$ defined by

[TABLE]

If the displacement $\bm{\chi}$ were to coincide with the true Lagrangian-to-Eulerian displacement $\bm{\psi}$, then the $\hat{\bm{q}}$ coordinates would equal the true Lagrangian ${\bm{q}}$ coordinates. In that case the density of particles in $\hat{\bm{q}}$ space would be uniform. If $\bm{\chi}$ is an approximate estimate of the true displacement $\bm{\psi}$, then $\hat{\bm{q}}$ are estimated approximate Lagrangian coordinates, and the density of particles in $\hat{\bm{q}}$ space is approximately uniform (the level to which they are uniform is shown in the upper panel of Fig. 1 that is based on progressively better $\bm{\chi}$).

Similar to above, we can use mass conservation to relate the density of particles in $\hat{\bm{q}}$ and $\bm{x}$ spaces,

[TABLE]

so that

[TABLE]

where $J_{\hat{q}}$ is the determinant of the Jacobian matrix

[TABLE]

We can now compute the Fourier transform of the density field in Fourier space $\Delta({\bm{k}})$:

[TABLE]

where in the last two expressions $\bm{\chi}$ should be thought of as a function of $\bm{x}$ and ${\bm{q}}$, respectively, through the relation

[TABLE]

Ideally the density in Lagrangian space would be uniform and thus $\Delta({\bm{k}})$ would be zero for nonzero ${\bm{k}}$. The density is slightly nonuniform if the estimated and true Lagrangian coordinates differ, $\hat{\bm{q}}\neq{\bm{q}}$, which is the case if the estimated and true displacements differ, $\bm{\chi}\neq\bm{\psi}$.

Given $\bm{\chi}$ one can construct the density in $\hat{\bm{q}}$ space and measure $\Delta({\bm{k}})$ directly from the data. The goal of this appendix is to compute the statistics of $\Delta({\bm{k}})$. One should keep in mind that $\bm{\chi}$ will also be computed from the data itself, it will be a function of $\delta$. In the EFT one has a perturbative expansion for $\bm{\psi}$ and $\delta$.

F.3 Displacement field

We can obtain an expression for $\bm{\chi}$ if we demand that the density field computed using the Zeldovich approximation starting with the displacement $\bm{\chi}$ equaled the filtered version of the density. The Zeldovich approximation has all the nonlinearities from the displacements but no dynamical interaction among the different modes. We are basically trying to solve for the nonlinear displacement, of course once one is inside the nonlinear regime, this is not really possible. That is to say we will find a displacement field that produces the final density from a uniform distribution, but this is not unique as we can always exchange particles after we are done in any given solution and create a new solution. Thus what we get on nonlinear scales is a bit random.

If we equated the Zeldovich density to the filtered one, we would obtain:

[TABLE]

where ${\bm{k}}={\bm{p}}_{1}+\cdots+{\bm{p}}_{n}$ and $F_{Z}$ are the kernels of the Zeldovich approximation:

[TABLE]

We then solve the equation perturbatively to obtain:

[TABLE]

In general we have:

[TABLE]

The perturbative version of the solution has

[TABLE]

where $[...]_{\mathrm{sym}}$ stands for symmetrized in the momenta. In standard reconstruction one is just using the linear version of this equation:

[TABLE]

F.4 Flavors of reconstruction

In standard reconstruction a uniform field in $\bm{x}$ space is generated which is then shifted with the same displacement $\bm{\chi}$. Let us call this field $\Delta_{S}$,

[TABLE]

where we have used the same expressions used to compute $\Delta({\bm{k}})$ but in this case the overdensity in $\bm{x}$ space is zero. In the standard reconstruction algorithm, one estimates $\delta_{0}$ as:

[TABLE]

If one has iteratively or perturbatively found a displacement field that recovers the density one can perhaps use directly that field as the estimate of the linear density field. This is what we do in this paper,

[TABLE]

We could also divide by $W(k)$ to undo the filtering so as to recover the linear density field at lowest order.

We compare with standard reconstruction in the main text as well as in Appendix B. Of course our new estimate of the initial density only makes sense if one has solved for $\delta_{\chi}({\bm{k}})$ iteratively; otherwise at linear order one has done nothing.

F.5 Perturbative results

In this section we will write down the expression for the reconstructed field in perturbation theory. We will first assume that the nonlinear density can be expressed as:

[TABLE]

We will need the first few kernels

[TABLE]

We will use these expressions to obtain a formula for the reconstructed fields in terms of the initial condition $\delta_{0}$,

[TABLE]

In standard reconstruction

[TABLE]

We can define:

[TABLE]

in terms of which we find:

[TABLE]

The solution for $\tilde{\delta}$ reads:

[TABLE]

For $J$ (which we only need to second order) we have:

[TABLE]

Finally for $\delta_{\chi}$ we have:

[TABLE]

The easiest way to get these formulas is to simply equate

[TABLE]

and solve for $\delta_{\chi}$.

One can easily use the above formulas to infer some properties of the reconstructed field. For example we can look at the quadratic kernels:

[TABLE]

It is easy to show that in the regime when $W\approx 1$ all the shift terms are canceled from both formulas; in fact,

[TABLE]

Indeed, this is true for the higher order kernels as well; the shift term are cancelled when $W\rightarrow 1$. Furthermore in the limit in which $p_{1},p_{2}\gg k$ both kernels have a nice UV limit; they scale as $(k/p)^{2}$. In fact, in this limit

[TABLE]

with $\mu={\bm{p}}\cdot{\bm{k}}/pk$. Thus the kernels have the correct UV limit regardless of the window. This is again still true for the higher order kernels. In fact, both reconstruction formulas look very much the same, they only differ in terms proportional to the gravitational interaction. From the fact that the second stage of the new and standard reconstructions differ only slightly in perturbation theory, we expect that for a given displacement, the second stage of the new and standard reconstruction should perform similarly. Therefore it makes sense that the extended standard reconstruction in Appendix B, which is using the eight-step displacement field $\bm{\chi}^{(8)}$, performs roughly as well as the new $\mathcal{O}(1)$ reconstruction based on the same displacement field (see Fig. 12). The key improvement from the new reconstruction comes from having a better displacement field (the second order correction and the reduced coefficients of the growth and tidal term play a much smaller role).

One could improve the reconstruction by removing the quadratic and cubic pieces. At one loop the cubic piece could be incorporated in a transfer function on the linear field. Thus one could estimate the initial density using

[TABLE]

and choose $t_{1}$ and $t_{2}$ to minimize the difference with $\delta_{0}$ in simulations as described in Appendix A. We tested this approach in the main text.