Bayesian weak lensing tomography: Reconstructing the 3D large-scale distribution of matter with a lognormal prior

TL;DR

This paper introduces a Bayesian method for 3D matter distribution reconstruction from weak lensing data, comparing Gaussian and lognormal priors, and demonstrating improved accuracy with the lognormal model in low-noise, small-scale scenarios.

Contribution

The work develops a Bayesian reconstruction algorithm using a lognormal prior for non-Gaussian density fields, outperforming Gaussian priors in certain realistic weak lensing scenarios.

Findings

Lognormal prior enforces non-negative densities.

Lognormal model better captures skewness and extremal values.

Higher pixel-wise correlation with true density in low-noise, small-scale cases.

Abstract

We present a Bayesian reconstruction algorithm that infers the three-dimensional large-scale matter distribution from the weak gravitational lensing effects measured in the image shapes of galaxies. The algorithm is designed to also work with non-Gaussian posterior distributions which arise, for example, from a non-Gaussian prior distribution. In this work, we use a lognormal prior and compare the reconstruction results to a Gaussian prior in a suite of increasingly realistic tests on mock data. We find that in cases of high noise levels (i.e. for low source galaxy densities and/or high shape measurement uncertainties), both normal and lognormal priors lead to reconstructions of comparable quality, but with the lognormal reconstruction being prone to mass-sheet degeneracy. In the low-noise regime and on small scales, the lognormal model produces better reconstructions than the normal model: The lognormal model 1) enforces non-negative densities, while negative densities are present when a normal prior is employed, 2) better traces the extremal values and the skewness of the true underlying distribution, and 3) yields a higher pixel-wise correlation between the reconstruction and the true density.

Tables

Table 1. Table 1: Quantitative comparison of reconstruction methods for tests A2 and A3. We compare the minimal and maximal values of the underlying and reconstructed fields ( δ 𝛿 \delta and δ ^ ^ 𝛿 \hat{\delta} , respectively).

redshifts	$\max \left(\hat{\delta }\right)$, $\max \left(\delta \right)$			$\min \left(\hat{\delta }\right)$, $\min \left(\delta \right)$
	lognormal	WF	N-body	lognormal	WF	N-body
1.	4.22	2.95	4.92	-0.63	-0.91	-0.79
0.25	7.38	5.00	8.49	-0.77	-1.28	-0.87
0.0	8.32	5.55	10.26	-0.99	-1.57	-0.89
light cone	32.90	20.81	70.20	-0.98	-4.28	-0.92

Table 2. Table 2: Quantitative comparison of reconstruction methods for tests A2 and A3. We compare the lognormal to the Gaussian prior in terms of the mean square pixel-wise difference (first column) and the Pearson-correlation coefficient (second column) between the reconstruction and the underlying field (which we denote δ ^ ^ 𝛿 \hat{\delta} and δ 𝛿 \delta , respectively).

redshifts	$⟨{\left(\hat{\delta }-\delta \right)}^2⟩/{\sigma }_{\delta }^2$		$\sqrt{⟨\hat{\delta }\delta ⟩/\left({\sigma }_{\delta }{\sigma }_{\hat{\delta }}\right)}$
	lognormal	WF	lognormal	WF
1.	0.66	0.69	0.58	0.56
0.25	0.73	0.77	0.52	0.49
0.0	0.73	0.77	0.53	0.48
light cone	0.65	0.73	0.61	0.51

Equations

$${\mathcal{C}}_{{\mathbf{x}}^{\prime }(\mathbf{x},\tau )}={\delta }_{\mathrm{D}}\left({\mathbf{x}}^{\prime }-\mathbf{x}\right){\delta }_{\mathrm{D}}\left(\tau -|\mathbf{x}|/c\right).$$

$$\rho (\mathbf{x},\tau )=\bar{\rho }\left[1+\delta (\mathbf{x},\tau )\right]=\bar{\rho }\left[1+{\int }_{\mathbf{y}}D(|\mathbf{x}-\mathbf{y}|,\tau ,{\tau }_0)\phi (\mathbf{y},{\tau }_0)\right]=\bar{\rho }\left[1+D_{\mathbf{xy}}(\tau ,{\tau }_0){\phi }_{\mathbf{y}}({\tau }_0)\right].$$

$$\phi \hookleftarrow \mathcal{N}(\phi |\bar{\phi },\Phi ).$$

$$\rho (\mathbf{x},\tau )=\bar{\rho }\left[1+\delta (\mathbf{x},\tau )\right]=\bar{\rho }\exp \left[D_{\mathbf{xy}}(\tau ,{\tau }_0){\phi }_{\mathbf{y}}({\tau }_0)\right].$$

$${\bar{\phi }}_{\mathbf{x}}(\tau )=-\frac{1}{2}{D}_{\mathbf{xy}}^{-1}(\tau ,{\tau }_0){\left[D(\tau ,{\tau }_0)\Phi ({\tau }_0)D(\tau ,{\tau }_0)\right]}_{\mathbf{yy}},$$

$${\bar{\phi }}_{\mathbf{x}}(\tau )=-\frac{1}{2}D(\tau ,{\tau }_0){\hat{\Phi }}_{\mathbf{x}}({\tau }_0),$$

$$\delta (\mathbf{k},\tau )=\sum _{n=1}{D}^{(n)}(\tau ){\delta }^{(n)}(\mathbf{k}),$$

$$D_{\mathbf{xy}}(\tau )$$

$${\delta }_{{\mathbf{x}}^{\prime }}={\mathcal{C}}_{{\mathbf{x}}^{\prime }(\mathbf{x},\tau )}\exp \left[D_{\mathbf{xy}}(\tau ){\phi }_{\mathbf{y}}(\tau )\right]-1.$$

$${\delta }_{{\mathbf{x}}^{\prime }}=\exp \left[D_{{\mathbf{x}}^{\prime }{\mathbf{y}}^{\prime }}{\phi }_{{\mathbf{y}}^{\prime }}\right]-1.$$

$${\mathbf{ϵ}}^{\mathrm{s}}\hookleftarrow \mathcal{N}({\mathbf{ϵ}}^{\mathrm{s}}|0,N^{\mathrm{s}}),N_{ij}^{\mathrm{s}}={\delta }_{ij}{\sigma }_{ϵ}^2.$$

$$ϵ=g+{ϵ}^{\mathrm{s}},$$

$$g=\frac{\gamma }{1-\kappa }\approx \gamma .$$

$${\gamma }_1(\mathbf{\theta })=\frac{1}{2}\left({\partial }_1^2-{\partial }_2^2\right)\psi (\mathbf{\theta });{\gamma }_2(\mathbf{\theta })={\partial }_1{\partial }_2\psi (\mathbf{\theta });\kappa (\mathbf{\theta })=\frac{1}{2}({\partial }_1^2+{\partial }_2^2)\psi (\mathbf{\theta }).$$

$$\psi (\mathbf{\theta })=\frac{2}{{\mathrm{c}}^2}\underset{0}{\overset{r^i}{\int }}\mathrm{d}r\frac{r^i-r}{r{r}^i}\varphi (r,r{\theta }_1,r{\theta }_2),$$

$$\partial_{k}\partial_{l}\psi(\boldsymbol{\theta})=\frac{2}{\mathrm{c}^{2}}\int\limits_{0}^{r^{i}}\mathrm{d}r\,W(r;r^{i})\partial_{r_{k}}\partial_{r_{l}}\phi\bigl{(}r,r\theta_{1},r\theta_{2}\bigr{)},$$

$$W(r;r^i)=\frac{r(r^i-r)}{r^i}.$$

$$T_{ij}({\mathbf{x}}^{\prime })={\partial }_{x_i^{\prime }}{\partial }_{x_j^{\prime }}\varphi ({\mathbf{x}}^{\prime }),i,j\in (0,1,2),$$

$${\nabla }^2\varphi ({\mathbf{x}}^{\prime })=\frac{3}{2}{\Omega }_{\mathrm{m}}{H}_0^2\frac{\delta ({\mathbf{x}}^{\prime })}{a(|{\mathbf{x}}^{\prime }|/\mathrm{c})},$$

