Isobaric Reconstruction of the Baryonic Acoustic Oscillation

TL;DR

This paper introduces an isobaric reconstruction algorithm that effectively recovers the linear BAO signature from non-linear matter density fields, significantly improving the accuracy of cosmological distance measurements.

Contribution

The paper presents a novel iterative reconstruction method that accurately recovers the linear BAO signal without detailed dynamical modeling, enhancing measurement precision.

Findings

Reconstruction error within a few percent at pixel level.

Reduces fractional error in sound horizon scale measurement by a factor of 2.7.

Reconstructed field aligns well with linear theory, immune to bulk-flow smearing.

Abstract

In this paper, we report a significant recovery of the linear baryonic acoustic oscillation (BAO) signature by applying the isobaric reconstruction algorithm to the non-linear matter density field. Assuming only the longitudinal component of the displacement being cosmologically relevant, this algorithm iteratively solves the coordinate transform between the Lagrangian and Eulerian frames without requiring any specific knowledge of the dynamics. For dark matter field, it produces the non-linear displacement potential with very high fidelity. The reconstruction error at the pixel level is within a few percent, and is caused only by the emergence of the transverse component after the shell-crossing. As it circumvents the strongest non-linearity of the density evolution, the reconstructed field is well-described by linear theory and immune from the bulk-flow smearing of the BAO signature. Therefore this algorithm could significantly improve the measurement accuracy of the sound horizon scale. For a perfect large-scale structure survey at redshift zero without Poisson or instrumental noise, the fractional error is reduced by a factor of 2.7, very close to the ideal limit with linear power spectrum and Gaussian covariance matrix.

Equations

$$\det \left(\frac{\partial x_i}{\partial q_j}\right)=\det \left({\delta }_{ij}+{\partial }_i{\Psi }_j\right)=\frac{\bar{\rho }}{\rho }=\frac{1}{1+{\delta }_{\rho }}.$$

$$d_{\tau }\left[\rho \left(\mathbf{x}\left(\mathbf{q}\right)\right)\det \left(\frac{\partial x_i}{\partial q_j}\right)\left(\mathbf{q}\right)\right]=0$$

$$x^i({\xi }^{\mu })={\xi }^{\mu }{\delta }_{\mu }^i+d{x}^i={\xi }^{\mu }{\delta }_{\mu }^i+\frac{\partial (d\varphi )}{\partial {\xi }^{\nu }}{\delta }^{i\nu },$$

$${\partial }_{\mu }\left(\rho \sqrt{g}{e}_i^{\mu }{\delta }^{i\nu }{\partial }_{\nu }\left(d_{\tau }\varphi \right)\right)={\partial }_{\mu }\left(\rho \sqrt{g}{e}_i^{\mu }\left(d_{\tau }{x}^i\right)\right)=d_{\tau }\rho ,$$

$${\partial }_{\mu }\left(\rho \sqrt{g}{e}_i^{\mu }{\delta }^{i\nu }{\partial }_{\nu }(d\varphi )\right)=S(-\Delta \rho +C+E).$$

$${\varphi }^{rec}(\mathbf{\xi })=\sum _i^{iters}d{\varphi }^{(i)}\left({\mathbf{\xi }}^{(i)}\left(\mathbf{\xi }\right)\right).$$

$$P_{\mathrm{nl}}(k)$$

$$\mathrm{Wig}(k)=\frac{P_{\mathrm{nl}}(k)}{P_{\mathrm{sm}}(k)}-1\approx \exp (-k^2{\sigma }_v^2){\mathrm{Wig}}_{\mathrm{ini}}(k),$$

$$F_{\Delta \ln s}=\frac{\partial \ln P(k_i)}{\partial \ln s}{C}_{ij}^{-1}\frac{\partial \ln P(k_j)}{\partial \ln s}.$$

$$C_{ij}=\frac{1}{n_s-1}\sum _m^{n_s}\left(\frac{P_m({\mathbf{k}}_i)}{⟨P({\mathbf{k}}_i)⟩}-1\right)\left(\frac{P_m({\mathbf{k}}_j)}{⟨P({\mathbf{k}}_j)⟩}-1\right),$$

$$C_{ij}^g=\frac{2}{N_k}{\delta }_{ij}^K=\frac{4{\pi }^2}{k_i^2\Delta k_i{V}_s}{\delta }_{ij}^K$$

Full text

1. Introduction

