The Reconstructed Power Spectrum in the Zeldovich Approximation

TL;DR

This paper develops an analytical model for density reconstruction using the Zeldovich approximation, accurately fitting real and redshift-space data and improving understanding of BAO signal sharpening.

Contribution

It introduces a comprehensive Zeldovich-based model including counterterms and bias terms up to quadratic order, fitting data across various scales and space configurations.

Findings

Sub-percent agreement with N-body data in real and redshift space

Effective modeling of BAO signal sharpening after reconstruction

Comparison with existing models highlights improvements and higher bias effects.

Abstract

Density-reconstruction sharpens the baryon acoustic oscillations signal by undoing some of the smoothing incurred by nonlinear structure formation. In this paper we present an analytical model for reconstruction based on the Zeldovich approximation, which for the first time includes a complete set of counterterms and bias terms up to quadratic order and can fit real and redshift-space data pre- and post-reconstruction data in both Fourier and configuration space over a wide range of scales. We compare our model to n-body data at $z = 0$ from the {\tt DarkSky} simulation \cite{skillman14}, finding sub-percent agreement in both real space and in the redshift-space power spectrum monopole out to $k = 0.4\ h$ Mpc$^{-1}$, and out to $k = 0.2\ h$ Mpc$^{-1}$ in the quadrupole, with comparable agreement in configuration space. We compare our model with several popular existing alternatives, updating existing theoretical results for exponential damping in wiggle/no-wiggle splits of the BAO signal and discuss the usually-ignored effect of higher bias contributions on the reconstructed signal. In the appendices, we re-derive the former within our formalism, present exploratory results on higher-order corrections due to nonlinearities inherent to reconstruction, and present numerical techniques with which to calculate the redshift-space power spectrum of biased tracers within the Zeldovich approximation.

Tables

Table 1. Table 1 : Number densities and bias values for the halo samples we use. Halo masses are log 10 subscript 10 \log_{10} of the mass in h − 1 ​ M ⊙ superscript ℎ 1 subscript 𝑀 direct-product h^{-1}M_{\odot} , number densities are times 10 − 3 ​ h 3 ​ Mpc − 3 superscript 10 3 superscript ℎ 3 superscript Mpc 3 10^{-3}\,h^{3}{\rm Mpc}^{-3} .

$\lg M$	Redshift	$\overline{n}$	$b$
$12.0-12.5$	0.0	$3.45$	0.87
$12.5-13.0$	0.0	$1.18$	1.05
$13.0-13.5$	0.0	$0.38$	1.30

Table 2. Table 2 : Table of power spectrum contributions in M ii .

${𝒜}_n$	$n=0$	$n=1$	$n=2$	…
$1$	$1$	0	0
$b_1$	0	$2KU\left(q\right)$	0
$b_1^2$	${\xi }_L\left(q\right)$	0	$K^{2}U{\left(q\right)}^2$
$b_2$	0	0	$K^{2}U{\left(q\right)}^2$
$b_1{b}_2$	0	$KU\left(q\right){\xi }_L\left(q\right)$	0
$b_2^2$	$\frac{1}{2}{\xi }_L{\left(q\right)}^2$	0	0
$b_s$	$-K^2{X}_{s^2}\left(q\right)$	0	$K^2{Y}_{s^2}\left(q\right)$
$b_1{b}_s$	0	$2KV\left(q\right)$	0
$b_2{b}_s$	${\chi }^{12}\left(q\right)$	0	0
$b_s^2$	$\zeta \left(q\right)$	0	0

Equations

$$1+{\delta }^a(\text{x},\tau )=\int d^3\text{q}{F}^a(\text{q}){\delta }_D\left(\text{x}-\text{q}-\mathbf{\Psi }(\text{q},\tau )\right).$$

$$P^{ab}(k)=\int d^3\text{q}{e}^{i\text{k}\cdot \text{q}}⟨F^a({\text{q}}_2)F^b({\text{q}}_1)e^{i\text{k}\cdot {\Delta }^{ab}}{⟩}_{q=|q_2-q_1|},{\Delta }^{ab}={\mathbf{\Psi }}^a({\text{q}}_2)-{\mathbf{\Psi }}^b({\text{q}}_1),$$

$$F^a(\text{q})=1+b_1^a{\delta }_0(\text{q})+\frac{1}{2}{b}_2^a\left({\delta }_0(\text{q}{)}^2-⟨{\delta }_0^2⟩\right)+b_s^a\left(s^2(\text{q})-⟨s^2⟩\right)+b_{{\nabla }^2}^a{\nabla }_q^2{\delta }_0(\text{q})+\cdots ,$$

$$P^{ab}(k)=$$

$$-\frac{1}{2}{k}_i{k}_j(b_2^b{U}_i^a{U}_j^a+b_2^a{U}_i^b{U}_j^b+2b_1^a{b}_1^b{U}_i^a{U}_j^b)+i{k}_i(b_2^b{b}_1^a{U}_i^a+b_1^b{b}_2^a{U}_i^b){\xi }_L$$

$$-\frac{1}{2}{k}_i{k}_j(b_s^a{{\rm Y}}_{ij}^b+b_s^a{{\rm Y}}_{ij}^b)+i{k}_i(b_1^a{b}_s^b{V}_i^{ab}+b_1^b{b}_s^a{V}_i^{ba})$$

$$\displaystyle+\frac{1}{2}(b_{2}^{a}b_{s}^{b}+b_{2}^{b}b_{s}^{a})\chi^{12}+b_{s}^{a}b_{s}^{b}\zeta+\cdots\Big{]}$$

$$A_{ij}^{ab}=⟨{\Delta }_i^{ab}{\Delta }_j^{ab}⟩,U_i^b=⟨{\Delta }_i^{ab}{\delta }_0({\text{q}}_2)⟩,{\xi }_L=⟨{\delta }_0({\text{q}}_2){\delta }_0({\text{q}}_1)⟩$$

$$\zeta =⟨s^2({\text{q}}_2)s^2({\text{q}}_1)⟩,{{\rm Y}}_{ij}^b=⟨{\Delta }_i^{ab}{\Delta }_j^{ab}{s}^2({\text{q}}_2)⟩,V_i^{ab}=⟨{\Delta }_i^{ab}{\delta }_0({\text{q}}_2)s^2({\text{q}}_1)⟩,{\chi }^{12}=⟨{\delta }_0^2({\text{q}}_1)s^2({\text{q}}_2)⟩.$$

$$k_i{k}_j⟨({\mathbf{\Psi }}_i+{\hat{n}}_i{\hat{n}}_l{\dot{\mathbf{\Psi }}}_l)({\mathbf{\Psi }}_j+{\hat{n}}_j{\hat{n}}_m{\dot{\mathbf{\Psi }}}_m)⟩\equiv k^2(X(0)+(2\dot{X}(0)+\ddot{X}(0)){\nu }^2),$$

$$\mathbf{\chi}_{\mathbf{k}}=-\frac{i\mathbf{k}}{k^{2}}\mathcal{S}(k)\ \Big{(}\frac{\delta_{g}(\mathbf{k})}{b+f\nu^{2}}\Big{)}\approx-\mathcal{S}(k)\boldsymbol{\Psi}^{(1)}(\textbf{k})$$

$$1+{\delta }_d(\mathbf{r})$$

$$=\int d^3\text{x}\int d^3\text{q}F(\text{q}){\delta }_D\left[\text{x}-\text{q}-\mathbf{\Psi }(\text{q})\right]{\delta }_D\left[\mathbf{r}-\text{x}-{\chi }_d(\text{x})\right]$$

$$=\int d^3\text{q}F(\text{q}){\delta }_D\left[\text{r}-\text{q}-\mathbf{\Psi }(\text{q})-{\chi }_d(\text{q}+\mathbf{\Psi }(\text{q}))\right],$$

$$(2\pi {)}^3{\delta }_D(\text{k})+{\delta }_d(\text{k})$$

$$(2\pi {)}^3{\delta }_D(\text{k})+{\delta }_s(\text{k})$$

$${\mathbf{\Psi }}^d=\mathbf{\Psi }+{\chi }_d,{\mathbf{\Psi }}^s={\chi }_s,$$

$${\mathbf{\Psi }}^d(\text{k})=\left[1-\mathcal{S}(k)\right]\mathbf{\Psi }(\text{k}),{\mathbf{\Psi }}^s(\text{k})=-\mathcal{S}(k)\mathbf{\Psi }(\text{k}).$$

$$P_L^{dd}(k)={\left[1-\mathcal{S}(k)\right]}^2{P}_L(k),P_L^{ds}(k)=-\mathcal{S}(k)\left[1-\mathcal{S}(k)\right]P_L(k),P_L^{ss}(k)=\mathcal{S}(k{)}^2{P}_L(k)$$

$$P_L^{dm}(k)=\left[1-\mathcal{S}(k)\right]P_L(k),P_L^{sm}(k)=-\mathcal{S}(k)P_L(k).$$

$$A_{ij}^{ds}(\text{q})\stackrel{q\to 0}{=}⟨{\mathbf{\Psi }}_i^d{\mathbf{\Psi }}_j^d⟩+⟨{\mathbf{\Psi }}_i^s{\mathbf{\Psi }}_j^s⟩-2⟨{\mathbf{\Psi }}_i^d{\mathbf{\Psi }}_j^s⟩\equiv {\Sigma }^2{\delta }_{ij},$$

$$A^{ds}_{ij}(\textbf{q})\stackrel{{\scriptstyle q\rightarrow\infty}}{{=}}\langle\boldsymbol{\Psi}^{d}_{i}\boldsymbol{\Psi}^{d}_{j}\rangle+\langle\boldsymbol{\Psi}^{s}_{i}\boldsymbol{\Psi}^{s}_{j}\rangle\equiv\Big{(}\Sigma_{dd}^{2}+\Sigma_{ss}^{2}\Big{)}\delta_{ij},$$

$$X^{dd}(q)$$

$$X^{ds}(q)$$

$$X^{ss}(q)$$

$$P_{\mathrm{sym}}(\text{k})=(b+f{\nu }^2{)}^2{P}_{\mathrm{L}}(k)+\mathcal{O}(P_{\mathrm{L}}^2)$$

$$\displaystyle P_{\rm iso}(\textbf{k})=\Big{[}(b+f\nu^{2})(1-\mathcal{S})+b\ \mathcal{S}\Big{]}^{2}P_{\rm L}(k)+\mathcal{O}(P_{\rm L}^{2}),$$

$${\mathbf{\Psi }}^d(\text{k})=\left[1-\mathcal{S}(k)\right]\mathbf{R}\mathbf{\Psi }(\text{k}),{\mathbf{\Psi }}^s(\text{k})=-\mathcal{S}(k)\mathbf{R}\mathbf{\Psi }(\text{k}).$$

$${\mathbf{\Psi }}^d(\text{k})=\left[1-\mathcal{S}(k)\right]{\mathbf{\Psi }}^R(\text{k})=\left[1-\mathcal{S}(k)\right]\mathbf{R}\mathbf{\Psi }(\text{k}),{\mathbf{\Psi }}^s(\text{k})=-\mathcal{S}(k)\mathbf{\Psi }(\text{k}).$$

$$A_{ij}^{ds,\mathrm{Iso}}$$

$$=R_{in}{R}_{jm}⟨{\mathbf{\Psi }}_n^d{\mathbf{\Psi }}_m^d{⟩}_{\mathrm{RealSpace}}+⟨{\mathbf{\Psi }}_i^s{\mathbf{\Psi }}_j^s{⟩}_{\mathrm{RealSpace}}-2R_{in}⟨{\mathbf{\Psi }}_n^d({\text{q}}_2){\mathbf{\Psi }}_j^s({\text{q}}_1){⟩}_{\mathrm{RealSpace}},$$