$$\exp \left(D_{{\mathbf{x}}^{\prime }{\mathbf{y}}^{\prime }}{\phi }_{{\mathbf{y}}^{\prime }}\right)-1=D_{{\mathbf{x}}^{\prime }{\mathbf{y}}^{\prime }}{\phi }_{{\mathbf{y}}^{\prime }}+\mathcal{O}({\phi }^2)\approx a(|{\mathbf{x}}^{\prime }|){\delta }_{\mathrm{D}}({\mathbf{x}}^{\prime }-{\mathbf{y}}^{\prime }){\phi }_{{\mathbf{y}}^{\prime }}+\mathcal{O}({\phi }^2),$$

$$d_i={\mathbf{ϵ}}_i=(R(\phi ){)}_i+{\mathbf{ϵ}}_i^{\mathrm{s}}+{\mathbf{n}}_i={\mathcal{G}}_{ij}{\mathcal{I}}_{j{\mathbf{r}}^{\prime }}{\mathcal{T}}_{{\mathbf{r}}^{\prime }{\mathbf{z}}^{\prime }}{\mathcal{P}}_{{\mathbf{z}}^{\prime }{\mathbf{y}}^{\prime }}\left[\exp \left(D_{{\mathbf{y}}^{\prime }{\mathbf{x}}^{\prime }}{\phi }_{{\mathbf{x}}^{\prime }}\right)-1_{{\mathbf{y}}^{\prime }}\right]+{\mathbf{ϵ}}_i^{\mathrm{s}}+{\mathbf{n}}_i.$$

$$P(\phi |d)=\frac{P(d|\phi )P(\phi )}{P(d)}.$$

$$P(d|\phi )=\int \mathcal{D}{\mathbf{ϵ}}^{\mathrm{s}}P\left(d|\phi ,{\mathbf{ϵ}}^{\mathrm{s}}\right)P\left({\mathbf{ϵ}}^{\mathrm{s}}\right)=\mathcal{N}\left[d|R(\phi ),N\right],$$

$$-\ln P(\phi |d)\widehat{=}\frac{1}{2}{\left[d-R(\phi )\right]}^{†}{N}^{-1}\left[d-R(\phi )\right]+\frac{1}{2}{\left(\phi -\bar{\phi }\right)}^{†}{\Phi }^{-1}\left(\phi -\bar{\phi }\right),$$

$$-\frac{\delta \ln P(\phi |d)}{\delta {\phi }_{\mathbf{u}}}={\Phi }_{\mathbf{uq}}^{-1}{\left(\phi -\bar{\phi }\right)}_{\mathbf{q}}+{\left[d-R(\phi )\right]}_i{N}_{ij}^{-1}{\mathcal{G}}_{jk}{\mathcal{I}}_{k{\mathbf{r}}^{\prime }}{\mathcal{T}}_{{\mathbf{r}}^{\prime }{\mathbf{z}}^{\prime }}{\mathcal{P}}_{{\mathbf{z}}^{\prime }{\mathbf{y}}^{\prime }}\left[\exp (D_{{\mathbf{y}}^{\prime }{\mathbf{x}}^{\prime }}{\phi }_{{\mathbf{x}}^{\prime }})\ast D_{{\mathbf{y}}^{\prime }\mathbf{u}}\right],$$

$$P(d|\phi )=\int \mathcal{D}{\mathbf{ϵ}}^{\mathrm{s}}\int \mathcal{D}{z}^{s}P(d|\phi ,{\mathbf{ϵ}}^{\mathrm{s}},z^s)P({\mathbf{ϵ}}^{\mathrm{s}})P(z^s|z_{\mathrm{photo}})=\int \mathcal{D}{z}^s\mathcal{N}[d|R^{z^s}(\phi ),N]P(z^s|z_{\mathrm{photo}}).$$

$$⟨d{⟩}_{P(d|\phi )}$$

$$⟨d_i{⟩}_{P(d|\phi )}$$

$$\equiv \left[\mathcal{G}\bar{\mathcal{I}}\bar{R}\right](\phi {)}_i\equiv \tilde{R}(\phi {)}_i,$$

$$⟨d{d}^{†}{⟩}_{P(d|\phi )}=N+⟨(R(\phi ))(R(\phi ){)}^{†}{⟩}_{P(z|z_{\mathrm{ph}})}.$$

$$⟨\left((R(\phi ){)}_i(R(\phi ){)}_j\right){⟩}_{P(z|z_{\mathrm{ph}})}=⟨{(R(\phi ))}_i{⟩}_{P(z|z_{\mathrm{ph}})}⟨{(R(\phi ))}_j{⟩}_{P(z|z_{\mathrm{ph}})}=(\tilde{R}\phi {)}_i(\tilde{R}\phi {)}_j$$

Full text

I Introduction

Weak gravitational lensing of galaxies offers a unique way to study the distribution of matter in the Universe (see  (Bartelmann and Schneider, 2001) for a review on weak gravitational lensing and (Kilbinger, 2015; Hoekstra and Jain, 2008; Refregier, 2003) for recent reviews on weak galaxy lensing). Lensing by structures along the line of sight causes distortions in the images of distant galaxies (which in this context are often referred to as sources), which leads to correlations between the apparent shapes of these galaxies. The dominant and most easily detectable image distortion that lensing induces is a shearing of the galaxy images. Because of this, the effect is often referred to as cosmic shear.

Galaxy shape measurements allow to constrain the clustering of matter at different scales and redshifts, which can then be translated into constraints on cosmological models and their parameters. The integrated lensing signal is mostly sensitive to a combination of the cosmic mean matter density $\Omega_{\mathrm{m}}$ and the matter power spectrum amplitude $\sigma_{8}$ (Bernardeau et al., 1997). Tomographic methods can yield additional constraints on all parameters guiding the growth rate of structure (Hu, 1999; Simon et al., 2004), notably also on the properties of dark energy (Hu, 2002; Huterer, 2002; Song and Knox, 2004), but rely on low photometric redshift uncertainties (Ma et al., 2006). Since lensing is a direct probe of the total matter, luminous and dark, it can be combined with measurements of the luminous matter distribution in order to learn about the relationship between baryons and dark matter. Another important feature of weak galaxy lensing lies in its ability to probe the matter distribution over a wide range of scales, from many tens of Mega parsec, where structure formation is still linear today and comparably easy to model, down to non-linear sub-Mega parsec scales. Due to its sensitivity to such a wide range of scales, lensing can provide a large amount of information to constrain models of non-linear structure formation and cosmology, in particular if also higher-order statistics are considered.

First firm statistical detections of cosmic shear were reported in 2000 by four different groups (Bacon et al., 2000; Van Waerbeke et al., 2000; Kaiser et al., 2000; Wittman et al., 2000). Since then, the field has seen a tremendous increase in the amount and quality of lensing data (Scoville et al., 2007; Heymans et al., 2012; de Jong et al., 2013; Tyson et al., 2001; Jarvis et al., 2016, e.g.) as well as a notable improvement in analysis techniques (e.g. Miller et al., 2007; Kitching et al., 2008; Miller et al., 2013; Zuntz et al., 2013; Gentile et al., 2012; Schneider et al., 2015; Benjamin et al., 2013; Heymans and Heavens, 2003; King and Schneider, 2003; Hildebrandt et al., 2012). Estimates of cosmological parameters have been inferred from these data by comparing the power spectrum of the fully projected 2D shear field (or related quantities) as well as auto- and cross-spectra in a number of redshift bins to theoretical predictions (Massey et al., 2007; Jee et al., 2013; Benjamin et al., 2007; Kilbinger et al., 2013; Schrabback et al., 2010; Benjamin et al., 2013; Hildebrandt et al., 2017; Jee et al., 2016; Heymans et al., 2013; Abbott et al., 2016). Several authors have also investigated the additional constraining power that can be achieved by incorporating third-order statistics and/or shear peak counts and correlations (Bernardeau et al., 1997; Takada and Jain, 2004; Semboloni et al., 2011; Dietrich and Hartlap, 2010; Marian et al., 2011, 2012, 2013; Hilbert et al., 2012; Liu et al., 2015; Kacprzak et al., 2016), which helps to break parameter degeneracies.