Large-scale structure (LSS) surveys provide a wealth of cosmological information about our Universe. Particularly, the wiggling feature of the baryonic acoustic oscillation imprinted on the matter power spectrum encodes the characteristic sound horizon scale of $s\sim 100\rm Mpc/h$, which serves as a standard ruler for cosmic distance measurements (Eisenstein, 2003; Blake & Glazebrook, 2003; Hu & Haiman, 2003; Seo & Eisenstein, 2003). Since its first detection in SDSS by Eisenstein et al. (2005), BAO has become essential for almost all LSS surveys. Despite the fact that $s$ is relatively large, as the Universe evolves, the accuracy of measurements of $s$ from the power spectrum (or two-point correlation function) starts to decrease rather rapidly due to non-linear gravitational effects. Statistically, the information leaks into higher-order bispectra or trispectra etc., as the fields become more and more non-Gaussian with increased variance. This could be seen from the measurement of the so-called Fisher information curve, as a plateau starts to appear in the quasi-linear regime (Rimes & Hamilton, 2005, 2006). Meanwhile, the same structure formation process moves galaxies away from their initial locations and erases the BAO feature (Crocce & Scoccimarro, 2008). Eventually, one would lose the constraining power on cosmology due to the combined effect of non-Gaussianity and the degraded BAO signal (Seo & Eisenstein, 2007; Ngan et al., 2012).

This situation has led to many efforts to reconstruct the BAO signature and improve the measurement precision of the sound horizon scale $s$ from LSS. For example, a logarithmic transformation or Gaussianization processes seems to improve the non-Gaussianity and bring back some information that leaks into higher order statistics (Weinberg, 1992; Neyrinck, 2011). On the other hand, Eisenstein et al. (2007) first demonstrated that it is possible to sharpen the BAO peak in the correlation function and improve the measurement accuracy. This simple yet powerful algorithm (hereafter ZA reconstruction) has since been improved and applied to real measurements (Padmanabhan et al., 2012; Alam et al., 2016), producing stringent cosmological constraints(Mehta et al., 2012; Alam et al., 2016). Meanwhile, as the physical picture starts to clarify, more algorithms have started to emerge (Tassev & Zaldarriaga, 2012b; Schmittfull et al., 2015; Schmittfull et al., 2017).

Historically, the attempt to undo non-linear structure formation dates far back before these efforts. Peebles (1989, 1990) first applied the principle of least action to reconstruct the trajectories of Local Group. This technique was further developed by other authors (Croft & Gaztanaga, 1997; Brenier et al., 2003) but without any application to cosmology because of its computational cost. On the other hand, various approaches have tried to reverse the dynamical system, particularly concentrating on the velocity potential, without worrying about any decaying modes so that observational errors will not overpower the real initial fluctuations (Nusser & Dekel, 1992; Gramann, 1993; Narayanan & Croft, 1999; Monaco & Efstathiou, 1999). In a method to ensure the uniqueness of the reconstruction, Frisch et al. (2002) then solved the optimal mass transportation problem using Monge-Ampère-Kantorovich (MAK) method (Brenier et al., 2003; Mohayaee et al., 2006; Mohayaee & Sobolevskiĭ, 2008). Meanwhile, as computational processing power increases steadily, there are also recent attempts to directly sample all possible initial states using Markov Chain Monte Carlo method with fast forward modelling (Wang et al., 2009; Jasche & Wandelt, 2013; Wang et al., 2013; Kitaura, 2013). Unlike most previous works, Eisenstein et al. (2007) tried to recover the linear power spectrum instead of the initial density distribution, which is almost uniform, $\delta_{\rho}\ll 1$, and therefore practically much more challenging.

From the perspective of Lagrangian dynamics, the strong non-linearity present in the density field $\delta_{\rho}$ is attributed to both the gravitational evolution of the particle movement, manifested by the non-linear displacement $\bm{\Psi}(\mathbf{q})=\bm{x}-\mathbf{q}$ of particles, as well as the coordinate transform between the Lagrangian and Eulerian frames, where $\bm{x}$ and $\mathbf{q}$ are the Eulerian and Lagrangian positions respectively. With the assumption of the mass conservation and a uniform initial matter distribution, this coordinate transform could be expressed as

[TABLE]

Recent developments suggest that this coordinate transform is the major source of various theoretical and practical difficulties encountered in large-scale structure measurements (Tassev & Zaldarriaga, 2012a, b). Therefore, simply by avoiding this non-linear mapping, the displacement field itself would be much easier to handle theoretically. Of course, by definition, $\bm{\Psi}(\mathbf{q})$ is not an observable as it involves the unknown initial location $\mathbf{q}$. As a result, there have been efforts concentrating on the statistical extraction of $\bm{\Psi}(\mathbf{q})$, e.g. the maximum likelihood method by Tassev & Zaldarriaga (2012b).