$$k_i{k}_j⟨{\mathbf{\Psi }}_i^d({\text{q}}_2){\mathbf{\Psi }}_j^s({\text{q}}_1)⟩=k_i{k}_j({\delta }_{ik}$$

$$=k^2(1+f{\nu }^2){\tilde{X}}^{ds}+k^2({\mu }^2+f\mu \nu (\hat{q}\cdot \hat{n})){\tilde{Y}}^{ds},$$

$$⟨{\mathbf{\Psi }}^d({\text{q}}_2){\mathbf{\Psi }}^s({\text{q}}_1){⟩}_{\mathrm{Real}\mathrm{Space}}={\tilde{X}}^{ds}(q){\delta }_{ij}+{\tilde{Y}}^{ds}(q){\hat{q}}_i{\hat{q}}_j.$$

$$\int \frac{d\varphi }{2\pi }{e}^{A\mu \sqrt{1-{\mu }^2}\cos \varphi }=\sum _{\ell =0}^{\infty }{H}_{\ell }^{(0)}(A)(A{\mu }^2{)}^{\ell },$$

$$H_{\ell }^{(0)}(A)=\sum _{m=0}^{\ell }\frac{(-1{)}^{\ell -m}{A}^{2m-\ell }\Gamma (m+\frac{1}{2})}{\sqrt{\pi }\Gamma (2m+1)\Gamma (2m-\ell +1)\Gamma (\ell -m+1)}.$$

Full text

1 Introduction

Density field reconstruction [2] is a means of improving the determination of the distance-redshift relation using baryon acoustic oscillations (BAO) [3]. The BAO method is a “standard ruler” test which seeks to measure the scale of a feature in the 2-point function whose physical size is known. Comparison with the observed size of this feature gives the angular diameter distance and Hubble parameter as a function of redshift. While the large size of the BAO feature ($100\,$Mpc) makes it relatively immune to systematic effects, nonlinear evolution erases the the oscillations on small scales, or broadens the peak in the correlation function, and reduces the accuracy with which the scale can be measured [4, 5, 6, 7, 8]. However much of the peak broadening comes from motions sourced by very long wavelength fluctuations [8] which are well measured by surveys aiming to measure BAO. This insight led ref. [2] to propose that density-field reconstruction could be applied to regain much of the information lost to non-linearities. It has been used in all recent BAO surveys to improve their constraints (e.g. see ref. [9] and references therein).

BAO reconstruction has been studied both numerically [10, 11, 12] and analytically [2, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27]. Our work builds upon these analytic calculations. Where earlier work made simplifications aimed at highlighting important physical effects, neglected complications such as redshift-space distortions, applied heuristics or otherwise simplified the calculations for explanatory effect, we aim to produce a consistent dynamical theory which can be compared directly to upcoming observational data. Hence we generalize these calculations to also consider the power spectrum and extend the model to include the complete set of quadratic bias terms. To our knowledge this is the first dynamical model with a full bias scheme that can produce consistent real and redshift-space results in both Fourier and configuration space, allowing it to be used for consistent fitting of upcoming data.

There has been significant theoretical work on reconstruction since the first algorithm [2] was suggested. Most recently, a variety of iterative or alternative reconstruction approaches have been developed [15, 22, 20, 23, 24, 25]. Though our calculations give some insights into these methods, for near-future experiments and for BAO scales these iterative methods do not lead to significant improvements and so we defer consideration of these more complex algorithms to future work.

The outline of this paper is as follows. Section 2 reviews the formalism of Lagrangian perturbation theory within which we work. Section 3 describes the reconstruction algorithm we seek to model, while Section 4 gives our results in Fourier space, comparing to the configuration-space results where appropriate (Section 5). We discuss alternative statistics in Section 6. To assess the range of validity of our models we compare to N-body simulations in Section 7. A comparison with earlier work is given in Section 8 before we conclude in Section 9. Some technical details are elaborated in the appendices.

2 Lagrangian Perturbation theory

The Lagrangian framework [28, 29, 30, 31, 32, 33, 34, 35] describes cosmological structure formation by tracking the displacements $\boldsymbol{\Psi}(\textbf{q})$ of infinitesimal parcels of the matter fluid from their initial (Lagrangian) positions $\textbf{q}.$ In this picture the present day matter over- and underdensities are a result of the clustering of the displaced Eulerian positions $\textbf{x}(\textbf{q},\tau)=\textbf{q}+\boldsymbol{\Psi}(\textbf{q},\tau)$. The displacements follow the equation of motion $\boldsymbol{\Psi}^{\prime\prime}(\textbf{q})+\mathcal{H}\boldsymbol{\Psi}^{\prime}(\textbf{q})=-\nabla_{x}\Phi(\textbf{x}),$ where $\Phi(\textbf{x})$ is the gravitational potential which is in turn sourced by the clustered matter fluid via Poisson’s equation $\nabla^{2}\Phi(\textbf{x},\tau)=\frac{3}{2}\Omega_{m}(\tau)\mathcal{H}^{2}(\tau)\delta(\textbf{x},\tau)$ with $\tau$ the conformal time. This set of equations can be solved perturbatively in terms of the linear overdensity, $\delta_{0}$, and the first order solution is given by $\boldsymbol{\Psi}=-D(\tau)\mathbf{\nabla}\nabla^{-2}\delta_{0},$ where $D(\tau)$ is the linear growth factor [28].

The Lagrangian picture treats tracer bias and advection separately. Given a biased tracer, $a$, with initial overdensity $F^{a}(\textbf{q})=F^{a}\left[\partial^{2}\Phi(\textbf{q}),...\right]$, the time-evolved tracer overdensity at conformal time $\tau$ is given by number conservation as [34]

[TABLE]

The cross power spectrum between two biased tracer populations $a$ and $b$ is then

[TABLE]

where we have used that the integrated expectation value can only depend on $\textbf{q}=\textbf{q}_{2}-\textbf{q}_{1}$, due to the translation invariance of the underlying theory. The bias functionals, $F^{a,b}$, can be Taylor expanded in terms of bias coefficients

[TABLE]

where $s^{2}=s_{ij}s_{ij}$ is the square of the shear field, i.e. the traceless part of $\partial\partial\Phi$. Following ref. [36], we also consider contributions from a “derivative bias” $b_{\nabla^{2}}$, i.e. corrections to the bias expansion at scales close to the halo radius $R_{h}$ proportional to $\nabla^{2}\delta_{0}$; such contributions will, however, be essentially degenerate with counterterms renormalizing nonlinearities in the Zeldovich power spectrum and we will therefore not enumerate them separately in the rest of this work unless otherwise stated.

In this work our focus will be on modelling reconstruction within the Zeldovich approximation [28, 37], which keeps only the linear order term in the dynamics of $\boldsymbol{\Psi}$ but re-sums the effects of the displacement to all orders in a Galilean-invariant manner (this is true for reconstruction also if we take it to mean that all displacements transform the same way). This is specifically accomplished by evaluating the exponential in Equation 2.2 via the cumulant expansion, and evaluating the bias expansion using functional derivatives (see e.g. refs. [34, 35, 36]). Following standard techniques, as outlined in the references above, the resulting expression for the cross spectrum is

[TABLE]

where we have defined111These are generalizations of the similar auto-spectrum quantities defined in refs. [35, 37, 36]. the quadratic two point functions

[TABLE]

and shear correlators

[TABLE]

Note that in the above calculations we have, without loss of generality, associated tracers $a$ and $b$ with Lagrangian positions $\textbf{q}_{2}$ and $\textbf{q}_{1}$, respectively. The quantities in Equation 2.4 with $a$ and $b$ swapped can also be calculated by swapping the positions $\textbf{q}_{1}\leftrightarrow\textbf{q}_{2}$. As an example, $U^{b}=-\langle\boldsymbol{\Psi}^{b}(\textbf{q}_{1})\delta_{0}(\textbf{q}_{2})\rangle$ is the two-point function between the displacement of tracer $b$ and the matter overdensity. The vector and tensor two point functions defined above can be decomposed via rotational symmetry into scalar components, e.g. $A_{ij}=X(q)\delta_{ij}+Y(q)\hat{q}_{i}\hat{q}_{j}$ and $U_{i}=U(q)\hat{q}_{i}$. Formulae for these functions, expressed as Hankel transforms of power spectra, are given in Appendix A. Finally, we include the contribution $\alpha_{0}k^{2}$ in the square brackets of Equation 2.4 as the lowest-order counterterm renormalizing sensitivities to small-scale power in $A_{ij}$ — in practice this simply modifies the matter contribution $P_{\rm Zel}(k)$ ($\propto 1$ in the square brackets) to $(1+\alpha_{0}k^{2})P_{\rm Zel}(k)$ (see e.g. refs. [38, 39]). Each term in Equation 2.4 can be evaluated as Hankel transforms (see e.g. ref. [36]) using the identities given at the end of [38], which we carry out using the mcfit package222https://github.com/eelregit/mcfit.

The Lagrangian formalism allows a straightforward translation between real and redshift space via a mapping of the Lagrangian displacements. In particular, assuming the plane-parallel approximation333This should be an excellent approximation on BAO scales [40], but if necessary the formalism can be modified to include “wide-angle” effects [41]. and working in the Zeldovich approximation, quantities in redshift space are given simply by substituting $\boldsymbol{\Psi}_{i}\rightarrow\boldsymbol{\Psi}^{R}_{i}=R_{ij}\boldsymbol{\Psi}_{j}$ [34]. Here $R_{ij}=\delta_{ij}+f\hat{n}_{i}\hat{n}_{j}$, where $\hat{n}$ denotes the line-of-sight direction and $f=d\ln D/d\ln a$ is the linear-theory growth rate. To lowest order, transforming into redshift space requires the inclusion of a second counterterm dependent on the line-of-sight angle $\nu=\hat{k}\cdot\hat{n}.$ We can see this explicitly, for example, in the UV-sensitive zero-lag term in $A_{ij}$, which gains an angular dependence

[TABLE]

where $\dot{\boldsymbol{\Psi}}$ is the velocity in Hubble units equal to $f\boldsymbol{\Psi}$ in the Zeldovich approximation444See refs. [42, 36] for a more detailed exposition of the “dot notation.”; roughly speaking, we need one angle-independent counterterm $\alpha_{0}k^{2}$ to absorb the UV dependence of $X(0)$ and another $\alpha_{2}k^{2}\nu^{2}$ to absorb the UV dependence of the velocities. The complete set of counterterms in redshift space thus makes a contribution of the form $(\alpha_{0}+\alpha_{2}\nu^{2})k^{2}P_{ZA}(\textbf{k});$ since $P_{ZA}(\textbf{k})$ is equal to $(1+f\nu^{2})^{2}P_{L}(k)$ to linear order, an equivalent viewpoint—which we will adopt in this work— is to have constant counterterms $\bar{\alpha}_{0}k^{2}$ and $\bar{\alpha}_{2}k^{2}$ for the monopole and quadrupole, respectively, where the barred counterterms are linear combinations of the unbarred quantities.

3 Reconstruction algorithm

In this section we describe two possible methods for reconstruction in redshift space, both built around the Zeldovich approximation. The standard procedure for reconstruction was developed in ref. [2] and involves displacing both observed galaxies and a spatially uniform distribution by a calculated shift field, $\chi$, then taking the relative density contrast between the two sets of particles as the reconstructed density field. For a suitably chosen $\chi$, this can reduce the effect of large scale (IR) bulk flows that “blur” the BAO feature. However there is no consensus in the community on the correct procedure for handling redshift-space distortions: the implementation in ref. [43] chose to multiply $\chi$ by $1+f$ in the line-of-sight direction for $\delta_{d}$ but not for $\delta_{s}$. This ‘undoes’ the supercluster infall effect [44] and reduces the $\ell>0$ moments of the 2-point function on large scales. Ref. [16] suggested a symmetric treatment of $\delta_{d}$ and $\delta_{s}$, which recovers linear theory on large scales. This is more natural from the point of view of perturbation theory and better behaved near the boundaries, but is less often implemented on data. A number of other choices were explored in ref. [18] but in this work we will restrict our attention to the two methods described above.

The reconstruction procedure consists of the following steps [2]:

Smooth the observed galaxy density field $\delta_{g}$ with a kernel $\mathcal{S}$ to filter out small scale (high $k$) modes, which are difficult to model. We use a Gaussian smoothing of scale $R_{s}$, specifically $\mathcal{S}(k)=\exp[-(kR_{s})^{2}/2]$, though none of our analytic results will depend specifically on this choice. For galaxy surveys Gaussian smoothing has been universally adopted (though with different conventions for $R_{s}$) but in other contexts it may be advantageous to implement a Wiener filter instead (e.g. ref. [19]). 2. 2.

Compute the shift, $\mathbf{\chi}$, by dividing the smoothed galaxy density field by a bias factor $b$ and linear RSD factor [44] and then take the inverse gradient. Assuming linear theory with scale-independent bias and supercluster infall holds on large scales, the calculated shift field should approximate the negative smoothed Zeldovich displacement. In a simulation with a periodic box, these first two steps can be implemented using FFTs as

[TABLE]

where the bias factor is related to the Lagrangian first-order bias by $b=1+b_{1}$ and we have defined the line-of-sight angle $\nu=\hat{n}\cdot\hat{k}$. For non-periodic data the relevant differential equation can be solved by multigrid555https://github.com/martinjameswhite/recon_code or by linear algebra techniques [43] or iteratively using FFTs [45]. 3. 3.

Move the galaxies by $\chi_{d}=\mathbf{R}\chi$ and compute the “displaced” density field, $\delta_{d}$. 4. 4.

Shift an initially spatially uniform distribution of particles by

• •

Rec-Sym: $\chi_{s}=\mathbf{R}\chi,$ i.e. the same amount as the observed galaxies, or,

• •

Rec-Iso: The un-redshifted $\chi_{s}=\chi$.

to form the “shifted” density field, $\delta_{s}$. Note that we have borrowed the nomenclature of ref. [18] for the latter, which “isotropizes” the reconstructed field on large scales. For the former we use “Rec-Sym” to indicate the symmetry of the treatment of $\delta_{d}$ and $\delta_{s}$. 5. 5.

The reconstructed density field is defined as $\delta_{r}\equiv\delta_{d}-\delta_{s}$ with power spectrum $P_{r}(k)\propto\langle\left|\delta_{r}^{2}\right|\rangle$.

Throughout we shall assume that the fiducial cosmology and halo bias are properly known during reconstruction (see e.g. refs. [27, 46] for relaxation of this assumption), and take the approximation in Eq. 3.1 to be exact. For further discussion of this point see refs. [17, 21]. The procedure in real space can be straightforwardly obtained by setting $f=0$, in which case Rec-Sym and Rec-Iso become equivalent. Taking the limit $\mathcal{S}\to 0$ or $\chi\to 0$ returns the ‘raw’ spectrum, before reconstruction.

4 Reconstructed power spectrum

There has been significant earlier work on modeling density-field reconstruction within perturbation theory [2, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27]. In particular ref. [16] presented a calculation of the configuration-space two-point function (the correlation function) under the assumption of Zeldovich dynamics and that $\chi=-\mathcal{S}\boldsymbol{\Psi}$. In this paper we generalize that calculation to a more complete bias model (see §5), including all terms allowed by symmetries up to quadratic order as well as a proper set of counterterms, and we show how to implement the model in Fourier space. We have explicitly checked that the Hankel transform of our Fourier-space expressions matches the direct configuration-space calculation to 1% in all terms, and we release code which makes consistent predictions for both statistics with a common set of parameters. To our knowledge this is the first calculation which provides self-consistent predictions in both spaces, uses a dynamical rather than a heuristic model, works in redshift space and has a full set of bias and counterterms.

Our focus in this section will be to model the reconstructed power spectrum using Lagrangian perturbation theory in both real and redshift space (the expression for the ‘propagator’ is given in Appendix B for completeness). Following the algorithm outlined above, the reconstructed power spectrum in real space is given by $P_{\rm recon}=P^{dd}+P^{ss}-2P^{ds}.$ Within the Lagrangian framework we can write the displaced density field as

[TABLE]

where we performed the x integral using the first $\delta$-function to go from the second to third lines. Importantly while the fluid displacement, $\boldsymbol{\Psi}$, is evaluated at the Lagrangian position, q, the shift field is evaluated at the shifted Eulerian position, $\textbf{q}+\boldsymbol{\Psi}$. The above equalities hold both when the pre-reconstruction coordinate, x, is in real or redshift space, with the implicit substitution of $\boldsymbol{\Psi}\rightarrow\mathbf{R}\boldsymbol{\Psi}$ in the latter case, as long as the appropriate shift field $\chi_{d}$ is chosen. The expression for the shifted density can be similarly derived or found by setting $\boldsymbol{\Psi}(\textbf{q})=0$ and $\chi_{d}\rightarrow\chi_{s}$ in the above expression. In Fourier space this translates to

[TABLE]

Below we will make the approximation $\chi(\textbf{q}+\boldsymbol{\Psi})\approx\chi(\textbf{q}).$ The nonlinearities from the Lagrangian-to-Eulerian mapping can be understood as a perturbation series in $\boldsymbol{\Psi}/R$, where $R$ is the smoothing scale, and we explore their consequences in Appendix E (see also refs. [17, 21] and the discussion in ref. [16]). Within this approximation we can treat the displaced and shifted field as tracers with displacements

[TABLE]

where the Zeldovich displacements should be understood as being in redshift space for the displaced field and in either redshift or real space for the shifted field depending on the method used. In this picture the “displaced” tracer has the same bias functional as the original galaxies ($F^{d}\equiv F^{g}$) while the “shifted” tracer is unbiased ($F^{s}\equiv 1$). A straightforward consequence of the reconstruction procedure is that, like that of any discrete tracer, the shift field autospectrum will contain an independent shot noise term $P^{ss}_{\rm SN}=1/n_{s},$ where $n_{s}$ is the number density of the uniform random particles. The full shot noise contribution to the reconstructed spectrum is the the sum of the galaxy and random particle shot noises.

4.1 Real space

In real space both the displaced and shifted fields are moved by the same, smoothed negative Zeldovich displacement, $\chi_{d}=\chi_{s}=-S\star\boldsymbol{\Psi}$, such that in Fourier space

[TABLE]

and the auto- and cross-spectra can be calculated using Equation 2.4 and the correlators in Appendix A, using linear theory spectra

[TABLE]

as well as tracer-matter power spectra

[TABLE]

Note that the shifted field is negatively correlated with both the matter and displaced fields by-construction, since the random particles are displaced in the opposite direction of the (smoothed) Zeldovich displacement.

The Lagrangian space two-point correlation functions required to calculate the pre- and post-reconstruction power spectra, normalized to their present-day values, are shown in Figure 1. For simplicity we have excluded the shear correlators and refer readers to Appendix A for further details. The components $X$ and $Y$ describe correlation functions of two displacements, while the $U$’s involve those with only one displacement, such that the former are Hankel transforms of the linear tracer-tracer spectra, while the latter involve the linear tracer-matter spectra. As expected, the $Y$’s and $U$’s for the displaced and shifted fields contain the behavior of the full matter contribution and small and large scales, respectively, and cross correlations between the shifted field and the displaced or matter fields is negative.

The $X(q)$ components however, especially the cross-correlation $X^{ds}$, display more subtle behavior. In particular, we have

[TABLE]

such that $X^{ds}(q)\rightarrow\Sigma^{2}$ as $q\rightarrow 0$. This is because, when evaluated at the same point, $\boldsymbol{\Psi}^{d}-\boldsymbol{\Psi}^{s}=\boldsymbol{\Psi},$ i.e. the difference between the displaced the shifted displacements is none other than the original Zeldovich displacement. This in turn implies that the cross spectrum is damped at small scales $\propto\exp[-k^{2}\Sigma^{2}/2]$ due to the nonzero displacement between the displaced and shifted fields. Similar behavior is seen in the evaluation of unequal-time correlation functions [47] and the baryon-cold dark matter cross-correlation [48, 49], though the physical mechanisms are of course different. At large scales, we similarly have

[TABLE]

such that $X^{ds}$ asymptotes to the average of $X^{dd}$ and $X^{ss}$ at large separations. For completeness, we give explicit expressions for the displaced and shifted $X^{ab}$ here:

[TABLE]

The corresponding expressions for $Y^{ab}$ can be directly obtained by calculating $-3$ times the $j_{2}$ components. As we shall discuss further in Section 8, the signs for the Bessel function coefficients in our expression for $X^{ds}$ differ from those in ref. [26]. We note also, as has been emphasized before [13], each of the three contributions to $P_{\rm recon}$ has a different damping factor which can only be roughly approximated by a single Gaussian term.

The lowest-order bias terms in the reconstructed real-space power spectrum at $z=0$ are shown Figure 2. The pure-matter piece (i.e. the “1” in Equation 2.4) is the only term that includes contributions from all three combinations of $d$ and $s$, while the $b_{1}^{2}$ piece consists of only the $dd$ contribution. While each piece individually differs from the linear power spectrum, compared to the pre-reconstruction power spectrum, the Zeldovich approximation predicts that the post-reconstruction power spectrum largely recovers the oscillatory features in the linear spectrum, as seen in the lower panels of Figure 2. We note that the structure of the breakdown into $P^{dd}$, $P^{ds}$ and $P^{ss}$ shown in Figure 2 proceeds similarly in the higher-order bias contributions: bias terms like $b_{1}^{2}$, that are products of two bias parameters (e.g. $b_{1}b_{2},b_{2}b_{s}$, …), do not involve any displacements ($\boldsymbol{\Psi}$) and can thus only enter in the autospectrum of the biased “$d$” tracer $P^{dd}$, while those like $b_{1}$ that involve only one bias parameter (e.g. $b_{2}$, $b_{s}$) involve two-point functions with one displacement contracted and thus contribute to the cross spectrum $P^{ds}$ as only one of the constituent tracers needs to be biased. The autospectrum $P^{ss}$ does not contain any bias terms.

Figure 3 shows all the contributions to the reconstructed galaxy power spectrum and correlation function, up to the quadratic bias and shear terms. As seen in the top panel, the reduced damping in the lowest-order bias term “wiggles,” barely visible in log-log plots of the reconstructed power spectrum, translate to significantly sharper and less shifted BAO features (right column). In the quadratic bias contributions (middle panels), reconstruction can be seen to dampen the amplitude of the BAO feature in the $b_{2}$ and $b_{1}b_{2}$ contributions, which “wiggle” in Fourier space, while leaving the spectrally smooth $b_{2}^{2}$ contribution essentially intact. Since the BAO feature in the quadratic bias contributions will tend to smear and shift the observed BAO peak from its linear theory position, reconstruction serves to remove these confounding nonlinearities as expected. The shear terms have less pronounced (i.e. smoother) features at the BAO scale—which we will show in Section 8 as being essentially in-phase with the linear theory oscillations— that are less affected by reconstruction.

Finally, as noted in the discussion below Equation 2.4, the exponentiated $A_{ij}$ in Zeldovich power spectra are assumed to be long wavelength, IR modes which can be resummed while contributions from the rest of the shorter modes are perturbatively expanded. These expanded modes thus carry also a UV (small-scale) sensitivity that should be renormalized by adding the appropriate counterterms, quadratic in wavenumber and proportional to the Zeldovich power spectrum: $\alpha_{ab}k^{2}P_{\rm Zel}^{ab}$. In principle, we expect such counterterms in all three pieces of our reconstructed power spectrum, however given that the $P^{ss}$ consists of mostly IR modes we expect its counterterm contribution to be suppressed relative to similar terms in $P^{dd}$ and $P^{ds}$, though it could still be non-vanishing due to contributions we neglected in approximating Equation 3.1. While the counterterms $\alpha_{dd}$ and $\alpha_{ds}$ are highly nondegenerate due to the different supports of $P^{dd}$ and $P^{ds}$ (see Figure 2) in $k$-space, in this work we will also explore modelling the reconstructed power spectrum using only one counterterm, $\propto k^{2}P_{\rm Zel}$, for both $P^{dd}$ and $P^{ds}$ contributions, since such a contribution would also be degenerate with any potential derivative biases (see e.g. ref. [36]). We will return to the difference between these options in Section 7.

4.2 Redshift space

In this section we develop analytic expressions for the redshift-space reconstructed power spectrum in both Rec-Sym and Rec-Iso. Methods recently developed in ref. [50] allow us to extend the LPT redshift-space power spectrum calculation to include bias and the specifics of reconstruction, which we summarise here and present in detail in Appendix C. As we will show shortly, Rec-Sym and Rec-Iso are not equivalent even to linear order. Specifically, we have

[TABLE]

i.e. while Rec-Sym restores supercluster infall at linear order, Rec-Iso removes redshift-space distortions at large scales while keeping them at small scales666We have amended Equation 4.11 to correct for a typo pointed out in ref. [51], wherein the $b\mathcal{S}$ in Equation 4.11 was missing a factor of $b$. No other results or conclusions are affected.. As we will see, this produces a smooth modulation in the broadband power nondegenerate with the BAO wiggles.

Since both the smoothed and displaced fields are uniformly multiplied by $R_{ij}$ in Rec-Sym, it is straightforward to calculate the reconstructed power spectrum using Equation 2.4 with

[TABLE]

In particular the angular structure of the q integral follows as in the calculation of the galaxy power spectrum without further modifications, and the set of bias terms in the $dd$, $ds$ and $ss$ spectra are identical to the real space case. The reconstructed power spectrum can then be calculated as one would the unreconstructed redshift space power spectrum. We develop the formalism to do the latter in Appendix C.3 and comment on the changes required to go to the reconstructed case therein.

The cross spectrum in Rec-Iso is slightly different since only the displaced field is multiplied by the redshift space transformation, $R_{ij}$. The displaced and shift fields in this case are thus instead

[TABLE]

Since the displaced and shift moves thus lie in redshift and real space, respectively, their auto spectra can also respectively be calculated as in Rec-Sym and real space reconstruction; however, the cross spectrum is only “half transformed” into redshift space and thus requires special attention. The exponentiated two-points displacements are given by

[TABLE]

such that the zero-lag piece due to the displaced-displaced correlation is fully transformed into redshift space, the zero-lag piece due to the shifted-shifted correlation is untransformed, and the coordinate dependent displaced-shifted correlation is“half transformed.” In particular, defining as usual $\textbf{q}=q\ \hat{q}$ and $\hat{k}\cdot\hat{q}=\mu$, the last piece is

[TABLE]

where we have defined the tilded quantities without the usual zero lag piece777For notational simplificity, the functions $X$ and $Y$ are always defined in real space.

[TABLE]

Note that $2\tilde{Y}=-Y$ since $Y$ does not posess a zero-lag piece. The azimuthal-angle dependence in $\hat{q}\cdot\hat{n}$ will require us to do the integral (Appendix C.2)

[TABLE]

where we have defined

[TABLE]

Note that the $\Gamma$ functions in the denominator will kill any terms in the sum for which $2m-\ell$ is negative, such that the sum really only contains $\ell/2$ terms and is always convergent in $A$. The full cross spectrum is then given by

[TABLE]

where $A=k^{2}f\nu\sqrt{1-\nu^{2}}\tilde{Y}^{(ds)}$ and we have defined $\Sigma^{2}=\tilde{X}(0)$ and $\alpha_{0}=1+f(f+2)\nu^{2}$. The remaining integrals can then be performed using the usual tricks for powers of $\mu$ using the series described in Appendix C.1, and are explicitly given at the end of Appendix C.2.

Figures 4 and 5 show the various bias contributions to the reconstructed redshift space power spectrum monopoles and quadrupoles within Rec-Sym and Rec-Iso, respectively. A significant difference between the two methods can be seen by comparing the matter (i.e. “1”) pieces in the top panels of the two figures. While all three linear bias contributions to the reconstructed power spectrum monopole ($\propto 1$, $b_{1}$, $b_{1}^{2}$) approach the Zeldovich monopole in the large scale limit in Rec-Sym, the matter contribution to the Rec-Iso monopole instead approaches the $b_{1}^{2}$ contribution, which does not receive redshift space distortions in the linear theory limit. This is because the power spectrum at the largest scales is dominated by the autospectrum of the un-redshifted shift field, $P^{ss}$. While the matter and $b_{1}$ contributions to the reconstructed quadrupole approach linear theory in Rec-Sym, they vanish on large scales in Rec-Iso. On the other hand, the majority of the higher bias contributions (excluding $b_{2}$ and $b_{s}$) are sourced only by $P^{dd}$ and are thus identical between the two methods, as can be seen by comparing the lower two rows of Figures 4 and 5. This corresponds to our intuition that redshift-space distortions are less prominent for highly biased tracers, and that the differences between Rec-Iso and Rec-Sym disappear if we remove RSD. In addition, the contributions enumerated above are supplemented by counterterms $(\alpha^{\ell}_{dd},\alpha^{\ell}_{ds},\alpha^{\ell}_{ss})$, where we need a separate counterterm for each pair and multipole as discussed below Equation 2.7, though as in the real space case we also explore the possibility of only fitting one counterterm each for the net reconstructed monopole and quadrupole.

5 Reconstructed correlation function

The configuration space two-point function (the correlation function) can be obtained from our Fourier-space results by Hankel transform. It is also possible to rewrite the $q$-dependent integrals to compute $\xi(r,\nu_{r})$ directly, where $\nu_{r}=\hat{n}\cdot\hat{r}$. Here we reprise the calculation of ref. [16], extending it to include the additional bias terms and commenting explicitly on several numerical issues which arise. We have checked that our Fourier and configuration space results agree numerically to significantly sub-percent levels in both real and redshift space for both Rec-Sym and Rec-Iso.

The general formula for the cross spectrum of two tracers $a$ and $b$ given in Equation 2.4 can be Fourier transformed to give [37, 36, 38, 42]

[TABLE]

where we have defined

[TABLE]

and placed the superscript $ab$ in $A^{ab}$ into the subscript for notational convenience. In the configuration space calculation above, the Lagrangian two-point functions (e.g. $A_{ij}$, $U_{i}$, $\Upsilon_{ij}$) can be computed using the formulae provided in Appendix A. The above formula can be translated into redshift space by multiplying the Lagrangian two-point functions with vector indices by the appropriate factors of $R_{ij}=\delta_{ij}+f\hat{n}_{i}\hat{n}_{j}$. Taking the line-of-sight to be in the $z$ direction without any loss of generality, this is equivalent to multiplying by the matrix $\text{diag}(1,1,1+f)$. When calculating the un-reconstructed redshift space correlation function, this multiplication is equivalent to multiplying each $z$ component index of vector and tensor quantities (e.g. $U^{a}_{z}$ or $A^{ab}_{yz}$) by $1+f$, and dividing the corresponding components in the matrix inverse, $A_{ab}^{-1}$, by the same factor. The redshift-space counterterm $\alpha_{2}k^{2}\nu^{2}$ can be included in the correlation function by adding $\alpha_{2}\hat{n}_{i}\hat{n}_{j}G_{ij},$ which similarly is equivalent to $\alpha_{2}G_{zz}$ when picking $z$ as the line-of-sight direction.

The reconstructed correlation function in real and redshift space can be calculated using Equation 5.1 by defining “displaced” and “shifted” tracers as in the case of the power spectrum (Sections 4.1 and 4.2) and calculating the combined quantity $\xi_{\rm recon}=\xi_{dd}+\xi_{ss}-2\xi_{ds}$. For reconstruction using Rec-Sym, the same shortcuts of multiplying by factors of $1+f$ in lieu of matrix multiplication and inversion apply, since all vector and tensor quantities undergo the same transformation by $R_{ij}$. The calculation for Rec-Iso is more complicated. As was the case in Fourier space, the displaced-displaced and shifted-shifted auto-correlation functions are equal to their counterparts in Rec-Sym and real-space reconstruction, while the displaced-shifted cross-correlation function contains a mix of real and redshift space factors. In particular, from Equation 4.14 we see that the two zero-lag pieces and one $q$-dependent piece of $A^{ds}_{ij}$ in Rec-Iso are independently transformed by different numbers of $R_{ij}$’s. For this reason, when calculating the correlation function in Rec-Iso, the matrix inverse of $A^{ds}$ in redshift space cannot be simply obtained by dividing the real space inverse by factors of $1+f$; rather, the uninverted matrix must be redshifted piece by piece as in Equation 4.14 and then inverted numerically (we use Cholesky decomposition).

6 Other statistics

While the correlation function and power spectrum are the most frequently considered 2-point functions, there are other variants that have some advantages. Since these can all be written in terms of the correlation function or power spectrum, our model provides a consistent prediction for them as well. Of particular interest for BAO is the $\omega_{\ell}$ statistic of ref. [52], which combines the scale-localization of the Fourier-space methods with the compactness and easy treatment of masks of the configuration-space methods.

In principle $\omega_{\ell}$ can be calculated from either the configuration-space or Fourier-space expressions given above, but we have found it more convenient to start from the Fourier expressions. Since these are computed using FFTlog they naturally cover a very wide range of $k$, making the transforms to $\omega_{\ell}$ easy to implement. For example

[TABLE]

with $\widetilde{W}_{0}$ given in ref. [52] (see their Fig. 1 and Appendix A). At large scales $\widetilde{W}_{0}\propto k^{2}$ while at small scales $\widetilde{W}_{0}\propto k^{-4}$. Our formalism naturally provides predictions for $\omega_{\ell}$ using the same set of bias and nuisance parameters as for $\xi_{\ell}$ and $P_{\ell}$.

7 Comparison to N-body

To look at the domain of validity of our analytic results we compare to the DarkSky N-body simulation suite888http://darksky.slac.stanford.edu, specifically simulation ds14_a [1]. This simulation used the 2HOT code [53] to evolve $10240^{3}$ particles in an $(8\,h^{-1}$Gpc$)^{3}$ volume to model the growth of structure in a $\Lambda$CDM cosmology with $\Omega_{M}=1-\Omega_{\Lambda}=0.295$, $h=0.688$, $n_{s}=0.968$ and $\sigma_{8}=0.835$. Initial conditions were generated from a glass using $2^{\rm nd}$ order Lagrangian perturbation theory at $z=93$. Halos were found using the Rockstar code [54]. We extracted the positions, velocities and masses of halos more massive than $M_{200b}=10^{12}\,h^{-1}M_{\odot}$ from the publicly available data at $z=0$ (data at higher $z$, which would have been a more relevant comparison, were not available). We computed the halo correlation functions and power spectra, in real and redshift space. For the redshift-space quantities we assumed the plane-parallel approximation with the line-of-sight being the $z$-axis. We also obtained the linear theory power spectrum used to generate the initial conditions, which we take as the input to our model.

We implemented the algorithm described in §3 using the periodicity of the box and FFTs to perform the smoothing and computation of the shifts. As for the power spectrum and correlation function, the plane-parallel approximation with line-of-sight the $z$-axis was assumed for the redshift-space quantities. The code takes as input an assumed large-scale bias, $b$, and growth parameter, $f$, in addition to a Gaussian smoothing length, $R$. We used the $b$ obtained from the ratio of the linear theory and real-space halo power spectra at low $k$ (see Table 1), and $f\simeq 0.508$ appropriate to the simulation cosmology at $z=0$, and note in passing that the goodness-of-fit of our results did not seem to be greatly improved by substituting the linear bias thus obtained with the value of $1+b_{1}$ obtained by fitting the pre-reconstruction data with our model up to quasi-nonlinear scales.

We computed the reconstructed field in both real and redshift space. In each case the shifted and displaced positions were computed using a $2048^{3}$ FFT, which resolves the (Gaussian) smoothing length by $2.5-5$ grid cells for $R\simeq 10-20\,h^{-1}$Mpc. We used as many “random” positions as halos in each case, for simplicity, and computed the power spectra and correlation functions for $dd$, $ds$ and $ss$ assuming periodic boundary conditions. The reconstructed power spectrum or correlation function can then be computed as $dd-2ds+ss$, and we can look at each of the contributions separately. Note that our choice of equal numbers of randoms and data points means the shot noise on the reconstructed power spectrum is twice that of the pre-reconstructed field.

We compare the N-body results to our model with $b_{1}$ and $b_{2}$ and include the minimal set of counterterms as described in the preceding sections (one and three pre- and post-reconstruction, respectively, in real space) as well as a constant shot noise component fit to the data. For brevity our discussion will focus on halos with masses between $12.5<\log_{10}(M/h^{-1}M_{\odot})<13.0.$, though we obtained qualitatively similar results in the lower and higher mass bin as well, and show fits of the reconstructed redshift space power spectrum in the latter at the end of this section. We have checked that including nonzero shear bias $b_{s}$ does not visibly improve the goodness-of-fit. The top-left pair of panels of Figure 6 compares the unreconstructed real-space power spectrum in our model with $(b_{1},b_{2})=(0.02,-0.8)$ with that in DarkSky. The quadratic bias, $b_{2}$, accounts for a non-negligible fraction of the total power at essentially all scales and significantly reduces the constant shot noise term in the fit. We find that with a counterterm $\alpha\approx 11\ h^{-2}\text{Mpc}^{2}$ our model agrees with the data at the percent level out to $k\simeq 0.4\,h\,{\rm Mpc}^{-1}$. The counterterm accounts for roughly a $10\%$ correction at $k=0.1\,h\,{\rm Mpc}^{-1}$, and it is worth noting that even in its absence our model accurately captures the BAO features in the power spectrum, as evidenced by the lack of oscillatory features in the fit residuals.

The remaining panels of Figure 6 show the fit for the reconstructed power spectra at three smoothing scales $R=10$, $15$, $20\,h^{-1}$Mpc. We have tested whether the data could be reproduced using only one counterterm, $\alpha$ (shown in orange), or equivalently from one derivative bias $b_{\nabla^{2}}$, and find that such a choice dramatically reduces the range-of-validity of the model compared to three counter terms. While we adopted a rather conservative approach in fitting these data, prioritizing the accuracy of our predictions at low $k$ rather than producing reasonable-looking fits to smaller scales, our model with three counterterms $(\alpha_{dd},\alpha_{ds},\alpha_{ss})$ nonetheless reproduces both the broadband power and oscillatory features of the reconstructed power spectrum out to $k=0.2\,h\,{\rm Mpc}^{-1}$ at the percent level for $R=15$ and $20h^{-1}$ Mpc. That each of the three constituent spectra in $P^{\rm recon}=P^{dd}+P^{ss}-2P^{ds}$ has distinct short-wavelength behavior and $k$-space supports underlies the success of our model with three counterterms—each of which has highly nondegenerate scale dependence—versus the one-counterterm alternative. We have found that setting $\alpha_{ss}=0$ does not qualitatively alter the degree to which our model fits the data; we have made this choice in all of our fits below, but note that as $\alpha_{ss}k^{2}$ vanishes quadratically towards low $k$, the data are also naturally rather insensitive to it. Indeed, since nonlinear corrections are typically of order $k^{2}\Sigma^{2}$, and the smoothing scale is chosen such that $(\Sigma/R)^{2}\ll 1$, the insensitivity of $P^{ss}$ to these corrections follows almost by construction. However, a bump-like feature around $k=0.1\,h\,{\rm Mpc}^{-1}$ is persistent across all the fits, peaking at less than half a percent when $R=20\,h^{-1}$Mpc and growing to a full percent at $R=10\,h^{-1}$Mpc. The appearance of such a feature, growing towards smaller smoothing scales, is consistent with our neglect of nonlinear corrections to the smoothed displacements, which should increase towards smaller smoothing scales roughly as $\boldsymbol{\Psi}/R$; we discuss one such nonlinearity in Appendix E. For sufficiently small smoothing scales, even the assumption that the smoothing of the BAO feature can be essentially captured with resummed linear displacements $\boldsymbol{\Psi}^{d,s}$ will break down, and indeed our fit residuals begin to show noticeable oscillatory behavior at the smallest smoothing scale shown ($R=10\,h^{-1}$Mpc). At $R=15\,h^{-1}$Mpc and in the sample variance limit with Gaussian errors, the feature at $k=0.1\,h\,{\rm Mpc}^{-1}$ should be detectable with $\chi^{2}=V_{\rm obs}/(2h^{-3}\ \text{Gpc}^{3})$, where $V_{\rm obs}$ is the total observed volume. If we were to instead smooth using the larger $R=20\,h^{-1}$Mpc, the $\chi^{2}$ is roughly halved. For such a smoothing this feature represents a $\chi^{2}$-penalty of $0.2$ for a sample variance limited survey of $14\,000\,\mathrm{deg}^{2}$ covering $0\leq z\leq 0.3$, and would be slightly smaller for finite number density.

The pre- and post-reconstruction real-space correlation functions can be directly compared by computing the Fourier transforms of the above fits. However, in comparing our theory with DarkSky we found that the $z=0$, pre-reconstruction halo power spectra all have significant excess power at low $k$ compared to the predictions of linear theory with scale-independent bias. The origin of this excess is unclear, and is not addressed in ref. [1]. It appears to arise from a significant number of low $k$ modes, and so is unlikely to be simply a statistical fluctuation in the initial conditions. It shows up in all of our halo samples, and is highly correlated among mass bins. This excess power is small for modes to the right of the power spectrum peak and probably has only a small impact on the dynamics on BAO scales. In Fourier space we simply confine our fitting and modeling to $k>0.01\,h\,\mathrm{Mpc}^{-1}$. In configuration space, however, the additional long-wavelength power slightly distorts the shape of the BAO peak, and to enable a fair comparison we have added appropriate long-wavelength modes to our theoretical predictions assuming linear theory; specifically, we find that the fitting form $P_{\rm lw}(k)=A\ (k/k_{0})^{n}$, where $A=3.5\times 10^{4}h^{-3}\ \text{Mpc}^{3}$, $k_{0}=10^{-3}\,h\,{\rm Mpc}^{-1}$ and $n=-1.7$, describes well both the long-wavelength excess seen in the power spectrum below $k<0.01\,h\,{\rm Mpc}^{-1}$ and dramatically improves the agreement between the unreconstructed correlation function in theory and DarkSky. The contribution to the pre- and post-reconstruction power spectra and correlation function of these long wavelength modes is shown in Figure 7. Without the long-wavelength correction, the DarkSky results do not agree with theory on the large scales to the right of the BAO peak, which should be well-described within linear theory, nor in the BAO “dip,” both pre- and post-reconstruction. Due to the ad-hoc nature of our correction, in the remainder of this section we will focus our comparisons on Fourier space, wherein long-wavelength modes must decouple. However, we caution that small, localized features in Fourier space can cause extended distortions in configuration space where data points are highly correlated. In Figure 8, we show the effect of the $k=0.1\,h\,{\rm Mpc}^{-1}$ bump described in the previous section by additively “filling” it with a small, localized Gaussian profile, as shown in the left panel. The effects of this bump, Fourier-transformed, are shown in the right panel: while sub-percent in Fourier space, the $k\simeq 0.1\,h\,{\rm Mpc}^{-1}$ feature gives rise to visible distortions to the BAO feature in configuration space.

Finally, fits for the pre- and post-reconstruction power spectra in redshift space are shown in Figure 9. We have chosen to summarize the angular dependence of the redshift-space power spectrum in terms of its monopole and quadrupole, though our model predicts the full $P(k,\mu)$ and higher multipoles as well. As in real space, we have fitted for the bias parameters $(b_{1},b_{2})$ using the unreconstructed data and applied the same set of bias parameters to predict the power spectra in both Rec-Sym and Rec-Iso. We adopt the full set of six counterterms, three each $\alpha^{\ell}_{dd},\alpha^{\ell}_{ds},\alpha^{\ell}_{ss}$ for the monopole ($\ell=0$) and quadrupole ($\ell=2$), but also explore the possibility of utilizing only one counterterm $\alpha^{\ell}$ per multipole (corresponding to a derivative bias for both the halo density and velocity). In all cases, our base model with six counterterms fits the data at the percent level or below past $k=0.2\,h\,{\rm Mpc}^{-1}$ in both the monopole and quadrupole moments. Notably the Zeldovich approximation produces oscillation-free residuals even in the absence of counterterms (green), with the counterterms providing a physics-based broadband model ($\sim\alpha_{ab}^{\ell}k^{2}P_{\ell,ab}$) that reproduces the N-body results at the percent level. Our fits do not explicitly include nonlinear redshift space distortions such as fingers-of-god, though such effects are perturbatively accounted for by velocity counterterms to lowest order. For completeness, in Figure 10 we show the same fits for the mass bin $13.0<\log(M/h^{-1}M_{\odot})<13.5,$ where our model fits the data at percent level over a similar range of scales using the parameters $(b_{1},b_{2})=(0.23,-1.0).$

Lastly, let us comment on the comparison fits in pre- and post- reconstructed cases. Given that our shift field, $\chi$, is constructed only from long-wavelength modes explicitly isolated from observed field, $\delta$, by filtering out the nonlinear scales larger than $k\gtrsim 1/R$, we have no reason to suppose that the perturbative structure of our results will significantly change. In other words, by performing the mapping in Equation (4.1), we have reconstructed only the long modes, thereby reducing nonlinear smoothing due to large scale (infrared) modes, while the bulk of the small-scale nonlinear modes, as well as FoG effects, should remain unreduced. In addition, Lagrangian perturbation theory (PT) conveniently separates nonlinearities due to long and short modes, exponentially resumming the former while expanding the latter order-by-order [39, 38]. Because of this, we do not expect dramatically different PT behavior in the pre- and post-reconstructed results. These arguments are also supported by Figures 9 and 10, which show our model exhibits quantitatively similar degrees of fit pre- and post reconstruction.

8 Comparison to earlier work

There has been significant theoretical activity in modeling post-reconstruction clustering (see references in the introduction). Our framework encompasses most of these previous perturbation theory expressions when appropriate approximations and phenomenological choices are accounted for. To the best of our knowledge, the framework presented here captures for the first time all of the relevant post-reconstruction effects and is unique in accurately handling both Fourier and configuration space results, in real and redshift space and includes all the bias operators to quadratic order.

Not all models are based on perturbation theory calculations however, and many phenomenological models have been introduced in order to describe the post-reconstruction statistics. Restricting ourselves just to models of the ‘standard’ reconstruction algorithm [2], §3.1 of ref. [10] discusses early models (which were of the form $P(k)=B(k)P_{\rm lin}(k)+A(k)$ with $B(k)$ and $A(k)$ smooth functions). Starting with the first applications to data in ref. [10] the form used to fit reconstructed power spectra is based upon a split between a “smooth” and “wiggle” contribution to $P(k)=P_{\rm nw}(k)+\Delta P_{\rm w}(k)$, with a phenomenological damping of the wiggle component motivated by perturbation theory [8]. In ref. [10] the parameters of the model were fit to N-body simulations, and this has become common. This approach has dominated the modeling of observations to date (e.g. refs. [55, 56, 57] for recent examples) though ref. [58] is an example of an analysis that did not take this approach. However, we note that the choice of the wiggle/no-wiggle split exhibits a certain amount of freedom in the separation of the wiggle and broadband part. This of course implies that, in order to extract accurate information from the e.g. BAO, either both wiggle and broadband part have to be modeled to the same level of accuracy, or the extracted wiggle part from the data needs to exactly correspond to the model (see also refs. [59, 60] for related discussion). The latter requirement, even though implicitly assumed in most of the current BAO treatments, is rarely subject to performance checks and scrutiny. In this context, it is also worth noting that the common choice of $P_{\rm nw}$ derived in ref. [61] does not fully capture the broadband linear power spectrum at the precision attained by modern Boltzmann codes. Figure 11 shows three possible linear wiggle power spectra, based on no-wiggle spectra computed using the fitting formula from ref. [61], B-splines [62] or a Savitsky-Golay filter; even the latter two, which agree asymptotically with the full linear theory power spectrum, exhibit noticeably different oscillatory behavior. This indicates that extracting the corresponding wiggle spectra from the data is a challenging and sensitive step which can, on the other hand, be avoided if the broadband is included in the theoretical framework. Models phenomenologically relying on a wide separation of scale, assuming scale-independent bias or sufficient smoothness that could be accounted for by nuisance parameters such as $A(k)$ above, might suffer from overall systematic offsets. Finally, it is also often the case that the nuisance parameters and BAO scaling parameters are not consistent between the configuration-space and Fourier-space analyses (i.e. the two do not form a Fourier transform pair) which could prove problematic if fits in both spaces are combined.

By contrast the Zeldovich calculation above gives a consistent framework for understanding the nonlinear smoothing of the BAO feature, both pre- and post-reconstruction, in both configuration and Fourier space. Roughly speaking, the Gaussian smoothing kernel in the empirical model is replaced by a Lagrangian coordinate-dependent kernel $\exp[-k_{i}k_{j}A_{ij}(\textbf{q})/2]$. One might thus hope to formally extract the model for wiggle-only part as an approximation to the calculation presented in the main body of this paper; indeed, such a calculation was performed in ref. [62] and extended to terms involving linear bias, redshift space distortions and reconstruction in ref. [26]999We note that redshift-space reconstruction model presented in ref. [26] contains phenomenological damping factors that do not capture the exact behaviour of $X_{ds}$ term given by in Equation (4.9). We repeat this calculation and derive the proper damping factors for the wiggle component in Appendix D.. Figure 12 compares the results of our full Zeldovich calculation in the Rec-Sym scheme, with the broadband subtracted out by calculating the corresponding Zeldovich power spectrum using the no-wiggle power spectrum, versus the resummed linear wiggle power spectrum (RWiggle; using the proper exponential damping dependencies given in Appendix D), for the same linear bias values and with all higher bias terms set to zero. The two are in excellent agreement, especially in the case of the reconstructed power spectrum, with RWiggle slighly underdamping the BAO wiggles towards small scales compared with the full Zeldovich calculation for the unreconstructed power spectrum.

However, even though RWiggle and the full Zeldovich calculation exhibit a high level of agreement on the shape of the wiggle component, the RWiggle method depends on the separation procedure of the wiggle and broadband components while the full Zeldovich calculation requires no such steps. Specifically, the Zeldovich calculation deals only with the combination $P^{\rm Zel}_{\rm w}+P^{\rm Zel}_{\rm nw}$, which is obviously invariant under the split, while RWiggle models only the split-dependent $P^{\rm Zel}_{\rm w}$. This implies that in order to use RWiggle in practical analyses either the broadband part needs to be modeled to equally high accuracy or a highly accurate wiggle extraction procedure is needed in order to guarantee feasible comparison of theoretical model and the data. The latter seems to be a challenging task, potentially subject to systematic offsets and bias. On the other hand, the resulting differences in the wiggle spectrum should still be broadband and could be fit away using nuisance parameters using sufficiently general broadband models.

Finally, our model differs from most in the literature in taking into account higher bias terms such as $b_{2}$ and $b_{s}$, allowing us to assess systematic effects introduced by assuming scale-independent bias. These higher biases can contribute both significant broadband power (e.g. the top-left panel of Figure 9) and modulate the phase and amplitude of BAO oscillations through mode-coupling effects [13]. However, explicit calculation shows that the latter effect is only noticeable at very high values of bias. Figure 13 shows the effects on the wiggle component of adding nonzero quadratic density and shear biases $b_{2}$, $b_{s}$, for bias values $(b_{1},b_{2})\approx(5,20)$ chosen according to the peak-background split (PBS) on a Press-Schechter mass function [63], and assuming $b_{s}\approx b_{2}$, as compared to RWiggle. The quadratic density bias, $b_{2}$, induces an apparent phase shift towards large $k$, and can be seen to be essentially out-of-phase with the linear BAO wiggles; however, these out-of-phase contributions are dramatically reduced by reconstruction. By contrast the shear bias, $b_{s}$, produces oscillatory features roughly in-phase with the linear theory contributions and is largely unaffected by reconstruction. For completeness, we have also plotted the potential oscillatory contribution of a derivative bias, $b_{\nabla^{2}}$, which modulates the overall amplitude of the power spectrum and is degenerate with the various counterterms, $\alpha_{ab}$.

To investigate the extent to which the broadband and oscillatory contributions of higher bias terms can be mitigated by a suitable broadband model, we conducted an exploratory “fit” of the redshift-space monopole and quadrupole pre- and post-reconstruction in the case where the truth is given by the Zeldovich approximation including nonzero $b_{2}$ and $b_{s}$ but fit by an empirical model with only $b_{1}$, an isotropic BAO scale paramter $\alpha_{\rm BAO}$ and polynomial broadband contributions of the form employed in ref. [64] before reconstruction. Specifically, we assume an empirical model of the form

[TABLE]

both pre- and post-reconstruction, where $P_{l,b_{1}}$ denotes redshift-space multipoles in the Zeldovich approximation with all higher biases set to zero. For this exercise we assumed a sample variance limited survey at $z=0$ and $z=1.2$ with Gaussian covariances between the monopole and quadrupole, fit up to $k_{\rm max}=0.25\,h\,{\rm Mpc}^{-1}$ and note that the results are independent of survey volume. In Figure 14 we have plotted the resulting shifts in the BAO scale assuming PBS values for $b_{1}$ and $b_{2}$, taking values for $b_{s}$ as a function of $b_{1}$ from ref. [65]. At $z=0$, we find that neglecting higher biases in favor of the empirical model induces shifts of less than half a percent in the BAO scale over a wide range of halo masses both pre- and post-reconstruction, though reconstruction more than halves the forecasted shift for essentially all values of bias surveyed (Figure 14). At $z=1.2$ the shifts are further reduced, amounting to less than a tenth of a percent across a wide range of bias values prior to reconstruction and essentially vanishing post reconstruction. These shifts would be well-within the margin of error of both current and next-generation surveys like DESI [66], especially post-reconstruction, suggesting that nonlinearities (e.g. higher bias) in the power spectrum should not hinder accurate recovery of the BAO signal. On the other hand, the value of the linear bias, $b_{1}$, was significantly affected by the choice of broadband model, with fits from the empirical model deviating from the true value by more than five percent in many cases.

9 Conclusions

Baryon acoustic oscillations (BAO) are an important probe of fundamental physics and a prime focus of upcoming surveys such as DESI [66] and EUCLID [67]. The BAO features act as a “standard ruler” whose cosmological evolution is largely immune to astrophysical effects but whose signal-to-noise ratio is lowered by nonlinear structure formation. BAO reconstruction attempts to sharpen the BAO signal by removing some of the nonlinear smearing due to large scale displacements [2]. In this paper we develop an analytical model, within the Lagrangian perturbation theory framework, to study the algorithm for density-field reconstruction proposed in ref. [2]. Linear Lagrangian perturbation theory (the Zeldovich approximation) provides an excellent description of these nearly linear displacements and BAO smoothing pre-reconstruction [28, 37], making LPT a promising arena within which to model the effects of reconstruction.

We develop a self-consistent framework with which to calculate the two-point statistics of galaxies, employing a consistent set of parameters to fit the power spectrum and correlation functions, pre- and post-reconstruction in real and redshift space. The broad validity of such LPT models allows for joint fits to the pre- and post-reconstruction two-point statistics enabling e.g. a fit for redshift-space distortions and the linear growth rate, $f\sigma_{8}$, simultaneously with Alcock-Paczynski distortions constrained by BAO analyses [16]. Based on ref. [50], we derive explicit formulae, to calculate the redshift-space power spectrum within the Zeldovich approximation, both pre- and post-reconstruction, as an infinite series of spherical Bessel transforms. Our model updates the developments for the reconstructed correlation function in ref. [16], and is – as far as we are aware – the first model of reconstruction to include a consistent set of bias terms up to quadratic order, including shear and derivative biases. We show that the oscillatory behavior induced by the quadratic density bias, $b_{2}$, are out of phase with the linear BAO feature and greatly reduced post-reconstruction, while those due to the quadratic shear bias, $b_{s}$, are in-phase and essentially unchanged. In addition, we show that each multipole moment of the reconstructed power spectrum should be, to lowest order, corrected for by a set of three counterterms each, which perturbatively correct both nonlinear smoothing and broadband power.

We compare our analytic predictions with N-body data from the DarkSky simulation [1] at $z=0$, focusing on halos between $12.5<\log_{10}(M/M_{\odot})<13.0$. Our base model, involving only $b_{1}$ and $b_{2}$ and appropriate counterterms, jointly fits the pre-reconstruction real-space power spectrum and redshift-space monopole out to $k=0.4\,h\,{\rm Mpc}^{-1},$ and the quadrupole out to $k=0.2\,h\,{\rm Mpc}^{-1},$ reproducing the oscillatory BAO wiggles in the data with high fidelity. Our model with the same bias parameters performs equally well in configuration space around the BAO scale, though we found it necessary to correct for a large excess in large-scale power encountered in the DarkSky data. Utilizing the same values for the bias parameters but allowing counterterms to vary, we find that our model performs similarly in real space post-reconstruction for smoothing scales $R=15$ and $20\,h^{-1}$Mpc, reproducing both the oscillatory features and broadband past $k=0.2\,h\,{\rm Mpc}^{-1}$, but fails to reproduce the oscillatory features when $R=10\,h^{-1}$Mpc, likely due to the fact that we have worked to lowest order and at $z=0$ displacements on that scale are significantly nonlinear. We point out a less severe feature in the residuals at $k=0.1\,h\,{\rm Mpc}^{-1}$ that diminishes with larger smoothing scales which we believe arise from higher order terms and caution that neither our calculation nor the standard reconstruction algorithm take these into account. A more complete, iterative reconstruction scheme (e.g. ref. [20]) may reduce these features. The modeling of these nonlinearities, and possible remedies, are beyond the scope of this paper, but as an exploratory example we calculate the effects of one possible nonlinearity due to the mapping between Eulerian and Lagrangian coordinates in Appendix E.

Our model also predicts the multipole moments of the redshift-space power spectrum and correlation functions in both of the redshift-space schemes (Rec-Sym and Rec-Iso) we consider. This is critical in order for it to be applied to data, since the most constraining BAO measurements are performed in redshift space. The model provides a good fit to the monopole and quadruople moments of $P(k)$ measured in DarkSky in both the Rec-Sym and Rec-Iso schemes for smoothing scales of $R=15\,h^{-1}$Mpc or larger (at $z=0$). Again, for smaller smoothing scales the Zeldovich model differs from the N-body results (as expected). These effects would be smaller at higher redshift, where the theory is more likely to be applied.

Finally, there exists an extensive literature studying the modeling of reconstruction and the BAO signal, and we compare our model to several existing alternatives. One popular technique, based on ref. [8], is to separate the power spectrum into a smooth “no wiggle” component and an oscillatory “wiggle component,” and to damp the latter by an exponential factor fit to simulations while supplementing the former with a polynomial in wavenumber to fit the broadband power. This technique can be more rigorously derived as a particular resummation of the nonlinear contribution of long-wavelength modes much like our Zeldovich calculation itself [62, 26], in which case the damping parameters can be derived theoretically. When the “wiggle” components are isolated we find that the latter is in excellent agreement with our Zeldovich calculation, particularly after reconstruction. In Appendix D we re-derive the IR-resummed “wiggle” power spectrum (RWiggle) directly within our Zeldovich framework, updating the exponential damping for the cross term $P^{ds}$. We highlight that our Zeldovich framework naturally encompasses broadband effects, while methods depending on wiggle/no-wiggle splitting might be subject to additional systematic offsets and biases. These could originate from the fact that the wiggle/no-wiggle splitting is not unique, and thus relies on correctly predicting the broadband or extracting the corresponding wiggle part from the data to high accuracy. On the other hand, the Zeldovich framework correctly captures broadband power over a large range of scales in addition to reproducing the oscillatory features in the reconstructed power spectrum. In fits to N-body data, we show how counterterms correct the sharpness of the BAO feature and broadband power simultaneously and consistently. Moreover, our model goes beyond linear bias to include quadratic density and shear bias, which we show contribute oscillatory terms to $P(k)$ that vary independently in amplitude and phase.

We close by noting a few avenues for future work. An obvious extension of our model is to include nonlinearities arising both from gravitational clustering and the reconstruction itself (e.g. Appendix E). The former may be most easily included in the context streaming models [36, 50], wherein the real-space modifications due to reconstruction and those proportional to the growth rate $f$ can be separately treated as modifications to the statistics of the galaxy density and galaxy density-weighted velocities, respectively, and which in addition have the advantage of resumming biased contributions to redshift-space distortions as well as nonlinear redshift-space phenomena like fingers-of-god. It is, however, not a-priori obvious which type of nonlinearity will present the most significant corrections. Other fruitful avenues would be to investigate the impact of wrong parameters on reconstruction (e.g. refs. [27, 46]) or to update the present treatment to newer reconstruction techniques. Finally, one could investigate the utility of our model for upcoming surveys like DESI [66] or Euclid [67]. These surveys will operate at higher redshifts where our calculations should perform even better, and our model will be a natural arena in which to understand the effects of highly biased tracers and the effects of cosmic evolution (e.g. evolving $b$ and $\sigma_{8}$) on the BAO feature measured in broad redshift bins.

We have publicly released our codes for configuration101010https://github.com/martinjameswhite/ZeldovichRecon and Fourier111111https://github.com/sfschen/ZeldovichReconPk space reconstruction, with the hope that they will be useful to other researchers. We have checked that the Hankel transform of the Fourier space code agrees, term by term, with the configuration space code to better than 1%, except very close to zero crossings.

Acknowledgments

We thank Hee-Jong Seo and Florian Beutler for useful discussions. SC is supported by the National Science Foundation Graduate Research Fellowship (Grant No. DGE 1106400) and by the UC Berkeley Theoretical Astrophysics Center Astronomy and Astrophysics Graduate Fellowship. M.W. is supported by the U.S. Department of Energy and by NSF grant number 1713791. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. This work made extensive use of the NASA Astrophysics Data System and of the astro-ph preprint archive at arXiv.org.

Appendix A Cross-spectra correlators

In this appendix we give analytic expressions for the two-point functions required to calculate cross-spectra, which are slightly different from those required to calculate the auto-spectra more commonly seen in the literature.

The two-point function for the Lagrangian displacement between two species separated by Lagrangian distance q is given by

[TABLE]

where

[TABLE]

Note that for cross spectra $X^{ab}(q)$ does not in general vanish as $q\rightarrow 0$. Similarly we have

[TABLE]

where

[TABLE]

and $P^{am}$ is the linear theory cross spectrum between tracer $a$ and matter, and the corresponding expression for $U^{b}$ follows by direct substitution.

Finally, the non-scalar shear correlators are given by

[TABLE]

where the functions of $q$ are given by

[TABLE]

and

[TABLE]

where following refs. [37, 36] we have defined

[TABLE]

The remaining scalar shear correlators, $\zeta$ and $\chi^{12}$, are identical to those found in evaluating the auto-spectrum, and we refer readers to refs. [37, 36].

Appendix B The pre- and post-reconstruction Zeldovich propagator

In this appendix we give expressions for the normalized cross-spectrum between the initial and final or reconstructed field. This is essentially a correlation coefficient, though it is also referred to as the propagator [68]. Specifically we define $G_{a}(\textbf{k})=\langle\delta_{0}(-\textbf{k})\delta_{a}(\textbf{k})\rangle/\langle\delta_{0}(-\textbf{k})\delta_{0}(\textbf{k})\rangle$, within the Zeldovich approximation, which quantifies the extent to which a tracer field $a$ is (de)correlated with the initial density $\delta_{0}$, and apply our results to derive the reconstructed field. Our results generalize those in ref. [14] to include halo bias.

As defined, the propagator $G_{a}$ is a special case of the cross spectrum and can be evaluated using Equation 2.4 by assuming that the linear field $\delta_{0}$ is a tracer $b$ with displacement $\boldsymbol{\Psi}^{b}=0$ and bias functional $F^{b}=\delta_{0}$, such that any Lagrangian two-point functions involving the displacement $\boldsymbol{\Psi}^{b}$ (e.g. $U^{b}$) or higher biases (e.g. $b_{2}^{b}$) vanish identically. Unlike in the conventional case, however, $F^{b}$ does not have a zero order piece equal to unity— we can thus compute our result directly by taking the derivative of Equation 2.4 with respect to $b_{1}^{b}$ with the above assumptions. This gives

[TABLE]

where we have used that

[TABLE]

only receives “half” of the zero-point contribution c.f. the power spectrum (where $\langle\boldsymbol{\Psi}^{b}\boldsymbol{\Psi}^{b}\rangle\neq 0$). Note that all higher bias contributions vanish. The generalization to the redshift space field can be straightforwardly accomplished by multiplying by appropriate factors of $R_{ij}$ in the numerator, though we will focus on real space in this appendix as RSD introduce an equally important but parallel form of decorrelation into the problem.

From the above results, the reconstructed-field ($\delta_{\rm rec}=\delta_{d}-\delta_{s}$) propagator can be written as $G_{\rm rec}=G_{d}-G_{s}$, where

[TABLE]

where the various linear spectra are defined as in Equation 4.6. The real-space post-reconstruction propagator is then

[TABLE]

The expression for $G_{\rm rec}$ helps to quantify how much of the decorrelation between the initial conditions and the final field arises due to bulk motions, and the manner in which this can be restored by the standard reconstruction algorithm. Roughly speaking, reconstruction reduces the decorrelation from the full matter $\Sigma^{2}$ to $\Sigma^{2}_{dd}$ past the smoothing scale for the matter piece, with the correlation at low $k$ close to unity assuming that the damping due to $\Sigma_{ss}^{2}$ there is negligible.

Appendix C Integrals for redshift space distortions via direct Lagrangian expansion

In this section we describe how to perform the three-dimensional integrals that occur when calculating redshift space power spectra in the Lagrangian formalism.

C.1 Angular Integrals

Calculations in Lagrangian perturbation theory frequently require evaluating integrals of the form

[TABLE]

These integrals can be conveniently expressed as infinite sums of spherical Bessel functions [69, 50], e.g.

[TABLE]

where $U(a,b,z)$ denotes the confluent hypergeometric function of the second kind.

C.2 Direct Lagrangian Expansion: Mi

Calculating the power spectrum in redshift space within the Zeldovich approximation requires a few extra steps when compared to the calculation in real space due to the line-of-sight dependence of RSD. In the following two sections we extend the two methods presented in ref. [50], Mi and Mii, to include bias terms up to quadratic order. These two methods correspond, roughly speaking, to active and passive transformations in Fourier space via $R_{ij}=\delta_{ij}+f\hat{n}_{i}\hat{n}_{j}$, respectively. For the general case, we found Mii to be somewhat more convenient, and will therefore offer only a cursory description of Mi, except for a special case involving the displaced-shifted field cross spectrum in Rec-Iso.

In Mi, the (halo auto-) power spectrum in redshift space is given by

[TABLE]

where the superscript $s$ denotes tensor quantities transformed by $R_{ij}$, e.g. $U^{s}_{i}=R_{ij}U_{j}$. Each quantity in the above integral can be written in terms of scalar functions and dot products between three unit vectors ($\hat{q}$, $\hat{k}$ and $\hat{n}$) whose angular structure underlies redshift space distortions; these are, in particular,

[TABLE]

where $\phi$ is the azimuthal angle in a polar coordinate system where the zenith is given by $\hat{k}$ and the plane $\phi=0$ is spanned by $\hat{k}$ and $\hat{q}$. The effect of the transformation $\mathbf{R}$ can then be captured by how the usual tensor basis $\hat{q}_{i},\delta_{ij},$ etc. is affected:

[TABLE]

where we have defined $\gamma(\mu,\nu)=\sqrt{1-\mu^{2}}\,\sqrt{1-\nu^{2}}/\mu\nu.$ The azimuthal dependence of $A^{s}_{ij}$ requires us to calculate polar-coordinate integrals of the form [50]

[TABLE]

where $\alpha_{2}=2f\nu^{2}(1+f\nu^{2})$ and $\alpha_{3}=f^{2}\nu^{4}$. These can be calculated as analytic power series in $\mu$ by taking derivatives of the identity

[TABLE]

where

[TABLE]

and $M(a,b,z)$ are hypergeometric functions of the first kind. Each term proportional to $\mu^{n}$ can in turn be evaluated as a spherical Bessel transform as in Equation C.2.

A simplified but demonstrative example of this calculation can be found in the calculation of the displaced-shifted cross spectrum in redshift space reconstruction via Rec-Iso. The $H^{0}_{\ell}$ expansion in Equation 4.16 is essentially Equation C.8 in the limit where $\beta\rightarrow 0.$ To proceed from Equation 4.16, we can use the identities in C.1 and refactor the resulting double sum over $n$ and $\ell$ to get

[TABLE]

where the redshift-space kernels are given by

[TABLE]

and, as before, $A=k^{2}f\nu\sqrt{1-\nu^{2}}\,\tilde{Y}^{ds}$ and $B=k^{2}(1+f\nu^{2})\,\tilde{Y}^{ds}$. Deriving these kernels for the other terms is entirely analagous121212The mater contribution was given in ref. [50]..

C.3 Direct Lagrangian Expansion: MII

An equivalent approach to rotating each Lagrangian displacement $\boldsymbol{\Psi}$ by $\mathbf{R}$ is to instead passively transform the Fourier basis by $\mathbf{R}^{\rm T}$, such that the wavenumber is given by $K_{i}\equiv R_{ij}k_{j}$. Defining $\mu=\hat{K}\cdot\hat{q}$ as the angle between the transformed wave vector and Lagrangian separation q, we have the dot products

[TABLE]

The power spectrum in this frame is then simply given by substituting $K$ for $k$ in the real space expression; e.g. the term in the exponential can be written in terms of the vector K as

[TABLE]

In redshift space, unlike in real space, the Fourier product $\textbf{k}\cdot\textbf{q}$ requires azimuthal angle integrals due to the appearance of $\cos\phi$ in the final line of C.10; however, in Mii the azimuthal dependence is factored entirely into the Fourier factor $i\textbf{k}\cdot\textbf{q}$ such that the $\phi$ integral can be performed analytically:

[TABLE]

The remaining $q$ and $\mu$ integrals can then be calculated using the usual combination of spherical-Bessel decompositions and Hankel transforms with the help of the identity [50]:

[TABLE]

where $\rho=\sqrt{A^{2}+C^{2}}$ and the function $\tilde{G}_{m}(A,B,\rho)$ is given by

[TABLE]

${}_{2}F_{1}$ is the ordinary hypergeometric function, and we have defined

[TABLE]

In our specific case, $\rho=kq$, $B=-K^{2}Y/2$ and $A=kqc.$ Defining

[TABLE]

such that $I_{n}=I_{0}^{(n)}(A)$ is the nth full derivative of $I_{0}$ with respect to $A$, we have recursively

[TABLE]

The first two derivatives of $\tilde{G}^{(0)}$ are given by

[TABLE]

[TABLE]

These are sufficient to calculate all terms up to quadratic order in the bias expansion. In our fiducial cosmology at $z=0$ and along the line of sight ($\nu=1$) we find that the sums in $\tilde{G}_{m}$ converge at the sub-percent level within thirty summands.

From the above, contributions to the redshift-space Zeldovich power spectrum can be calculated using spherical Bessel transforms of the specific form

[TABLE]

where the scalar function $\mathcal{A}_{n}$ are tabulated in Table 2.

Appendix D Wiggle/No-Wiggle split

Most analyses of BAO data to date have employed empirical models for the post-reconstruction power spectrum or correlation function often motivated by theoretical calculations and calibrated to N-body simulations. Refs. [62, 26] showed that the analytic form of these empirical models can be interpreted within perturbation theory as a resummation of bulk displacements at the BAO scale. In this appendix we re-derive their results within our Zeldovich calculation, updating the scale dependences and redshift-space factors where appropriate.

Let us first examine the displaced-displaced cross spectrum in redshift space. Following refs. [62, 26] we split the displacement two-point function into $A^{dd}_{ij}=A^{dd,\rm nw}_{ij}+\Delta A^{dd,\rm w}_{ij}$, where the no-wiggle and wiggle pieces are calculated by substituting $P_{\rm nw}$ and $\Delta P_{\rm w}$ into Equation A.2. Making the assumption that the latter, $\Delta A^{dd,\rm w}_{ij}$, is small enough as to be perturbative131313Taking the nonlinear scale to be given by $k_{\rm nl}^{2}\Sigma^{2}(z)\sim 1$, we have $k_{\rm nl}\propto D^{-1}(z)$, such that the peak magnitude of $k_{i}k_{j}\Delta A^{\rm w}_{ij}$ is roughly in the few tenths of a percent range for our reference cosmology independent of redshift., we can Taylor expand the exponential in the Zeldovich integrand to get

[TABLE]

where we have used the transformed $K_{i}=R_{ij}k_{j}$ to encode redshift-space effects. Given that the no-wiggle spectrum reproduces the broadband scale dependence of the linear theory power spectrum, we can think of the no-wiggle exponential as resumming the non-BAO component of large scale bulk flows. Since the wiggle component contributes negligibly to the displacement power in the perturbative limit, keeping only one power of the wiggle power spectrum in our calculations serves to distinguish the effect of the IR bulk flows from BAO phenomena. The two-point functions entering the bias terms can likewise be split into no-wiggle and wiggle pieces, e.g. $U(q)=U^{\rm nw}+\Delta U^{\rm w}$, where again, roughly speaking, the former will contribute only to to the broadband power while the latter will give rise to oscillatory behavior. Keeping the above expression to order141414Note, however, that terms involving more powers of the wiggle displacement will be more suppressed than those involving no-wiggle displacements. $\mathcal{O}(k^{2}\Sigma^{2})$, and discarding terms that don’t contain any no-wiggle pieces, we then have

[TABLE]

where in the final line we have used the fact that the wiggle contributions will be confined in support around the BAO scale ($q\sim 100\,$Mpc) and the non-wiggle pieces vary smoothly at this Lagrangian separation, so we can pull the exponentiated no-wiggle contribution out of the integral as an average. Following ref. [62], we have defined the quantity $\bar{A}^{dd,\rm nw}_{ij}$ as the “average” of the un-barred quantity over the support of the wiggle component; to zeroth order in the approximation this is equivalent to evaluating $A_{ij}$ at the peak $q_{\rm max}$ of the support of the wiggle feature. Neglecting any angular effects in $\nu=\hat{q}\cdot\hat{k}$, which will enter at higher order in the wave number, we further have that $\bar{A}^{dd,\rm nw}_{ij}\simeq(X^{dd,\rm nw}+\frac{1}{3}Y^{dd,\rm nw})\delta_{ij}$151515The factor of a third, included also in ref. [62] but not in ref. [26], comes from the angular average $\langle\hat{q}_{i}\hat{q}_{j}\rangle=\delta_{ij}/3$. This can be justified by noting that the integral

$\frac{1}{2}\int\ d\mu\ e^{ikq\mu-\frac{1}{2}k^{2}\mu^{2}Y/2}=e^{-k^{2}Y/6}j_{0}(kq)+\mathcal{O}(k^{4}\Sigma^{4})$

(D.1)

We note, however, that this prescription is only approximate; for example, the same integral with an additional factor of $\mu$ in the integrand, relevant for the $b_{1}$ contribution, would instead yield $\exp[-3k^{2}Y/10]\,j_{1}(kq)$ at leading order. In general, bias contributions with more angular dependence will be damped more. This effect is automatically included in the full Zeldovich calculation.. Plugging in for the expression of $\textbf{K}=\textbf{R}^{\textbf{T}}\textbf{k}$, with $K^{2}=(1+f(f+2)\nu^{2})k^{2}$, the wiggle contribution to the power spectrum is then approximately

[TABLE]

where in the penultimate equality we have used the definition of the displaced field and defined $\Sigma_{dd}^{2}=(X^{dd,\rm nw}+\frac{1}{3}Y^{dd,\rm nw})(q_{max})$ to be evaluated at the peak of the wiggle displacements. This recovers the form of the empirical model in ref. [26] when we take the Eulerian bias to be $b_{1}^{E}=1+b_{1}$, and stick to the damping expansion approximation introduced in ref. [62]. Explicit expressions for $X$ and $Y$ are given in Equation 4.9. Taking $S\rightarrow 0$ in the above expression gives the unreconstructed power spectrum within this approximation.

We can now derive the analytical form of the reconstructed power spectrum for Rec-Sym in this approximation. Explicitly, we have

[TABLE]

The $\Sigma^{2}_{ab}$ are defined as in the $dd$ case. These are the same expressions as derived in ref. [26], though we differ on the expressions for the $\Sigma_{ab}^{2}$ that are involved. Our expressions also agree with those in refs. [19, 27] in the limit that $q_{max}\rightarrow\infty$, though we note that this limit doesn’t as accurately capture the damping of the feature since it resums the IR displacements at $q$ beyond the BAO scale. Adding the three spectra together, we recover the Kaiser limit as $k\rightarrow 0$, with different damping factors entering at different scales via $\Sigma^{2}_{ab}$. Note that in Rec-Sym the angular dependence of the damping is identical in each piece and is encoded within the $\nu$ dependence of $K^{2}$.

The reconstructed power spectrum with Rec-Iso requires a few additional modifications. The displaced-displaced auto spectrum is unchanged, and the shifted-shifted auto spectrum can be calculated by setting $f=0$ in all formulae, as noted in the main body of the text. However, the $ds$ cross spectrum requires more care, since the zero lag pieces do not transform equally. In particular, direct inspection of the exact expression in Equation 4.16 shows that we should instead define

[TABLE]

Note this expression differs in detail from that in ref. [26]. The cross spectrum is then instead

[TABLE]

where the angular dependence is subsumed into the defintion of $\Sigma_{ds,\rm iso}$. Unlike Rec-Sym, the damping factor in Rec-Iso is not captured by a single angular dependence.

We end this section with a discussion of the inclusion of higher bias terms and other corrections. As seen in the main body of the text, higher bias terms $b_{2}$ and $b_{s}$, incorporated in our Zeldovich calculation, contribute not only to the broadband but also serve to shift and smear the BAO feature itself. It might thus be of interest to extend the above approximation to include also these higher bias contributions. A potential avenue has been highlighted in ref. [26], although an approach closer to our perturbative bias expansion could also be explored.

Finally, the calculation in ref. [26] included a derivative bias, $b_{\nabla^{2}}\nabla^{2}\delta$, as a proxy to estimate the contributions of the higher bias operators. These derivative bias terms can easily be included in the above expressions by substituting $b_{1}\rightarrow b_{1}+k^{2}b_{\nabla^{2}}$. However, there is another context in which such a term might arise in which it would differ across the three pieces $dd$, $ds$ and $ss$: if the smoothing due to the $\Sigma^{2}_{ab}$’s as defined above do not accurately capture the IR bulk flows – for example if the broadband properities of $P_{\rm nw}$ are slightly off – but differ by some perturbatively small $k^{2}\delta\Sigma^{2}_{ab}$, the resulting correction could be corrected for by terms of the form $c_{ab}^{2}k^{2}P^{ab,\rm w}(k)$, where $c_{ab}^{2}$ would constants fit individually to $dd$, $ds$ and $ss$. Such corrections are essentially identical to the EFT corrections described in the text for the full Zeldovich calculation.

Appendix E Nonlinearities from the Lagrangian to Eulerian mapping

In standard density field reconstruction, each galaxy is shifted by a smoothed displacement field $\chi$ evaluated at the galaxy’s current Eulerian position $\textbf{x}=\textbf{q}+\boldsymbol{\Psi}(\textbf{q})$ (Equation 4.2). In the main body of the text, we worked in the approximation that $\chi(\textbf{x})\approx\chi(\textbf{q})$, with the understanding that nonlinear corrections would be suppressed by the smoothing scale $\sim\boldsymbol{\Psi}/R$. The goal of this appendix is to flesh out this statement by explicitly computing the leading order corrections to the reconstructed matter power spectrum due to the mapping nonlinearity in real space. For the sake of brevity we will defer the effects of other nonlinearities, such as those arising from dynamics or from translating between displacements and densities, to future work. Earlier treatments of such effects in Eulerian perturbation theory can be found in refs. [17, 21].

Assuming that the shift field $\chi(\textbf{q})$ defined in Lagrangian space is Gaussian, the displaced field with mapping nonlinearities unsuppressed is given by

[TABLE]

where we have kept the convention used in the main text to refer to the linear piece as $\boldsymbol{\Psi}^{d}=(1-\mathcal{S})\ast\boldsymbol{\Psi}$, referring to the nonlinear field as $\tilde{\boldsymbol{\Psi}}^{d}=\boldsymbol{\Psi}(\textbf{q})+\chi(\textbf{q}+\boldsymbol{\Psi})$. For the remainder of this appendix we will focus on corrections due to $\boldsymbol{\Psi}^{(d,2)}$.161616At one loop order all corrections due to $\boldsymbol{\Psi}^{(d,3)}$ are degenerate with the counterterms in our model. To see this, note that such corrections contractions with linear displacements, e.g.

$\langle\boldsymbol{\Psi}^{a}_{1,i}\boldsymbol{\Psi}^{(d,3)}_{2,j}\rangle=\langle\boldsymbol{\Psi}_{m}\partial_{m}\partial_{n}\chi_{j}\rangle\langle\boldsymbol{\Psi}^{a}_{1,i}\boldsymbol{\Psi}_{2,n}\rangle+\frac{1}{2}\langle\boldsymbol{\Psi}_{n}\boldsymbol{\Psi}_{m}\rangle\langle\boldsymbol{\Psi}^{a}_{1,i}\partial_{n}\partial_{m}\chi_{2,j}\rangle,\;a=d,s.$

(E.2)

Multiplied by the appropriate factor of $-\frac{1}{2}k_{i}k_{j}$, the two pieces on the right hand side Fourier transform into $\sim k^{2}P^{am}_{L}(k)$ and $\sim k^{2}P^{as}_{L}(k)$, respectively, thus taking the form of our counterterms $\sim k^{2}P^{ab}(k)$.

From the above,we can write the nonlinear displaced-displaced autospectrum as

[TABLE]

where we have defined

[TABLE]

as in the case of the nonlinear matter power spectrum (e.g. [38]). To calculate these we note that

[TABLE]

and

[TABLE]

where numerical subscripts refer to coordinates $\textbf{q}_{1,2}$ as usual.

The mapping corrections to the cross spectrum $P^{ds}$ can be similarly calculated. In this case we need the displacement correlators

[TABLE]

where all expectation values are evaluated at a point since $\boldsymbol{\Psi}^{s}$ receives no nonlinear corrections from the Eulerian-Lagrangian mapping and similarly

[TABLE]

to one loop order.

As might have been expected, the mapping corrections above all take the form of products of displacement two point functions and their derivatives. Roughly speaking these corrections each have amplitudes given by powers of the Zeldovich displacement $\Sigma^{2}$ and wavenumber $k$ capped at $R^{-1}$ by the smoothing filter; we can thus expect the corrections to enter at order $(\Sigma/R)^{4}$. For a $15h$ Mpc*-1* filter at $z=0$ this amounts to a percent-level effect, with smaller effects at higher $z$. Concretely, these two point functions can be calculated using

[TABLE]

where the ellipses denote all distinct permutations and the scalar functions are given by

[TABLE]

where we have used the identification $\chi=\boldsymbol{\Psi}^{s}.$ The remaining correlator $\left\langle\,\boldsymbol{\Psi}_{2,i}\boldsymbol{\Psi}_{1,j}\,\right\rangle$ is simply minus the non-zero lag piece of $A_{ij}$. Finally, when some or all of the displacement correlators in each product are contracted at the same point, as for example in the first term in Equation E.4 and Equation E.5, the resulting contribution becomes proportional to $k^{2}P_{ab,\rm Zel}$ and degenerate with the counterterms included in our model. The above corrections from the nonlinear Eulerian-Lagrangian mapping are plotted at $z=0$ for the smoothing scale $R=15h^{-1}$ Mpc in Figure 15. As expected, even at $z=0$ they are never more than a few percent of the total reconstructed power, though interestingly they can become comparable or larger than the Zeldovich $P^{dd}$ and $P^{ds}$ individually on scales where the Zeldovich spectra lose support. We caution that these curves do not include comparable corrections due to bias or dynamical nonlinearities.

We close with some general comments about nonlinearities in reconstruction. Firstly, the mapping corrections enumerated above are not the only ones at one-loop order; by focusing only on corrections due to $\boldsymbol{\Psi}^{(d,2)}$ we have explicitly avoided the $(13)$ contributions due to third-order mapping corrections. Moreover, as this was an exploratory exercise with which to evaluate the magnitude of mapping nonlinearities, we chose not to include the effects of bias, which would require the inclusion of terms such as $\langle\delta_{1}^{2}\boldsymbol{\Psi}^{(d,2)}_{i}\rangle$, though these will be in general decomposable into components much like those in Equation E.7. Finally, in addition to mapping nonlinearities, by approximating the shift vector $\chi$ with the smoothed Zeldovich displacement we have ignored nonlinearities induced by translating between the density field and displacements. We expect these will be of similar importance to the mapping corrections but defer their evaluation for future work, noting that only that both nonlinearities can be trivially reduced by pushing the smoothing scale $R$ deeper into the linear regime. That these effects are expansions in $\Sigma/R$ distinguishes them from nonlinear bias or dynamics.