The full information content, however, lies in the three-dimensional non-linear shear field. 3D weak shear analysis methods have been proposed by a number of authors (Heavens, 2003; Castro et al., 2005; Kitching et al., 2011) and were recently applied to data from the Canada France Hawaii Telescope Lensing Survey (CFHTLenS)  (Kitching et al., 2014). Furthermore, the measured shear field can be used for 3D reconstructions of the underlying density field. This can then be directly compared to models of structure formation and simplifies the cross correlation with other tracers of matter. Algorithms that invert the lens equation to obtain the density have been worked out by a number of authors (Kaiser and Squires, 1993; Hu and Keeton, 2002; Taylor, 2001). Since this inversion is under-constrained, it requires some regularization method or choice of prior on the density field. Most of the algorithms employ a Wiener filter, which corresponds to a normal (Gaussian) prior, that can be complemented with information about galaxy clustering (Simon, 2013; Szepietowski et al., 2014).

We aim to extend the work on tomographic reconstruction of the 3D matter distribution by designing a fully Bayesian reconstruction algorithm that uses a lognormal prior on the density field. The algorithm is designed to reconstruct the 3D cosmic density fluctuation field $\delta(\mathbf{x},\tau)$ from weak galaxy lensing data, i.e. a measurement of galaxy ellipticities at different (photometrically measured) redshifts. Its derivation is based on the language of information field theory (Enßlin et al., 2009), which has already been used to address similar tomography problems (Greiner et al., 2016). We do not make use of the flat-sky approximation i.e. lines of sight are allowed to be non-perpendicular to a fixed 3D coordinate grid. Further, we do not bin the data into pixels but take each galaxy into account as an individual contribution to the likelihood. This allows us to incorporate distance uncertainties of individual galaxies instead of sample redshift distributions.

In contrast to a normal prior for the density field, as it has often been assumed before, a lognormal prior automatically enforces the strict positivity of the field and allows to capture some of the non-Gaussian features that are imprinted on the density distribution by non-linear structure formation. Hubble was the first to notice that galaxy number counts could be well approximated by a lognormal distribution (Hubble, 1934). Characterizing the matter overdensities in the Universe as a lognormal field was first assessed by Coles & Jones in 1991 (Coles and Jones, 1991). Subsequent studies showed that a logarithmic mapping of the nonlinear matter distribution can partly re-gaussianize the field, and that non-linear features in the matter power spectrum can be reproduced by a lognormal transformation of the linear matter power spectrum (Neyrinck et al., 2009; Greiner and Enßlin, 2015). A lognormal prior has already been used and shown to be superior to a Gaussian one in Bayesian algorithms that reconstruct the large-scale matter distribution from the observed galaxy distribution (Enßlin et al., 2009; Kitaura et al., 2010; Jasche et al., 2010). Lognormal distributions have also already been considered in the context of weak lensing: Analyses of ray-tracing simulations and the Dark Energy Survey (DES) Science Verification data showed that the 1-point distribution function of the lensing convergence is better described by a lognormal than a Gaussian model (Taruya et al., 2002; Clerkin et al., 2016). Also the cosmic shear covariance can be modeled to better accuracy under the assumption that the underlying convergence field follows a lognormal distribution instead of a Gaussian one (Hilbert et al., 2011).

Bayesian inference methods are widely used in weak shear analyses, most prominently in the context of shear measurements from galaxy images (Miller et al., 2007; Bernstein and Armstrong, 2014; Schneider et al., 2015). Recently, notable effort has been put into developing a fully Bayesian analysis pipeline that propagates all uncertainties consistently from the raw image to the inferred cosmological parameters (Alsing et al., 2016; Heavens et al., 2016).

This paper is organized as follows: This introduction is followed by a short section, Sec. II, in which we briefly introduce the notations and coordinate systems that will be used in the derivation of the formalism. Our lognormal prior model for the density is described in detail in Sec. III. In Sec. IV, we present the data model, i.e. the lensing formalism that connects the data from a cosmic shear measurement to the underlying density field and give a brief overview over its implementation in Sec. V. The maximum a posteriori estimator that is used to infer the matter distribution is introduced in Sec. VI and extended to include redshift uncertainties of individual sources in Sec. VII. In Sec. VIII, we show results of the density reconstruction on increasingly realistic mock data. We conclude this work with a summary and discussion in Sec. IX.

II Coordinate systems and notational conventions

In the derivation of the formalism we work with three different types of coordinate systems. First, we use three-dimensional purely spatial comoving coordinates $\mathbf{x}=(x_{0},x_{1},x_{2})$ at fixed comoving lookback time $\tau$, that, combined with the time coordinate, form the four-dimensional coordinate system $(\mathbf{x},\tau)$.

Second, we use a 3D comoving coordinate system on the light cone of an observer at the origin. Vectors on the light cone are marked by a prime, e.g. $\mathbf{x}^{\prime}$. Since $\mathbf{x}^{\prime}$-coordinates implicitly define a comoving lookback time $\tau(\mathbf{x}^{\prime})=|\mathbf{x}^{\prime}|/\mathrm{c}$ (where $\mathrm{c}$ denotes the speed of light), we omit spelling out the time explicitly and write $A(\mathbf{x}^{\prime})=A(\mathbf{x}^{\prime},\tau)$ for any quantity $A$ that is defined on the light cone. The operation that links quantities on the light cone to their corresponding quantities in 4D spacetime can be encoded in a projection operator with kernel

[TABLE]

Third, we employ a set of coordinate systems on the light cone, in which each system is orientated such that one axis points into the direction of a source galaxy. These line of sight coordinate systems are centered on the observer and spanned by the vectors $(\hat{\mathbf{r}}_{0}^{i},\hat{\mathbf{r}}_{1}^{i},\hat{\mathbf{r}}_{2}^{i})$, where $\hat{\mathbf{r}}_{0}^{i}(\mathbf{x}^{\prime})$ points into the direction of the $i$th source galaxy and the normal vectors $(\hat{\mathbf{r}}_{1}^{i},\hat{\mathbf{r}}_{2}^{i})$ span the two-dimensional plane perpendicular to $\hat{\mathbf{r}}_{0}^{i}$. The radial comoving distance of each galaxy from the observer is denoted $r^{i}=|\mathbf{r}^{i}|\equiv|\mathbf{r}_{0}^{i}|$. The transformation from the light cone system into the line of sight (LOS) system of a source galaxy $i$ is achieved by a rotation of the $x^{\prime}_{0}$-axis into the $i$th line of sight, while $x^{\prime}_{1}$ and $x^{\prime}_{2}$ get aligned with the image coordinates of the galaxy observation. We label the corresponding transformation operators ${\mathcal{R}^{i}}_{\mathbf{r}\mathbf{x}^{\prime}}$.

III Prior model

The comoving density $\rho(\mathbf{x},\tau)$ can be split into its time-independent spatial mean $\bar{\rho}$, and a perturbation $\delta\rho(\mathbf{x},\tau)=\bar{\rho}\,\delta(\mathbf{x},\tau)$. The fractional overdensity $\delta(\mathbf{x},\tau)$ is commonly modeled as a homogeneous, isotropic Gaussian field with zero mean and power spectrum $P_{\delta}(k,\tau)$. This is an excellent description at early times where fluctuations are very small, as e.g., shown by observations of the cosmic microwave background radiation (CMB). At linear level, valid for $\delta\ll 1$, the time evolution of the density field can be described by

[TABLE]

In the last expression we have introduced a short-hand notation that will be used in the rest of the paper: repeated indices are integrated over if they do not appear on both sides of the equation. In Eq. (2), $D(|\mathbf{x}-\mathbf{y}|,\tau,\tau_{0})$ is the integration kernel of a linear homogeneous and isotropic, but possibly scale-dependent, growth operator. The field $\varphi$ is an isotropic and homogeneous random field whose values are drawn from a multivariate normal distribution with mean $\bar{\varphi}$ and covariance $\Phi$,