However in principle, given a non-linear density map $\delta_{\rho}$, and assuming only the longitudinal component of the displacement field, i.e. $\bm{\Psi}=-\nabla\phi$, is cosmologically relevant, equation (1) reduces to the so-called elliptical Monge-Ampère equation and becomes solvable with appropriate boundary conditions111 The existence and uniqueness of the solution to the Monge-Ampère equation is out the scope of this paper. Theorems were proven by Bakelman (1957); Aleksandrov (1958) for the generalized solution with Dirichlet boundary conditions and Lions et al. (1986) for the Neumann problem. One could also see e.g. Trudinger & Wang (2008) for more detailed discussion. . Of course, this assumption only holds if the input map $\delta_{\rho}$ is generated from particle movement without any stream-crossing or rotation. In reality, these transverse components do emerge, but only at smaller scales (Chan, 2014). Therefore the induced reconstruction error in the BAO regime would be negligible.

Recently, Zhu et al. (2016) first proposed the potential isobaric reconstruction algorithm to reverse equation (1) and solve for the displacement potential $\phi$, and demonstrated that the resulting reconstruction of displacement potential is highly correlated with the initial density field. Pan et al. (2016) further showed that the Fisher information of the reconstructed field is significantly increased especially in the quasi-linear regime, and Yu et al. (2017) performed a similar study for the halo density field. Shortly after these papers were released, Schmittfull et al. (2017) also proposed a particle-based iterative reconstruction algorithm that performs similarly.

In this paper, we apply the same algorithm and particularly concentrate on cosmological constraints from measuring baryonic acoustic oscillations. In section 2., we first review this algorithm and demonstrate its ability to recover the nonlinear displacement field. We focus on the BAO signature in section 3., and show how the algorithm can improve the accuracy of the recovered wiggles and thereby improve cosmological measurements. Finally, we conclude and discuss in section 4. Throughout this paper, the N-body simulations were run with numerical code CUBE, with the box size of $1024~{}\rm Mpc/h$, and $512^{3}$ particles. We adopt the Planck 2015 results (Planck Collaboration et al., 2016) as our fiducial cosmological parameters.

2. The Potential Isobaric Reconstruction

As a non-linear partial differential equation, the numerical solution to equation (1) is nontrivial. Various numerical strategies, e.g. Zheligovsky et al. (2010); Liu et al. (2014), exist in the literature, but one caveat of directly applying these algorithms in cosmology is the frame mismatch between the Lagrangian gradient of the displacement $\partial_{i}\Psi_{j}=\partial\Psi_{i}/\partial q_{j}$ and the Eulerian density field $\rho(\bm{x})$ 222It has been demonstrated that with the assumption of no stream-crossing, there exists a potential field $\theta(\bm{x})$ that satisfies $\mathbf{q}=\bm{x}-\nabla_{x}\theta(\bm{x})$ (Arnold, 1979; Frisch et al., 2002), and therefore equation (1) is indeed well-defined and solvable with standard algorithms.. In this paper, we follow the multi-grid moving mesh algorithm proposed by Pen (1995, 1998) to iteratively solve the differential form of this mass conservation equation in a moving curvilinear coordinate system, which we name as the potential isobaric coordinate.

To present our method, we first take the differentiation of equation (1) with respect to some time-like iteration parameter $\tau$, which gives the conservation equation

[TABLE]

of a unit Lagrangian volume. For given Eulerian cell, it further simplifies to $d_{\tau}\rho=0$. Notice that both equations (1) and (2) are expressed locally, and therefore do not explain how the conservation would be maintained as long as the density and volume of the cell vary consistently. But from the perspective of the field, equation (2) is nothing but the continuity equation with respect to $\tau$, which is the equation we will solve. To trace the deformation of the fluid elements, we work in a moving numerical grid with curvilinear coordinate $\bm{\xi}=\{\xi^{\mu}\}=\{\xi^{1},\xi^{2},\xi^{3}\}$. To determine the physical position of the lattice, one needs to specify the coordinate transform relation $x^{i}(\xi^{\mu})$. Following Pen (1995, 1998), we borrow the notation from general relativity to describe a general curvilinear coordinate system. In this paper, we use Greek indices for curvilinear grid coordinates, and Latin indices for Cartesian coordinates. The coordinate transform is described by the triad $e^{i}_{\mu}=\partial x^{i}/\partial\xi^{\mu}$. The flat metric in Cartesian coordinates $x^{i}$ is simply the Kronecker delta $g_{ij}=\delta_{ij}$, while the curvilinear metric is given by $g_{\mu\nu}=e^{i}_{\mu}e^{j}_{\nu}\delta_{ij}$, both of which could be used to raise or lower the indices. Particularly we adopt the following ansatz at each iteration step