[TABLE]

Here, it describes the three dimensional density fluctuations at time $\tau_{0}$ and is translated to other times $\tau$ by the growth operation $D_{\mathbf{x}\mathbf{y}}(\tau,\tau_{0})$. This implies $\bar{\varphi}=0$ for this linear Gaussian model.

The description Eq. (2) breaks down when $\delta\not\ll 1$ such that non-linearities become important. A possible way to account for non-linearities is to include higher-order terms from a perturbation series expansion of the full non-linear evolution equations. However, a further shortcoming of the model (2) is that it allows arbitrarily negative density contrasts, which physically can not be smaller than -1. To obtain a strictly positive density field, we instead modify Eq. (2) by an exponential:

[TABLE]

Since the expectation value of the density $\rho(\mathbf{x},\tau)$ must equal $\bar{\rho}$, the mean of $\varphi$ must be set to

[TABLE]

for every time $\tau$. Note that we integrate over the index $\mathbf{y}$ in this expression, i.e. the diagonal of the composite operator in square brackets is treated as a field.

For a local growth operator, $D_{\mathbf{x}\mathbf{y}}(\tau,\tau_{0})=D(\tau,\tau_{0})\delta_{\mathrm{D}}(\mathbf{x}-\mathbf{y})$, this mean correction simplifies to

[TABLE]

where we defined $\hat{\Phi}_{\mathbf{x}}(\tau_{0})\equiv\Phi_{\mathbf{x}\mathbf{x}}(\tau_{0})$.

The Gaussian field $\varphi$ and the growth operator $D$ can be related to known quantities. To see this, consider the expansion of the Fourier modes of $\delta$ in Eulerian perturbation theory (see e.g. Bernardeau et al. (2002)),

[TABLE]

where $\delta^{(n)}$ are convolutions of $n$ initial fields $\delta(\mathbf{k},\tau_{0})$ with an integration kernel that changes from order to order. The first term in this series is $D^{(1)}(\tau)/D^{(1)}(\tau_{0})\delta_{0}(\mathbf{k},\tau_{0})$, where $D^{(1)}(\tau)$ is the growing solution to the linearized growth equation (Linder and Jenkins, 2003).

We use this analogy to motivate the simplest possible form of the growth operation in the lognormal model and write

[TABLE]

where we have set $D^{(1)}(\tau_{0})=1$. This approximation erases any a-priori assumption of scale-dependent growth and mode-coupling of the log field $\varphi$. Such a simplification is viable since the model in Eq. (8) describes only our prior assumptions about the density field $\rho$. The algorithm will find the most probable realization of $\varphi$ for a fixed growth operator $D$ given the data. If a scale-dependence is favored by the data, it will be recovered, at least partially, in the estimate of $\varphi$.

Our algorithm also allows to incorporate a more general growth operation at the expense of computation time and memory usage. The application of the most general $D_{\mathbf{x}\mathbf{y}}(\tau)$ generates a four-dimensional field: three-dimensional spatial comoving volumes for every time-slice $\tau$. This very large volume is then restricted to a three-dimensional cut by application of the light cone operator [Eq. (1)]. The prior model for the matter overdensity on the light cone then becomes

[TABLE]

For the purpose of comparing the lognormal to the Gaussian density model we choose the simplest form of $D$, Eq. (8), since it can be applied to the three-dimensional field on the light cone

[TABLE]

More complicated models, e.g. motivated by perturbation theories for large-scale structure, could be envisioned but lie beyond the scope of this work. A growth model that is better informed about the spatial correlations in the density field could potentially improve the reconstruction in regions with low signal to noise (see also the discussion in Sec. IX.

IV Data model

In this section we establish the analytic relation between the signal field $\varphi$, the field of overdensities $\delta$ that we aim to reconstruct, and the data that is obtained from a weak lensing measurement.

Weak galaxy lensing surveys produce galaxy image ellipticities that can be quantified, e.g., by a complex number $\epsilon=\epsilon_{1}+i\epsilon_{2}$ such that $\epsilon=(a-b)/(a+b)\exp(2\mathrm{i}\eta)$ for an ellipse with major axis $a$, minor axis $b$, and position angle $\eta$ (e.g. Blandford et al., 1991). We use the common approximation that the components of the intrinsic source galaxy ellipticity $\boldsymbol{\epsilon}^{\mathrm{s}}=(\epsilon^{\mathrm{s}}_{1},\epsilon^{\mathrm{s}}_{2})$, which define the shape that would be observed in the absence of lensing111We use vector notation, e.g. $\boldsymbol{\epsilon}$, to denote the tuple of real and imaginary part $(\epsilon_{1},\epsilon_{2})$ of a complex number $\epsilon{=}{\epsilon_{1}}+i{\epsilon}_{2}$., follow a global bivariate Gaussian distribution with zero mean and variance $\sigma_{\epsilon}^{2}$ per component:

[TABLE]

This approximation has shortcomings (see e.g. (Miller et al., 2013)), but serves for the proof of our concept, since we create the mock data on which we test the algorithm with exactly this shape noise model. In the future, more elaborated (Bayesian hierarchical) shear estimators, that e.g. take into account galaxy properties, can be incorporated into the algorithm (Miller et al., 2007; Schneider et al., 2015; Bernstein and Armstrong, 2014).

Lensing distorts the galaxy images in shape and size (Gunn, 1967; Kochanek, 1990; Blandford et al., 1991; Seitz et al., 1994; Kaiser, 1995). If the distortion is small, i.e. in the limit of weak lensing, the relation between intrinsic source ellipticity and observed ellipticity can be linearized and simplifies to (Seitz and Schneider, 1997; Krause and Hirata, 2010)

[TABLE]

where $g$ is the reduced shear. The reduced shear combines the effect of anisotropic lensing distortions, encoded in the shear $\gamma=\gamma_{1}+i\gamma_{2}$, and the isotropic distortion, encoded in the convergence $\kappa$

[TABLE]

If $\kappa\ll 1$, which is often the case for galaxy lensing, the reduced shear can be approximated by the shear itself $g\approx\gamma$.

The shear and convergence at angular position $\boldsymbol{\theta}$ are related to the lensing potential $\psi$ by

[TABLE]

The lensing potential $\psi$ is an integrated measure of scalar perturbations to the background metric along the perturbed photon geodesic. Integrating the perturbations along the unperturbed geodesic turns out to be an excellent approximation Shapiro and Cooray (2006); Krause and Hirata (2010). Working in this so-called Born approximation, choosing the Newtonian gauge (Mukhanov et al., 1992) and assuming no anisotropic stress, the lensing potential can be written as a weighted projection of the peculiar Newtonian gravitational potential $\phi$ along the line of sight. For a source at LOS distance $r^{i}$, this integration reads

[TABLE]

where we have assumed a spatially flat Universe. Applying the angular derivatives in Eq. (14) to the expression for the lensing potential in Eq. (15), we get

[TABLE]

where $k,l\in(1,2)$, and the lensing efficiency

[TABLE]

In practice the distance to the source $r^{i}$ cannot be determined directly but follows from the photometrically measured redshift $z^{i}$. Photometrically measured redshifts are associated with a relatively high error, $\sigma_{z}/(1+z)\approx 0.04-0.06$ (Hildebrandt et al., 2012). In its most simple form the algorithm ignores this uncertainty. We will use this simplified model to validate the functionality of the algorithm in terms of reconstructing non-linear structures in the lognormal approximation. Redshift uncertainties will be included later in Section VII.

The lensing shear is completely determined by the second derivatives of the lensing potential perpendicular to the LOS. The tidal tensor $\partial_{r_{k}}\partial_{r_{l}}\phi(\mathbf{r})$ along the LOS of the $i$th source galaxy is obtained by rotating the tidal tensor on the global coordinate grid $\mathbf{x}^{\prime}$

[TABLE]

into the specific coordinate system (with coordinates $\mathbf{r}^{i}$) that points into the direction of this $i$th galaxy and projecting it onto the $(r^{i}_{1}-r^{i}_{2})$-plane perpendicular to the LOS.

The last relation required to connect the data to the density fluctuations is Poisson’s equation. It relates the potential $\phi(\mathbf{x}^{\prime})$ to the density fluctuations $\delta(\mathbf{x}^{\prime})$,

[TABLE]

where $H_{0}$ denotes the Hubble constant (which will be parametrized by $h=H_{0}/(100\,\mathrm{km\,s}^{-1}\mathrm{Mpc}^{-1})$ in our test simulations), and $a(|\mathbf{x}^{\prime}|/\mathrm{c})=a(\tau)$ denotes the scale factor at the time $\tau$ corresponding to LOS distance $r^{i}(x^{\prime})$.

V Implementation

The implementation, not only of the data model, but of the entire algorithm, is based on NIFTy (Selig et al., 2013), a versatile software package for the development of inference algorithms. We further compute cosmology-dependent quantities, like power spectra and distance-redshift relations, with the publicly available CLASS code (Blas et al., 2011). To summarize the data model we introduce short-hand notations for each operation in terms of operators.

In its most general form, the prior and data model, that connect the Gaussian field $\varphi$ with a data vector of $N_{s}$ measured source ellipticities, are as follows: The growth operator $D_{\mathbf{x}\mathbf{y}}(\tau)$ imprints a growth structure on the Gaussian field. The resulting four-dimensional field is plugged into the exponential of $E(\cdot)=\exp(\cdot)-1$ [see Eq. (4)] to obtain the fractional overdensity $\delta_{\mathbf{x}}(\tau)$. The overdensity induces the potential $\phi_{\mathbf{y}}(\tau)$ by the Poisson equation, encoded in the operator $\mathcal{P}_{\mathbf{y}\mathbf{x}}(\tau)=\Delta^{-1}\frac{3}{2}\Omega_{\mathrm{m}}H_{0}^{2}\delta(\mathbf{x}^{\prime})/a(|\mathbf{x}^{\prime}|/\mathrm{c})$ [Eq. (19)]. The potential can be computed efficiently in Fourier space. The gravitational potential is then restricted to the light cone of the observer by the light cone operator $\mathcal{C}_{\mathbf{y}^{\prime}(\mathbf{y},\tau)}$ [Eq. (1)]. We compute the tidal tensor of the resulting 3D field [(Eq. (18)] by application of a global differential operator, which we denote $\mathcal{T}_{\mathbf{z}^{\prime}\mathbf{y}^{\prime}}$. The resulting tidal tensor is then rotated into each galaxy’s LOS coordinate system by a rotation operator, $\mathcal{R}^{i}_{\mathbf{r}\mathbf{z}^{\prime}}$. An integration operator $\mathcal{I}_{j\mathbf{r}}$, which applies the integration in Eq. (16), integrates the components of each of the resulting $N_{s}$ tidal tensors along the unperturbed photon geodesic. For this operator, we use an adapted version of the implementation that was already successfully used in a similar reconstruction method (Greiner et al., 2016). The application of $\mathcal{I}_{j\mathbf{r}}$ yields derivatives of the lensing potential for each galaxy location. From this we can compute the shear components by a linear operator $\mathcal{L}_{ij}$ that comprises the equations in Eq. (14). Rotation and integration map the three dimensional continuous signal space into the discrete space (one point for every galaxy) of the data. The shear components are thus automatically computed at the locations of the galaxies.

In the simplified implementation that we use for this work, we avoid the 4D coordinate grid $(\mathbf{x},\tau)$ and work on the three-dimensional light cone from the beginning. In the prior, we model the Gaussian log-density $\varphi$ with the power spectrum of matter fluctuations today $P_{\mathrm{lin}}(k,a{=}1)$, where $a$ denotes the scale factor. The growth operator is diagonal in configuration space and only a function of comoving distance to the observer $D_{\mathbf{x}^{\prime}\mathbf{y}^{\prime}}=D^{(1)}(\tau)\delta_{\mathrm{D}}(\mathbf{x}^{\prime}-\mathbf{y}^{\prime})$, where $D^{(1)}(\tau)$ is the growing solution of the linearized growth equation. The Poisson operator is split into two parts. First, a multiplication with $3/2\,\Omega_{\mathrm{m}}H_{0}^{2}/a(|\mathbf{x}^{\prime}|/\mathrm{c})$, i.e. an operation that is diagonal in configuration space. Second, the inverse Laplace operation $\Delta^{-1}$, which is diagonal in Fourier space. The inverse Laplacian in the Poisson equation is a non-local operation that should strictly be applied to 3D spatial volumes at fixed time. Here, we apply it on the light cone noting that the induced error in radial direction should be small (roughly of order $a^{2}$) if we apply the first part of the Poisson operator first and if $D^{(1)}(\tau)$ is approximately proportional to the scale factor. In this case, the first order term of the exponential expression

[TABLE]

and therefore the first order time dependence of the overdensity field will partly be canceled by the $1/a$ in the Poisson equation before the inverse Laplacian rescales the field. This cancellation corresponds to the commonly known fact that the comoving gravitational potential is constant in a matter dominated Universe.

After computing the six independent entries of the tidal shear tensor, we integrate each of its components along the LOS. Only after the integration we rotate the resulting tensors into their LOS coordinate systems and project out all entries which are not in the plane perpendicular to $\mathbf{r}$. This change of order allows to efficiently combine three operations, that is the rotation into the LOS, projecting out non-perpendicular components of the tidal tensor, and the computation of shear components and convergence, in a single linear operation. We denote the corresponding operator $\mathcal{G}$ (Gamma-Projection-Rotation).

For one data point the implemented data model in operator notation reads

[TABLE]

The total data vector ${d}$ has dimensions $2\times N_{s}$, i.e. two ellipticity components for each of the $N_{s}$ source galaxies. We use the letter ${R}$ to encode the total response of the data to a signal $\varphi$, i.e. the composite action of all operators. We have also added an experimental noise $\mathbf{n}$ here for completeness. In general this will be subdominant to the shape noise $\boldsymbol{\epsilon}^{\mathrm{s}}$ and we ignore it in the following. Note, however, that the formalism allows for an incorporation of several independent noise sources.

VI MAP estimator

Our aim is to obtain a maximum a posteriori (MAP) estimate of the signal field $\varphi$. The posterior distribution is related to the likelihood and the prior by Bayes’ theorem

[TABLE]

The prior probability $P(\varphi)$ is modeled as a Gaussian distribution with covariance $\Phi$. In a first simple approximation $\Phi$ can be taken to be diagonal in Fourier space with the usual power spectrum of matter overdensities $P_{\mathrm{lin}}(k,a)$.

To obtain the likelihood, we marginalize over all realizations of the shape noise

[TABLE]

where the covariance $N$ of the shape noise was defined in Eq. (11) and we neglect any other sources of measurement noise. With this, the negative log-posterior becomes

[TABLE]

where we have dropped most terms that are independent of the field of interest, $\varphi$.

The maximum of the posterior distribution is found by minimizing the expression in Eq. (24). Note that the posterior distribution is not Gaussian. Due to the exponential in the response ${R}$ it is not quadratic in $\varphi$. To find the minimum of the negative log-posterior we apply a Newton-like minimization scheme (Byrd et al., 1995). This requires the derivative of the negative log-posterior with respect to $\varphi$

[TABLE]

where the star denotes a pointwise product in position space, i.e. $(\alpha*\beta)_{\mathbf{x}}=\alpha_{\mathbf{x}}\beta_{\mathbf{x}}$.

VII Redshift-marginalized likelihood and posterior

We can take into account the source redshift uncertainty by generalizing the marginalized likelihood in Eq. (23) to include the probability of the source redshifts $z_{s}$ given their measured photometric redshifts $z_{\mathrm{photo}}$. This probability is given by the posterior redshift distribution function $P(z^{s}|z_{\mathrm{photo}})$, where $z^{s}$ and $z_{\mathrm{photo}}$ are vectors of the redshifts of all sources. Including this redshift posterior in the likelihood yields,

[TABLE]

In most cases, this integration cannot be done analytically and the resulting distribution will in general not be Gaussian. In the absence of further information about the general shape of this distribution, a Gaussian characterized by the first and second central moment of the full distribution is the most conservative approximation one can use (it is the maximum entropy approximation and therefore adds the least additional information). This motivates the use of a Gaussian to model the redshift-marginalized likelihood.

To obtain the first and second moment, we need to calculate $\langle d\rangle_{P(d|\varphi)}$ and $\langle dd^{\dagger}\rangle_{P(d|\varphi)}$. For the expectation value of the data, we obtain

[TABLE]

where $\bar{I}$ denotes the redshift averaged integration operator defined in the square brackets in the first line of Eq. (28) and we have introduced $\bar{R}(\cdot)\equiv\mathcal{T}\mathcal{P}\left(\mathrm{exp}\left[D(\cdot)\right]-1\right)$ to summarize the action of all source-redshift independent operators.

The second moment is

[TABLE]

Non-diagonal elements of the second term in Eq. (29) read

[TABLE]

and diagonal elements are

[TABLE]

where the new operator $\tilde{\mathcal{I}}$ denotes the squared average of the integration operator, i.e. the square brackets in the first line of Eq. (31).

The Gaussian approximation to the likelihood is then $\mathcal{N}[d|\tilde{R}(\varphi),\tilde{N}]$, where $\tilde{N}=N+Q$ and $Q$ is

[TABLE]

This expression is still signal-dependent and we approximate it further by replacing $\varphi$ by its posterior mean $\langle\varphi\rangle_{P(\varphi|d)}$. Since this mean depends on $\tilde{N}$, the resulting set of equations must be solved iteratively.

To estimate the effect of the additional redshift-marginalization, we compare the effective lensing efficiency that appears in the redshift-marginalized data model (Eq. 28),

[TABLE]

to the exact lensing efficiency $W(r,r^{s})$.

In the left panel Fig. 1 we plot both kernels for a source at redshift $z^{s}=0.8$ and a Gaussian redshift posterior with variance $\sigma_{z}=0.04(1+z^{s})$. For this particular configuration we find that both integration kernels are very similar over most of the distance to the source, i.e. their relative difference never exceeds 4% for $r<0.9r^{s}$. Only very close to the galaxy the kernels start to differ significantly. At this point the amplitude of the kernels has already dropped by more than $80\%$. The contribution to the lensing signal from structures with $r>0.9r^{s}$ is therefore strongly suppressed.

In the right panel of Fig. 1 we plot the relative difference of the integration volumes of the kernels assuming a uniform density contrast along the line of sight and the same Gaussian redshift posterior with $\sigma_{z}=0.04(1+z^{s})$. We find that the difference between the kernels is sub-percent for sources with $z^{s}>0.43$.

These tests suggest that the redshift uncertainty should change the lensing signal of individual galaxies on average at the percent level.

We caution, however, that these tests are very simplified and only serve to get a rough estimate of the effect of redshift uncertainties, especially since we have ignored catastrophic outliers or the fact that redshift posteriors can be bimodial. For the proof of the concept of our reconstruction algorithm, we ignore redshift uncertainties. Following a test of their exact impact (e.g. by comparing a reconstruction with exact spectroscopic redshifts to one with redshifts obtained from the corresponding photometric redshift posteriors), they could be included in the way that has been outlined in this section.

VIII Validation and Tests

To validate the implementation and assess the goodness of the tomographic reconstruction, we perform a number of increasingly realistic tests, which we denote tests A, B and C. In tests of type A we employ an even source distribution over the entire box. These tests are the least realistic ones and serve to validate the correctness of the implementation. Test B uses a survey-like source-redshift distribution and test C adds realistic shape noise.

In all of these tests, we place the observer in the center of the bottom of the computational box222This position is not fixed by the algorithm. For a reconstruction from a survey that covers a significant fraction of the sky, the observer can be placed in the center of the box, for example. and resolve the underlying and reconstructed overdensity fields with $128^{3}$ pixels. Depending on the test, we allow the physical sizes of the box to differ. The current pixel resolution is limited due to computation time and memory usage. A higher resolution will be accessible after parallelization and adaption of the code for the usage on a high-performance cluster.

The minimization is achieved through 300 steps of an LBFGS algorithm Byrd et al. (1995) followed by a steepest descent algorithm. The latter is set to have reached convergence if the maximal pixel-wise relative difference between two subsequent field estimates is smaller than $10^{-4}$ for three iteration in a row. If the steepest descent does not converge within 200 iterations, the same combination of LBFGS and steepest descent minimizers is repeated. We find that, at latest in this second run, all of our example reconstructions meet the convergence criterium. In the current serial implementation, the entire minimization takes between 1 day (for idealized tests at fixed redshifts) and 6 days (for a realistic source distribution and shape noise) on a 2.3 Ghz core of an Intel Xeon E5 processor.

Most of the tests are based on mock data that we create by applying the data model described in Sec. IV to non-linear density fields obtained from N-body simulations. For tests at fixed redshifts we use snapshots of the Millennium-XXL simulation Angulo et al. (2012) from which we take the entire volume of size $[3h^{-1}\,\mathrm{Gpc}]^{3}$. For more realistic test cases we construct a light cone of size $[500h^{-1}\,\mathrm{Mpc},500h^{-1}\,\mathrm{Mpc},4000h^{-1}\,\mathrm{Mpc}]$ by joining snapshots of the Millennium Run Springel et al. (2005). The Millennium box measures $500h^{-1}\,\mathrm{Mpc}$ along each side and we shift it every time $500h^{-1}\,\mathrm{Mpc}$ in the $z$-direction have been constructed. This procedure ensures that a LOS is unlikely to hit the same structure repeatedly. In the resulting light cone we achieve a resolution of $3.9h^{-1}\,\mathrm{Mpc}$ in the $x$- and $y$-directions corresponding to mildly non-linear scales and reach a redshift of $z=2.2$ in the $z$-direction. Note that the physical size of the pixels is eight times longer along the $z$-axis, meaning that we obtain a poorer resolution in radial direction.

We map the simulations onto our pixel grid by a nearest grid point interpolation. This procedure suppresses power on small scales and modifies the power spectrum by a kernel consisting of cardinal sine functions (see e.g.Jing (2005)). To be consistent with this modification, we apply the same kernel to the theory power spectra which we use in the prior distributions.

Both the Millennium-XXL and the Millennium simulation use a flat $\Lambda$CDM cosmology with $\Omega_{\mathrm{m}}=0.25$ and $h=0.73$.

VIII.1 Fidelity of the Reconstruction

The BFGS algorithm is a quasi-Newton method. While Newton methods rely on an exact evaluation of the inverse Hessian at each step, quasi-Newton methods estimate the curvature from previous evaluations of the gradient. The ’L’ in L-BFGS refers to the limited memory version of this algorithm in which the inverse Hessian is approximated by a sparse matrix instead of a dense one. This additional approximation becomes necessary for problems with many degrees of freedom such as this. In principle the final estimate of the inverse Hessian could be used as an approximation for the pixel-pixel covariance, but it depends on the correlated previous steps of the algorithm and only probes the curvature along the search directions. The covariance, $C$, of a Gaussian posterior distribution, $\mathcal{G}(\varphi|d)$, with maximum $m$ can be estimated by sampling. A straightforward algorithm to obtain such an estimate is outlined in the following.

Draw random realizations of signal and noise from the Gaussian prior and noise distributions.

[TABLE] 2. 2.

Use the resulting realizations to model a data vector $d^{\prime}$

[TABLE] 3. 3.

From this random data vector reconstruct the Wiener filter estimate $m^{\prime}$, i.e. the maximum of the Gaussian posterior distribution,

[TABLE] 4. 4.

Repeat steps 1.-3. $N$ times and estimate $C$ from

[TABLE]

Computing the full covariance is computationally infeasible since it would require $\sim 9$ Tb to store a full covariance matrix of size $1/2n(n+1)$ for a map with $n=128^{3}$ pixels. Numerically feasible are, e.g., estimates of its diagonal (the pixel-wise covariance), or of the full covariance of 1D slices of the reconstruction box.

The above algorithm yields an unbiased estimator of the posterior covariance only if the data model is linear in the density field. The lognormal model is inherently non-linear and therefore the posterior non-Gaussian. However, for sufficiently small perturbation $\Delta m$ around the MAP estimate $m$ the Laplace approximation should be applicable, meaning that we can approximate the negative log posterior as a Gaussian,

[TABLE]

To obtain Eq. (38), we have linearized the response operator

[TABLE]

redefined the data vector $\tilde{d}=d-R(m)+R^{\prime}(m-\bar{\varphi})$ and shifted the signal $\tilde{\varphi}=m+\Delta m-\bar{\varphi}$. The covariance $C^{\prime}$ of this Gaussian approximation to the posterior is

[TABLE]

This approximate covariance can be evaluated as outlined above. As an example, we estimate $\mathrm{diag}({C^{\prime}})$, which can be interpreted as an uncertainty map for the reconstruction, for one of our test cases in the next section.

VIII.2 Tests A: Simple geometry

In a first series of tests we distribute 500 000 sources evenly in the box and slightly beyond ($\pm 100h^{-1}\,\mathrm{Mpc}$). We place sources outside of the reconstructed volume in order to increase the area in which we can recover the underlying density field with high resolution. The quality of the reconstruction decreases with increasing distance to the observer because of increasing distances between lines of sight and a decreasing number of background sources that help to break the LOS degeneracy of the lensing kernel. We further only add unphysically low, negligible shape noise in this first test set.

We perform three different tests of type A for which we create data from increasingly realistic input fields:

1. A self-consistency test, where we create an overdensity field from a random realization of the Gaussian field $\varphi$ [Eq. (3)] in a cubic box of side length $1000h^{-1}\,\mathrm{Mpc}$. We then apply the implemented data model and use the algorithm to recover this input field. We do not show any results from this test, since they would not provide any additional insights compared to tests A2 and A3, but simply state that we can recover the input field with high fidelity which means that the implementation is in itself consistent.

2. A test in which we create shears from overdensity fields taken from an N-body simulation at fixed redshifts. We use three different snapshots at individual redshifts $z=1$, $z=0.25$ and $z=0$. This allows to assess the algorithm’s ability to recover increasingly non-linear fields. In each case, we compare the lognormal reconstruction to a Wiener filter (WF) reconstruction, which uses the same data model and source distribution but a Gaussian prior on the overdensity field $\delta$. Comparisons between input and reconstructed fields in each case are shown in Fig. 2. The resolution of the reconstruction decreases with distance to the observer. This is because 1) the density of lines of sight decreases 2) there are less sources behind the point we want to reconstruct 3) the information from these sources is suppressed by the shape of the integration kernel. In Table 1 we compare the minimal and maximal values in the reconstructions to the extremal values in the true density field. The lognormal reconstruction is superior in capturing the highest values of the density field and avoids unphysically low density contrasts below -1. Table 2 summarizes pixel-wise quadratic differences and the pixel-wise Pearson correlation between underlying and reconstructed fields. For all redshifts tested the lognormal prior yields higher correlations with the original field than the Gaussian prior. The difference between lognormal and WF reconstruction increases as the input field becomes more non-linear. This is also reflected in the 1-point probability distribution functions (PDFs), which we show in Fig. 3. While the Wiener Filter PDF is closer to symmetric in all cases, we can capture part of the skewness of the input field by applying the lognormal prior. An notable feature of the 1-point PDFs is that the maximum value is slightly biased in both reconstructions. The distributions of the reconstructions peak at zero, i.e. at their mean, while the underlying field peaks below. This is a feature of both priors since they prefer the mean density (or a zero density contrast) if the data does not contain enough information on the density.333To be precise one should note that the maximum of the lognormal posterior lies slightly below the mean density. However, the difference to the mean turns out to be small in all test cases. For the light cone box we have calculated that the prior distribution peaks at $\delta=-0.1$. This value is consistent with the lognormal reconstructions that we obtain in our test cases. Closely inspecting the right panels of Figs. 5 and  10, one can see that the bulge of the reconstructed densities lies slightly below $1$, at $\delta+1\approx 0.9$).