[TABLE]

where $\delta^{i}_{\mu}$ is the Kronecker delta function, and we use Einstein summation convention. Here $d\phi(\xi^{\mu})$ is the increment of the displacement potential; since we neglect the transverse component, this ansatz will help to minimize the mesh twisting.

As demonstrated by Pen (1998), in the isobaric frame, in which the mass per grid cell is kept constant, the continuity equation could be expressed as a linear elliptic equation

[TABLE]

where the dreibein $e^{\mu}_{i}$ is the inverse of the triad, or the Jacobian matrix, $e^{i}_{\mu}=\partial x^{i}/\partial\xi^{\mu}$, the volume element $\sqrt{g}=\det(e^{i}_{\mu})$. For given non-linear density map $\rho(\bm{x})$, we would like to reverse the frame deformation, and solve for the increment of potential $d\phi$, where $d\phi=d_{\tau}\phi d\tau=d_{\tau}\phi$ with the assumption that $d\tau=1$ since we are free to parametrize $\tau$. As shown in Zhu et al. (2016), this is done by substituting $d\rho$ in equation (4) with the negative non-linear density contrast $-\Delta\rho=\bar{\rho}-\sqrt{g}\rho$ together with a spatial smoothing $S$ to ensure the resultant $d\phi$ is small enough, i.e.

[TABLE]

Here $C$ and $E$ are the compression and expansion limiters for improving the numerical stability (Pen, 1998). Consequently, the triad $e^{i}_{\mu}$ is always kept positive definite. Then after moving the grid by equation (3), we estimate the density increment $d\rho$ from equation (4) and then repeat equation (5) with re-evaluated density contrast $\rho^{(i+1)}=\rho^{(i)}+d\rho^{(i)}$.

Eventually, the reconstructed displacement potential is simply a summation of all iterations,

[TABLE]

Here we have explicitly denoted the one-to-one mapping between the intermediate and final coordinate $\bm{\xi}^{(i)}\left(\bm{\xi}\right)$, which is only possible because of the positive definiteness of the transform matrix $e^{i}_{\mu}$ we achieved in solving equation (5). As a result, our final product $\phi^{rec}$ is then charted in a frame $\bm{\xi}$ that is almost identical to the theoretically defined Lagrangian coordinate $\mathbf{q}$, which means we have $\phi^{rec}(\bm{\xi})\approx\phi^{sim}_{long}(\mathbf{q})$ at the pixel level.

In the first column of Figure 1., i.e. panel $a1$ and $a2$ (with the same colorbar), we compare the simulation displacement potential $\phi^{sim}_{long}$ (panel $a2$) with the reconstructed potential $\phi^{rec}$ (panel $a1$). Here, $\phi^{sim}_{long}$ is calculated simply by taking the longitudinal component of the real displacement. The visual difference is not very obvious even when it is compared with the full displacement vector field, i.e. $\Psi^{sim}_{x}$ (panel $b1$) vs. $\Psi^{rec}_{x}=-\nabla_{x}\phi^{rec}$ (panel $b2$), when transverse contribution starts to emerge after the shell-crossing. As shown in the histogram of $\Delta\phi=\phi^{sim}_{long}(\mathbf{q})-\phi^{rec}(\bm{\xi})$ (panel $e1$), the reconstruction error is within a percent for $\phi^{rec}$, and is around ten-percent for the displacement $\bm{\Psi}$ itself (panel $e2$).

In panels $c1$ and $c2$, we also visually compare the input density map and its residual after $1400$ iterations. As shown, the non-vanishing residuals are basically around high density regions, where the displacement field $\bm{\Psi}$ would no longer be captured by longitudinal component $-\nabla\phi$ alone anymore. From the power spectra of the residual map in panel $d$, this starts to occur around the comoving Fourier scale $k\gtrsim 0.6h/\rm Mpc$ and peaks around $0.8h/\rm Mpc$. Even though the large scale power has been reduced to about six orders of magnitude smaller than the input, the small-scale power is still around unity.

3. Measuring the Baryonic Acoustic Oscillation