3. A test similar to A2 but with a redshift dependent density field constructed from different snapshots. The physical size of the box spans $[500h^{-1}\,\mathrm{Mpc}]^{2}$ in the $(x-y)$-plane and $4000h^{-1}\,\mathrm{Mpc}$ into the $z$-direction corresponding to a maximum redshift of $z=2.2$.444The light cone is constructed by merging planes from different snapshots along the $z$-direction of the box, but the algorithm assumes that the distance and redshift increase in radial direction from the observer. This leads to a mismatch between the distance to a (source) position assumed by the algorithm and the distance or rather redshift that this position corresponds to in the simulation. The distance in the simulation is $\cos\theta$ times the distance assumed by the algorithm, where $\theta$ is the angle between the LOS and the $z$-axis. This mismatch should only become relevant for relatively large angles $\theta>\pi/2$. Lines of sight with such angles only cover a minor fraction of the box and we therefore expect the test to not be severely affected by this approximation. Results of test A3 are shown in Fig. 4 and Fig. 5. The quantitative comparison between the lognormal and the corresponding WF run is listed in the last rows of Tables 1 and 2, labeled “light cone”. We find a similar match between the underlying and reconstructed fields as in tests A2 and again a superiority of the lognormal prior over the Wiener Filter reconstruction in terms of correlation coefficients, pixel-wise differences and extremal values. The distributions of the pixel-wise matching between input and reconstructions in the left panel of Fig. 5 imply that the lognormal prior is more likely to reconstruct overdensities correctly, while the Wiener Filter is better in tracing the underdensities. The full distributions in the right panel of the same Figure suggest that the better reconstruction of low densities with the Wiener Filter is a consequence of the symmetry of the Gaussian prior: The Gaussian prior suppresses both positive and negative offsets from the mean density by an equal amount. Because of this property it also allows negative densities. The lognormal prior strongly suppresses deviations towards low densities, while it is less sensitive to deviations towards higher densities (since it suppresses deviations equally on a logarithmic scale). As a result, the lognormal model refrains from assigning very low densities to regions with small signal-to-noise.

In Fig. 7 we show estimates of the reconstruction uncertainties in this test which were obtained by following the prescription in Sec. VIII.1. We find that the estimator in Eq. 37 converges very quickly and use 20 samples (from the 15th sample on the average change per pixel from adding another sample is less than 0.05, adding the 20th sample changes the uncertainty estimate on average only by 0.03 per pixel). For both priors, but especially for the lognormal prior, a comparison to the input field shows that overdensities are reconstructed with higher fidelity than low density regions. The Wiener filter reconstruction is less sensitive to overdensities but slightly more sensitive to underdensities than the lognormal reconstruction.

In Fig. 6 we compare the correlation coefficient between input field and reconstructions on different scales. The Gaussian model provides better fits on the largest scales. This could be due to the aforementioned missing power in low density regions in the lognormal reconstruction, which can be best seen by eye in the lower row of Fig. 4. On intermediate and small scales the power spectrum of the lognormal reconstruction shows higher correlations with the original field. We note that the estimator in this work was designed for the map level and is not expected to be an optimal estimator for the power spectrum.

VIII.3 Test B: Realistic survey geometry

In the test of type B, we employ a realistic survey geometry and source distribution. We use input overdensities from the same light cone as in test A3 to generate mock data and place sources in a cone that spans $7.15^{o}$. The sources are distributed according to a distribution function of form

[TABLE]

where we fit $\alpha$, $\beta$, and $z_{0}$ to the publicly available source distribution of the CFHTLenS survey555http://www.cfhtlens.org/astronomers/cosmological-data-products. The resulting source distribution function is shown in the left panel of Fig. 8. The fit does not exactly match the survey distribution, but the similarity is sufficient for this test. We then draw $\Omega[\mathrm{armcin}^{2}]\rho_{s}[\mathrm{gal}/\mathrm{arcmin}^{2}]$ source positions from this distribution, where $\Omega$ is the angular opening area of the cone in $\mathrm{arcmin}^{2}$ and $\rho_{s}$ the source density of the survey, which we choose to be $\rho_{s}=11\,\mathrm{gal}/\mathrm{arcmin}^{2}$ corresponding to the source density in CFHTLenS. A scatter plot of the resulting spatial source distribution is shown in the right panel of Fig. 8. The reconstructions are depicted in Fig. 9.