We are particularly interested in applying this isobaric reconstruction algorithm to detect the BAO signature in the power spectrum. As demonstrated by various authors, e.g. Crocce & Scoccimarro (2008), the sharp BAO features from the linear $P_{\rm ini}(k)$ are smeared out by the bulk flow of the galaxies. This could be seen by simply expressing the nonlinear density $\delta_{\rho}$ from the continuity equation (1) in Fourier space, and then taking the power spectrum (Crocce & Scoccimarro, 2006, 2008)

[TABLE]

Here $\mathbf{r}=\mathbf{q}-\mathbf{q}^{\prime}$, $\Delta\bm{\Psi}=\bm{\Psi}(\mathbf{q})-\bm{\Psi}(\mathbf{q}^{\prime})$, $G(k)$ is the so-called non-linear propagator $G^{2}(k)\approx\langle e^{i\mathbf{k}\cdot\Delta\bm{\Psi}}\rangle\approx\exp(-k^{2}\sigma^{2}_{v})$, $P_{\rm mc}(k)$ is the model-coupling power spectrum and $\sigma_{v}^{2}$ is the variance of the velocity. Since $P_{\rm mc}(k)$ dominates at large $k$ and is quite smooth, the BAO feature which could be defined as dividing $P_{{\rm nl}}(k)$ by some no-wiggle template $P_{\rm sm}(k)$, i.e.

[TABLE]

would be smeared out by the bulk flow $\Delta\bm{\Psi}$. Clearly, the reconstructed displacement field will no longer be affected by this smearing.

To examine the reconstructed BAO signature, we run a few N-body simulations with different transfer functions (with/without BAO, and with different sound horizon scale $s$) but the same initial random seeds. This enables us to reveal the signal with little numerical efforts but is free from the sample variance. As shown in the upper panel of Figure 2., at redshift $z=0$, BAO signature becomes basically invisible after the fourth peak (red-dashed line), whereas the non-linear displacement potential $\phi^{sim}_{long}$ (blue-dashed line) is almost identical to the initial density field (red-dotted line) in this $k$ range. The BAO of the reconstructed field $\phi^{rec}$ (black-solid line), are recovered with very high fidelity, although there seems to exist a percent-level change of broad band power. With an appropriate template, the sound horizon scale $s$ would be recovered at a level that would be very close to the linear theory. This is clear from the numerical derivative $d\ln P(k)/d\ln s$ in the lower panel of Figure 2, which is derived from the reconstruction of simulations with slightly different $s$. As a comparison, we also display the performance of ZA reconstruction method (Eisenstein et al., 2007) assuming a Gaussian smoothing length $15\rm Mpc/h$.

To evaluate the impact on the constraints of the cosmic distance, we estimate the Fisher information on the sound horizon scale $s$ for a noiseless survey. By definition, the Fisher matrix is

[TABLE]

Here $C_{ij}$ is the normalized covariance matrix

[TABLE]

which is estimated from $n_{s}=200$ N-body simulations. We then rescaled the large-scale component ($k<0.1$) of matrix $C_{ij}$ to coincide with a Gaussian covariance matrix $C^{g}_{ij}$

[TABLE]

assuming a survey volume of $V_{s}=1~{}({\rm Gpc}/h)^{3}$, where $N_{k}$ is the number of $k$ modes in such a survey.

In Figure 3., we show that the fractional error on the sound horizon $\Delta s/s$ measured from our reconstructed field is roughly $2.7$ times smaller than from the non-linear density field itself (red dashed line). This is very close to the accuracy one could ever achieve with the linear BAO (black-dotted line) feature and Gaussian covariance matrix. As a comparison, the same number is roughly $\sim 1.8$ for the ZA reconstruction.

4. Discussion

As the method of circumventing the strongest non-linearity in the structure formation process, we believe our reconstruction algorithm will significantly affect the way we measure large-scale structure in the future. However, in this paper, we did not take into account various complications, e.g. redshift space distortions (RSD) and biasing. In redshift space, the reconstructed field is simply a combination of the original displacement with the peculiar velocity, $\bm{\Psi}_{rsd}=\bm{\Psi}+\mathbf{v}/\mathcal{H}$. With some processing, e.g. a Wiener filter, one could then recover $\bm{\Psi}$. The bias, on the other hand, is probably more challenging. Without any preprocessing, one could perform the reconstruction for the biased density field itself, but the output $\phi^{rec}_{b}$ will then contain the bias information, and will be mapped on a grid other than $\mathbf{q}$. Alternatively, one could simply fix the bias parameters from the unconstructed density field, e.g. via HOD fitting or simply linear bias (Padmanabhan et al., 2012).

XW would like to thank Simon Foreman for proofreading the manuscript.