For a realistic source distribution the lognormal model seems to better reconstruct the underlying field at lower redshifts, while the Gaussian prior better traces distant low density regions. We interpret this is again a consequence of the symmetry of the Gaussian distribution, the same property which leads to the abundant reconstruction of negative densities, nicely visible in the 1-point PDF (Fig. 10).

VIII.4 Test C: Realistic shape noise

In the last test of our test series, we add realistic shape noise with $\sigma_{\epsilon}=0.25$ to the mock data on the light cone. This value is close to the value in CFHTLenS, $\sigma_{\epsilon}=0.279$ Heymans et al. (2013). Compared to test B, we use the same opening angle of $7.15^{o}$ but a slightly higher galaxy density of $15\,\mathrm{gal/arcmin^{2}}$. We show results of this final test in Figures 11 and 12. The left panel of Figure 11 shows that the lognormal model is more likely to reconstruct the correct density in a pixel. However, the PDF in the right panel seems to have a spurious tail of high density values. This tail does not significantly deviate from the high density tail of the underlying reconstruction and is therefore hardly suppressed by the prior. Figure  11 reveals that these high density values are due to an overestimation of densities in slices normal to the direction of observation. We have verified that these overdense sheets are not due to a wrong minimum or saddle-point by starting the algorithm at either the WF solution or the true input field and changing convergence criteria. We interpret these features as follows: Adding a constant density offset to a slice normal to the observational direction does not significantly change the likelihood, because of the well known mass-sheet degeneracy - the fact that adding a homogeneous mass sheet does not change the lensing signal. The lognormal prior does not punish or even encourages these overdensities since it expects a skewed density distribution. The Gaussian prior suppresses strong deviations from a symmetric density distribution, which is why these features are not present in the Gaussian case.

IX Dicussion and Future Improvements

We have presented a novel, fully Bayesian, lensing tomography algorithm that reconstructs the three-dimensional matter distribution from measurements of individual galaxy shapes. Given a data vector of measured galaxy shears, the galaxies’ angular positions and their redshifts, this algorithm finds the maximum a posteriori estimate of the density contrast between the lensed galaxies and the observer. Since galaxy shapes provide only limited information about the intervening structure, this reconstruction problem is in general ill-constrained and requires regularization. In case of a maximum a posteriori estimate this regularization is provided by the prior probability distribution - the distribution which ideally contains all available information about the underlying field before the lensing measurement.

This ideal distribution, however, is in reality unaccessible. Similar previous works have used a Gaussian prior in this context. A Gaussian prior (together with a Gaussian noise model) results in a Gaussian posterior distribution which has a distinct minimum and no saddle points. Maximizing this posterior is comparably straightforward. The algorithm presented in this work was especially designed to also handle non-Gaussian posteriors. We show this ability by using a lognormal prior which should be an improvement over the Gaussian approximation in two ways: First, it enforces the strict positivity of the density while a Gaussian prior allows unphysical, negative densities. Second, it incorporates an a priori knowledge of the presence of odd moments in the matter distribution, which arise as a consequence of non-linear structure formation. These corrections are relevant since the cosmic shear signal probes structures down to scales that lie well within the non-linear regime.

We note that neither of these priors describe the exact distribution of the non-linearly evolved density contrast. In regions with high signal-to-noise this is not problematic since the prior is updated by the information contained in the data. In less well constrained regions the prior will paint in missing information. This can lead to systematic deviations from the true distribution. However, non optimal priors can be sufficient depending on the quantity and scales one is interested in reconstructing. For example a Gaussian prior will be a good choice for large scales, where the field is still Gaussian. On small scales the lognormal model will yield better reconstructions since it is informed about the skewness of the distribution and avoids negative densities.

To test the algorithm and assess the differences between a reconstruction with a lognormal and a Gaussian prior (in the latter case the reconstruction simplifies to a Wiener filter), we have applied it in both configurations to increasingly realistic mock data sets. Those data sets were produced by applying the algorithm’s inherent data model to non-linearly evolved densities from N-body simulations. In tests with negligible shape noise, but a realistic CFHTLenS-like source distribution, we find that the reconstruction with a lognormal prior avoids negative densities that are abundantly present in the Wiener filter reconstruction. It is also better in reconstructing the highest local peaks in the density field, and it leads to a higher point-wise correlation between the true underlying density and its reconstruction. The Gaussian model seems better in tracing low density regions and capturing the largest scales. In regions with low signal-to-noise both priors tend towards the mean density.

For realistic shape noise, the data contains too little information to break the line-of-sight degeneracy of the tomographic problem with either prior. This is expected for current data sets, which are commonly analyzed in a limited number of broad redshift bins in order to yield sufficient signal to noise per bin. In addition, the lognormal model is prone to mass-sheet degeneracy: since the prior encourages a skewed distribution with more high density pixels, the reconstruction tends towards inserting sheets of higher density normal to the observational direction. The situation might slightly improve for future surveys with a higher source density. But since the shape noise scales only with the square root of the source density, even future data is unlikely to contain sufficient information for a 3D reconstruction without this degeneracy.

The fidelity of the reconstructions in this work might also be overestimated since we assumed precisely known source redshifts and point spread function, and ignored intrinsic alignments. In this work we show how photometric redshift uncertainties can be included, but leave the implementation for future work. Other sources of uncertainty could be included in a similar fashion. Another approximation we made concerns the source distribution. In reality the source distribution is correlated with the underlying density field, which we ignore in our test cases. A more clustered distribution of sources would render the fidelity of the reconstruction more inhomogeneous since structures in front of high density regions would be detected with higher signal to noise than structures in front of low density regions. Since the volume fraction of high density regions is smaller than the volume fraction of low density regions, this could decrease the overall fidelity of the reconstruction. To get an estimate of the effect of the clustering, future tests of the algorithm should account for the correlation between galaxies and density field by using the halo positions identified in the simulation by a halo finder.

Possible improvements of the algorithm could go into two directions, 1) improving the prior or data model, 2) including additional sources of information. An improved data model could be obtained by forward modeling structure formation numerically for an initially Gaussian distributed field. While this is conceptually straightforward to integrate in the algorithm, it is numerically very expensive. In general, we expect any prior model that is informed about the relative orientation of structures (as e.g. the Zeldovich approximation) to partly break the line-of-sight degeneracy or inhibit spurious mass sheets. Additional information could come from the clustering of galaxies. This would require the inference of at least one additional bias parameter (depending on the bias model), but would add information about the relative orientation of structure, also possibly breaking the line-of-sight and mass-sheet degeneracy.

In the current implementation the algorithm assumes a fixed set of cosmological parameters for both models. For future applications the parameter estimation could be made part of the algorithm by analyzing the joint posterior of the density distribution and cosmological parameters. The cosmology enters the density estimate through the theoretical matter power spectrum, the growth factor and the background evolution (that determines e.g. the redshift-distance relation). Lensing measurements are mostly sensitive to a combination of $\sigma_{8}$ and $\Omega_{m}$. For two cosmological parameters only, sampling of the joint posterior would be feasible already with the current serialized implementation, for one could sample the $\sigma_{8}$-$\Omega_{m}$ plane (e.g. on a grid) and evaluate the minimum of $P(\delta|\sigma_{8},\Omega_{m},d)$ for each parameter combination in parallel. Since the total information content of the data does not depend on the analysis method, a joint reconstruction of the 3D density field and cosmological parameters should yield competitive constraints to traditional analyses that make use of redshift binning.

Acknowledgements.

We thank Jakob Knollmüller for helpful insights on covariance probing and Reimar Leike and Steffen Hagstotz for useful comments on the draft.