Bayesian weak lensing tomography: Reconstructing the 3D large-scale distribution of matter with a lognormal prior Vanessa B¨ohm,1 Stefan Hilbert,2,3 Maksim Greiner,1,2 and Torsten A. Enßlin1,3,2 1Max-Planck-Institut fu¨r Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany 2Exzellenzcluster Universe, Boltzmannstr. 2, 85748 Garching, Germany 3Ludwig-Maximilians-Universita¨t, Universita¨ts-Sternwarte, Scheinerstr. 1, 81679 Mu¨nchen, Germany (Dated: January 10, 2017) WepresentaBayesianreconstructionalgorithmthatinfersthethree-dimensionallarge-scalemat- terdistributionfromtheweakgravitationallensingeffectsmeasuredintheimageshapesofgalaxies. Thealgorithmassumesthatthepriorprobabilitydistributionofthematterdensityislognormal,in contrast to many existing methods that assume normal (Gaussian) distributed density fields. We 7 compare the reconstruction results for both priors in a suite of increasingly realistic tests on mock 1 data. Wefindthatincasesofhighnoiselevels(i.e. forlowsourcegalaxydensitiesand/orhighshape 0 measurement uncertainties), both normal and lognormal priors lead to reconstructions of compa- 2 rable quality. In the low-noise regime, however, the lognormal model produces significantly better reconstructions than the normal model: The lognormal model 1) enforces non-negative densities, n whilenegativedensitiesarepresentwhenanormalpriorisemployed, 2)bettertracestheextremal a J values and the skewness of the true underlying distribution, and 3) yields a higher correlation be- tween the reconstruction and the true density distribution. Hence, the lognormal model is to be 7 preferredoverthenormalmodel,inparticularsincethesedifferencesbecomerelevantfordatafrom current and futures surveys. ] O C . I. INTRODUCTION h p - WeakgravitationallensingofgalaxiesoffersauniquewaytostudythedistributionofmatterintheUniverse(see [1] o forareviewonweakgravitationallensingand[2–4]forrecentreviewsonweakgalaxylensing). Lensingbystructures r t along the line of sight causes distortions in the images of distant galaxies (which in this context are often referred to s a assources), which leadstocorrelations betweentheapparentshapesofthese galaxies. The dominantandmosteasily [ 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. 1 v Galaxy shape measurements allow to constrain the clustering of matter at different scales and redshifts, which can 6 then be translated into constraints on cosmological models and their parameters. The integrated lensing signal is 8 mostlysensitivetoacombinationofthecosmicmeanmatterdensityΩ andthematterpowerspectrumamplitudeσ . m 8 8 Tomographic methods yield additional constraints on the properties of dark energy and on modified-gravity theories. 1 Since lensing is a direct probe of the total matter, luminous and dark, it can be combined with measurements of the 0 luminous matter distribution in order to learn about the relationship between baryons and dark matter. Another . 1 importantfeatureofweakgalaxylensingliesinitsabilitytoprobethematterdistributionoverawiderangeofscales, 0 from many tens of Mega parsec, where structure formation is still linear today and comparably easy to model, down 7 to non-linear sub-Mega parsec scales. Due to its sensitivity to such a wide range of scales, lensing can provide a 1 largeamountofinformationtoconstrainmodelsofnon-linearstructureformationandcosmology,inparticularifalso : v higher-order statistics are considered. i X Firstfirmstatisticaldetectionsofcosmicshearwerereportedin2000byfourdifferentgroups[5–8]. Sincethen, the field has seen a tremendous increase in the amount and quality of lensing data as well as a notable improvement in r a analysistechniques. Acommonmethodtoinfercosmologicalparametersfromweakgalaxylensingdataistocompare thepowerspectrumofthefullyprojected2Dshearfieldtotheoreticalpredictions. Strongerconstraintscanbeobtained by measuring this power spectrum in a number of redshift bins [9]. This compressed tomographic procedure allows to better assess the time evolution of structure formation and thus to tighten constraints on parameters affecting this evolution. Several authors have investigated the additional constraining power that can be achieved by incorporating third-orderstatisticsand/orshearpeakcountsandcorrelations[10–19],whichhelpstobreakparameterdegeneracies. 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 [20–22] and were recently applied to data from the Canada France Hawaii Telescope Lensing Survey (CFHTLenS) [23]. 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 toobtainthedensityhavebeenworkedoutbyanumberofauthors[24–26]. Sincethisinversionisunder-constrained, itrequiressomeregularizationmethodorchoiceofprioronthedensityfield. MostofthealgorithmsemployaWiener 2 filter, which corresponds to a normal (Gaussian) prior, that can be complemented with information about galaxy clustering [27, 28]. Weaimtoextendtheworkontomographicreconstructionofthe3DmatterdistributionbydesigningafullyBayesian reconstructionalgorithmthatusesalognormalprioronthedensityfield. Thealgorithmisdesignedtoreconstructthe 3Dcosmicdensityfluctuation fieldδ(x,τ)fromweakgalaxy lensing data, i.e. ameasurementofgalaxyellipticitiesat different(photometricallymeasured)redshifts. Itsderivationisbasedonthelanguageofinformationfieldtheory[29], which has already been used to address similar tomography problems [30]. 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 auto- matically 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 [31]. Characterizing the matter overdensities in the Universe as a lognormal field was first assessed by Coles & Jones in 1991 [32]. 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 [33, 34]. 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 [29, 35, 36]. Lognormal distributions have also already been considered in the context of weak lensing: Analyses of ray-tracing simulationsandtheDarkEnergySurvey(DES)ScienceVerificationdatashowedthatthe1-pointdistributionfunction of the lensing convergence is better described by a lognormal than a Gaussian model [37, 38]. 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 [39]. Bayesian inference methods are widely used in weak shear analyses, most prominently in the context of shear measurements from galaxy images [40–42]. 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 pa- rameters [43, 44]. 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. InSec.VIII,weshowresultsofthedensityreconstructiononincreasinglyrealisticmockdata. Weconcludethiswork 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 x=(x ,x ,x ) at fixed comoving lookback time τ, that, combined 0 1 2 with the time coordinate, form the four-dimensional coordinate system (x,τ). 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. x(cid:48). Since x(cid:48)-coordinates implicitly define a comoving lookback time τ(x(cid:48))= x(cid:48) /c(wherecdenotesthespeedoflight),weomitspellingoutthetimeexplicitlyandwriteA(x(cid:48))=A(x(cid:48),τ) | | 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 C =δ (x(cid:48) x)δ (τ x/c). (1) x(cid:48)(x,τ) D D − −| | 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 (ˆri,ˆri,ˆri), whereˆri(x(cid:48)) points into the direction of the ith source galaxy and the normal 0 1 2 0 vectors (ˆri,ˆri) span the two-dimensional plane perpendicular to ˆri. The radial comoving distance of each galaxy 1 2 0 from the observer is denoted ri = ri ri . The transformation from the light cone system into the line of sight (LOS)systemofasourcegalaxyiis|ac|h≡iev|ed0|byarotationofthex(cid:48)-axisintotheithlineofsight,whilex(cid:48) andx(cid:48) get 0 1 2 aligned with the image coordinates of the galaxy observation. We label the corresponding transformation operators irx(cid:48). R 3 III. PRIOR MODEL The comoving density ρ(x,τ) can be split into its time-independent spatial mean ρ¯, and a perturbation δρ(x,τ)= ρ¯δ(x,τ). The fractional overdensity δ(x,τ) is commonly modeled as a homogeneous, isotropic Gaussian field with zero mean and power spectrum P (k,τ). 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 δ 1, the time evolution of the density field can be described by (cid:28) (cid:20) (cid:90) (cid:21) ρ(x,τ)=ρ¯[1+δ(x,τ)]=ρ¯ 1+ D(x y,τ,τ )ϕ(y,τ ) =ρ¯[1+D (τ,τ )ϕ (τ )]. (2) 0 0 xy 0 y 0 | − | y 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(x y,τ,τ ) is the 0 | − | integration kernel of a linear homogeneous and isotropic, but possibly scale-dependent, growth operator. The field ϕ is an isotropic and homogeneous random field whose values are drawn from a multivariate normal distribution with mean ϕ¯ and covariance Φ, ϕ (cid:45) (ϕ ϕ¯,Φ). (3) ← N | Here, it describes the three dimensional density fluctuations at time τ and is translated to other times τ by the 0 growth operation D (τ,τ ). This implies ϕ¯=0 for this linear Gaussian model. xy 0 The description Eq. (2) breaks down when δ 1 such that non-linearities become important. A possible way to (cid:54)(cid:28) 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: ρ(x,τ)=ρ¯[1+δ(x,τ)]=ρ¯exp[D (τ,τ )ϕ (τ )]. (4) xy 0 y 0 Since the expectation value of the density ρ(x,τ) must equal ρ¯, the mean of ϕ must be set to 1 ϕ¯ (τ)= D−1(τ,τ )[D(τ,τ )Φ(τ )D(τ,τ )] , (5) x −2 xy 0 0 0 0 yy foreverytimeτ. Notethatweintegrateovertheindexyinthisexpression,i.e. thediagonalofthecompositeoperator in square brackets is treated as a field. For a local growth operator, D (τ,τ )=D(τ,τ )δ (x y), this mean correction simplifies to xy 0 0 D − 1 ϕ¯ (τ)= D(τ,τ )Φˆ (τ ), (6) x 0 x 0 −2 where we defined Φˆ (τ ) Φ (τ ). x 0 xx 0 ≡ The Gaussian field ϕ and the growth operator D can be related to known quantities. To see this, consider the expansion of the Fourier modes of δ in Eulerian perturbation theory (see e.g. [45]), (cid:88) δ(k,τ)= D(n)(τ)δ(n)(k), (7) n=1 whereδ(n) areconvolutionsofninitialfieldsδ(k,τ )withanintegrationkernelthatchangesfromordertoorder. The 0 first term in this series is D(1)(τ)/D(1)(τ )δ (k,τ ), where D(1)(τ) is the growing solution to the linearized growth 0 0 0 equation [46]. We use this analogy to motivate the simplest possible form of the growth operation in the lognormal model and write D (τ)=δ (x y)D(1)(τ) (8) xy D − where we have set D(1)(τ ) = 1. This approximation erases any a-priori assumption of scale-dependent growth and 0 mode-coupling of the log field ϕ. Such a simplification is viable since the model in Eq. (8) describes only our prior assumptions about the density field ρ. The algorithm will find the most probable realization of ϕ 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 ϕ. 4 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 (τ) generates a four-dimensional field: three-dimensional xy spatial comoving volumes for every time-slice τ. 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 δx(cid:48) = x(cid:48)(x,τ)exp[Dxy(τ)ϕy(τ)] 1. (9) C − Applying a complicated, non-diagonal growth operation on all time slices separately before constructing the light cone is numerically not feasable for a high resolution reconstruction. The simplest form of D, Eq. (8), that depends only on time can be applied to the light cone directly δx(cid:48) =exp[Dx(cid:48)y(cid:48)ϕy(cid:48)] 1. (10) − For this work, we use the growth operator in Eq. (8) with the usual linear growth factor and leave extension of this model for future work. IV. DATA MODEL In this section we establish the analytic relation between the signal field ϕ, the field of overdensities δ 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 (cid:15)=(cid:15) +i(cid:15) such that (cid:15)=(a b)/(a+b)exp(2iη) for an ellipse with major axis a, minor axis b, and position angle η. 1 2 We use the common approxim−ation that the components of the intrinsic source galaxy ellipticity (cid:15)s =((cid:15)s,(cid:15)s), which 1 2 define the shape that would be observed in the absence of lensing1, follow a global bivariate Gaussian distribution with zero mean and variance σ2 per component: (cid:15) (cid:15)s (cid:45) ((cid:15)s 0,Ns), Ns =δ σ2. (11) ← N | ij ij (cid:15) This approximation has shortcomings (see e.g. [47]), 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 [40–42]. Lensing distorts the galaxy images in shape and size. If the distortion is small, i.e. in the limit of weak lensing, the relation between intrinsic source ellipticity and observed ellipticity can be linearised and simplifies to (cid:15)=g+(cid:15)s, (12) where g is the reduced shear. The reduced shear combines the effect of anisotropic lensing distortions, encoded in the shear γ =γ +iγ , and the isotropic distortion, encoded in the convergence κ 1 2 γ g = γ. (13) 1 κ ≈ − If κ 1, which is often the case for galaxy lensing, the reduced shear can be approximated by the shear itself g γ. (cid:28) ≈ The shear and convergence at angular position θ are related to the lensing potential ψ by γ (θ)= 1(cid:0)∂2 ∂2(cid:1)ψ(θ); γ (θ)=∂ ∂ ψ(θ); κ(θ)= 1(∂2+∂2)ψ(θ). (14) 1 2 1 − 2 2 1 2 2 1 2 The lensing potential is a weighted projection of the peculiar Newtonian gravitational potential φ along the line of sight. For a source at LOS distance ri, ri 2 (cid:90) ri r ψ(θ)= dr − φ(r,rθ ,rθ ) (15) c2 rri 1 2 0 1 Weusevectornotation,e.g. (cid:15),todenotethetupleofrealandimaginarypart((cid:15)1,(cid:15)2)ofacomplexnumber(cid:15)=(cid:15)1+i(cid:15)2. 5 in a spatially flat Universe. Applying the angular derivatives in Eq. (14) to the expression for the lensing potential in Eq. (15), we get ri (cid:90) ∂ ∂ ψ(θ)= 2 drW(r;ri)∂ ∂ φ(cid:0)r,rθ ,rθ (cid:1), (16) k l c2 rk rl 1 2 0 where k,l (1,2), and the lensing efficiency ∈ r(ri r) W(r;ri)= − . (17) ri Inpracticethedistancetothesourceri cannotbedetermineddirectlybutfollowsfromthephotometricallymeasured redshift zi. Photometrically measured redshifts are associated with a relatively high error, σ /(1+z) 0.03. In z ≈ its most simple form the algorithm ignores this uncertainty. We will use this simplified model to validate the func- tionality 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 ∂ ∂ φ(r) along the LOS of the ith source galaxy is obtained by rotating the tidal tensor on rk rl the global coordinate grid x(cid:48) Tij(x(cid:48))=∂x(cid:48)∂x(cid:48)φ(x(cid:48)), i,j (0,1,2), (18) i j ∈ intothespecificcoordinatesystem(withcoordinatesri)thatpointsintothedirectionofthisithgalaxyandprojecting it onto the (ri ri)-plane perpendicular to the LOS. 1− 2 The last relation required to connect the data to the density fluctuations is Poisson’s equation. It relates the potential φ(x(cid:48)) to the density fluctuations δ(x(cid:48)), 3 δ(x(cid:48)) 2φ(x(cid:48))= Ω H2 , (19) ∇ 2 m 0a(x(cid:48) /c) | | where H denotes the Hubble constant (which will be parametrized by h = H /(100kms−1Mpc−1) in our test 0 0 simulations), and a(x(cid:48) /c)=a(τ) denotes the scale factor at the time τ corresponding to LOS distance ri(x(cid:48)). | | V. IMPLEMENTATION The implementation, not only of the data model, but of the entire algorithm, is based on NIFTy [48], 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 [49]. 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 ϕ with a data vector of N s measuredsourceellipticities,areasfollows: ThegrowthoperatorD (τ)imprintsagrowthstructureontheGaussian xy field. The resulting four-dimensional field is plugged into the exponential of E()=exp() 1 [see Eq. (4)] to obtain · · − thefractionaloverdensityδ (τ). Theoverdensityinducesthepotentialφ (τ)bythePoissonequation, encodedinthe x y operator (τ)=∆−13Ω H2δ(x(cid:48))/a(x(cid:48) /c) [Eq. (19)]. The potential can be computed efficiently in Fourier space. Pyx 2 m 0 | | The gravitational potential is then restricted to the light cone of the observer by the light cone operator y(cid:48)(y,τ) C [Eq. (1)]. We compute the tidal tensor of the resulting 3D field [(Eq. (18)] by application of a global differential operator, which we denote z(cid:48)y(cid:48). The resulting tidal tensor is then rotated into each galaxy’s LOS coordinate system by a rotation operator, iT. An integration operator , which applies the integration in Eq. (16), integrates the Rrz(cid:48) Ijr componentsofeachoftheresultingN tidaltensorsalongtheunperturbedphotongeodesic. Forthisoperator,weuse s an adapted version of the implementation that was already successfully used in a similar reconstruction method [30]. The application of yields derivatives of the lensing potential for each galaxy location. From this we can compute jr I the shear components by a linear operator that comprises the equations in Eq. (14). Rotation and integration ij L 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 (x,τ) and work on the three-dimensional light cone from the beginning. In the prior, we model the Gaussian log-density ϕ with the power spectrum of matter fluctuations today P (k,a=1), where a denotes the scale factor. The growth operator is lin 6 diagonal in configuration space and only a function of comoving distance to the observer Dx(cid:48)y(cid:48) =D(1)(τ)δD(x(cid:48) y(cid:48)), − where D(1)(τ) is the growing solution of the linearized growth equation. The Poisson operator is split into two parts. First, a multiplication with 3/2Ω H2/a(x(cid:48) /c), i.e. an operation that is diagonal in configuration space. Second, the inverse Laplace operation ∆−1m, wh0ich|is d|iagonal 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 a2) if we apply the first part of the Poisson operator first and if D(1)(τ) is approximately proportional to the scale factor. In this case, the first order term of the exponential expression exp(Dx(cid:48)y(cid:48)ϕy(cid:48)) 1=Dx(cid:48)y(cid:48)ϕy(cid:48) + (ϕ2) a(x(cid:48) )δD(x(cid:48) y(cid:48))ϕy(cid:48) + (ϕ2), (20) − O ≈ | | − O andthereforethefirstordertimedependenceoftheoverdensityfieldwillpartlybecanceledbythe1/ainthePoisson 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 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 (Gamma-Projection-Rotation). G For one data point the implemented data model in operator notation reads di =(cid:15)i =(R(ϕ))i+(cid:15)si +ni =GijIjr(cid:48)Tr(cid:48)z(cid:48)Pz(cid:48)y(cid:48)[exp(Dy(cid:48)x(cid:48)ϕx(cid:48))−1y(cid:48)]+(cid:15)si +ni. (21) The total data vector d has dimensions 2 N , i.e. two ellipticity components for each of the N source galaxies. s s × We use the letter R to encode the total response of the data to a signal ϕ, i.e. the composite action of all operators. We have also added an experimental noise n here for completeness. In general this will be subdominant to the shape noise (cid:15)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 ϕ. The posterior distribution is related to the likelihood and the prior by Bayes’ theorem (dϕ) (ϕ) (ϕd)= P | P . (22) P | (d) P The prior probability (ϕ) is modeled as a Gaussian distribution with covariance Φ. In a first simple approximation P Φ can be taken to be diagonal in Fourier space with the usual power spectrum of matter overdensities P (k,a). lin To obtain the likelihood, we marginalize over all realizations of the shape noise (cid:90) (dϕ)= (cid:15)s (dϕ,(cid:15)s) ((cid:15)s)= [dR(ϕ),N], (23) P | D P | P N | 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 1 1 ln (ϕd)= [d R(ϕ)]†N−1[d R(ϕ)]+ (ϕ ϕ¯)†Φ−1(ϕ ϕ¯), (24) (cid:98) − P | 2 − − 2 − − where we have dropped most terms that are independent of the field of interest, ϕ. 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 ϕ. To find the minimum of thenegativelog-posterior we apply aNewton-like minimization scheme [50]. This requiresthe derivative of the negative log-posterior with respect to ϕ δlnP(ϕd) − δϕ | =Φ−uq1(ϕ−ϕ¯)q+[d−R(ϕ)]iNi−j1GjkIkr(cid:48)Tr(cid:48)z(cid:48)Pz(cid:48)y(cid:48)[exp(Dy(cid:48)x(cid:48)ϕx(cid:48))∗Dy(cid:48)u], (25) u where the star denotes a pointwise product in position space, i.e. (α β) =α β . x x x ∗ 7 VII. REDSHIFT-MARGINALIZED LIKELIHOOD AND POSTERIOR We can take into account the source redshift uncertainty by extending the marginalized likelihood in Eq. (23) by the redshift distribution function (zs) P (cid:90) (cid:90) (cid:90) (dϕ)= n zs (dϕ,n,zs) (n) (zs)= zs [dRzs(ϕ),N] (zs), (26) P | D D P | P P D N | P wherezs denotesthevectorofredshiftsofallsources. Inmostcases, thisintegrationcannotbedoneanalyticallyand theresultingdistributionwillingeneralnotbeGaussian. However,togoodapproximationwecanfitthemarginalized likelihood with a Gaussian that we characterize by the first and second central moment of the full distribution. To obtain them, we need to calculate d and dd† . For the expectation value of the data, we obtain (P(d|ϕ)) (P(d|ϕ)) (cid:104) (cid:105) (cid:104) (cid:105) (cid:90) (cid:90) (cid:90) d = zs dd [dRzs(ϕ),N] (zs)= zsRzs(ϕ) (zs) (27) P(d|ϕ) (cid:104) (cid:105) D D N | P D P  ∞ ∞  (cid:20)(cid:90) (cid:21) (cid:90) (cid:90) dzs (cid:104)di(cid:105)P(d|ϕ) =Gij dzsP(zs)Ijzrs R¯(ϕ)r =Gij dr drjW(r,rj)P(zjs)drjjR¯(ϕ)r (28) 0 r (cid:2) ¯R¯(cid:3)(ϕ) R˜(ϕ) , i i ≡ GI ≡ whereI¯denotestheredshiftaveragedintegrationoperatordefinedinthesquarebracketsinthefirstlineofEq.(28)and we have introduced R¯() (exp[D()] 1) to summarize the action of all source-redshift independent operators. · ≡TP · − The second moment is dd† =N + (R(ϕ))(R(ϕ))† . (29) (P(d|ϕ)) P(z) (cid:104) (cid:105) (cid:104) (cid:105) Non-diagonal elements of the second term in Eq. (29) read ((R(ϕ)) (R(ϕ)) ) = (R(ϕ)) (R(ϕ)) =(R˜ϕ) (R˜ϕ) (30) (cid:104) i j (cid:105)P(z) (cid:104) i(cid:105)P(z)(cid:104) j(cid:105)P(z) i j and diagonal elements are   ∞ ∞ ∞ (cid:90) (cid:90) (cid:90) dzj (cid:104)(R(ϕ))i(R(ϕ))i(cid:105)P(z) =GijGij dr dr(cid:48) drjW(r,rj)W(r(cid:48),rj)P(zjs)drjR¯(ϕ)rR¯(ϕ)r(cid:48) (31) 0 0 max(r,r(cid:48)) ij ij˜jr,r(cid:48)R¯(ϕ)rR¯(ϕ)r(cid:48), ≡G G I wherethenewoperator ˜ denotesthesquaredaverageoftheintegrationoperator, i.e. thesquarebracketsinthefirst I line of Eq. (31). The Gaussian approximation to the likelihood is then [dR˜(ϕ),N˜], where N˜ =N +Q and Q is N | Q= (R(ϕ)) (R(ϕ)) (R˜ϕ) (R˜ϕ) = (R(ϕ)) (R(ϕ)) (R˜ϕ) (R˜ϕ) . (32) i j P(z) i P(z) j P(z) i i P(z) i i (cid:104) (cid:105) −(cid:104) (cid:105) (cid:104) (cid:105) (cid:104) (cid:105) − Thisexpressionisstillsignal-dependentandweapproximateitfurtherbyreplacingϕbyitsposteriormean ϕ . P(ϕ|d) (cid:104) (cid:105) Since this mean depends on N˜, the resulting set of equations must be solved iteratively. VIII. VALIDATION AND TESTS Tovalidatetheimplementationandassessthegoodnessofthetomographicreconstruction,weperformanumberof increasinglyrealistictests,whichwedenotetestsA,BandC.IntestsoftypeAweemployanevensourcedistribution overtheentirebox. Thesetestsaretheleastrealisticonesandservetovalidatethecorrectnessoftheimplementation. Test B uses a realistic source 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 box2 and resolve the underlying and reconstructed overdensity fields with 1283 pixels. Depending on the test, we allow the physical sizes 2 Thispositionisnotfixedbythealgorithm. Forareconstructionfromasurveythatcoversasignificantfractionofthesky,theobserver canbeplacedinthecenterofthebox,forexample. 8 TABLE I. Quantitative comparison of reconstruction methods for tests A2 and A3. We compare the mean square pixel-wise difference between the underlying and reconstructed fields (δ and δˆ, respectively), their minimal and maximal values, and the Pearson-correlation coefficient between the true density and the reconstructions. (cid:113) redshifts (cid:104)(δˆ−δ)2(cid:105)/σ2 max(δˆ), max(δ) min(δˆ), min(δ) (cid:104)δˆδ(cid:105)/(σ σ ) δ δ δˆ lognormal WF lognormal WF N-body lognormal WF N-body lognormal WF 1. 0.70 0.71 4.07 1.69 4.92 -0.64 -1.04 -0.79 0.56 0.54 0.25 0.67 0.72 7.30 2.05 8.49 -0.87 -5.02 -0.87 0.58 0.54 0.0 0.66 0.71 7.83 2.27 10.26 -1.00 -5.96 -0.89 0.59 0.55 light cone 0.71 0.83 33.03 3.56 70.20 -1.00 -12.03 -0.92 0.56 0.44 of the box to differ. The current pixel resolution is limited due to computation time and memory usage. A higher resolutionwillbeaccessibleafterparallelizationandadaptionofthecodefortheusageonahigh-performancecluster. 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 [51] from which we take sub-volumes of size [1h−1Gpc]3 and smooth the density fields to fit the desired resolution of 1283 pixels. For more realistic test cases we construct a light cone of size [500h−1Mpc,500h−1Mpc,4000h−1Mpc] by joining snapshots of a sub-volume of the Millennium Run [52]. This sub- volume measures 500h−1Mpc along each side and we shift it every time 500h−1Mpc in the z-direction have been constructed in order to avoid that a LOS hits the same structure repeatedly. In the resulting light cone we achieve a resolution of 3.9h−1Mpc in the x- and y-directions corresponding to mildly non-linear scales and reach a redshift of z = 2.2 in the z-direction. Since we use 1283 pixel to resolve the light cone box, the physical size of the pixels is longer along the z-axis, meaning that we get a poorer resolution in this direction. Both the Millennium-XXL and the Millennium simulation, use a flat ΛCDM cosmology with Ω = 0.25 and m h=0.73. A. Tests A: Simple geometry In this series of tests we distribute 500 000 sources evenly in the box and slightly beyond ( 100h−1Mpc). We ± placesourcesoutsideofthereconstructedvolumeinordertoincreasetheareainwhichwecanrecovertheunderlying densityfieldwithhighresolution. Thequalityofthereconstructiondecreaseswithincreasingdistancetotheobserver 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 ϕ [Eq.(3)]inacubicboxofsidelength1000h−1Mpc. Wethenapplytheimplementeddatamodelandusethealgorithm 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 but a Gaussian prior on the overdensity field δ. Comparisons between input and reconstructed fields in each case are shown in Fig. 1. Table I summarizes pixelwise quadratic differences between underlying and reconstructed fields, minimal and maximal values and the pixel-wise Pearson correlation between the true density and its reconstructions. In all these parameters and for all redshifts tested the lognormal prior yields better results 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 distributionfunctions(PDFs),whichweshowinFig.2. WhiletheWienerFilterPDFisalmostsymmetricinallcases, we can capture some 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 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. 9 1 1 − ln[1+δ] FIG.1. QualitativecomparisonofreconstructionmethodsinthetestseriesA2,forwhichwecreatemockdatafromindividual simulation snapshots, apply an even source distribution and add negligible shape noise. We show central slices through the 3-dimensional fields. The observer is located in the center of the (x−y)-plane at the origin of the z-axis. The redshift of the snapshotdecreasesfromtoptobottom,z=[1, 0.25,0.]. Inthetwoupperrowsweshowthe(x−y)planeatz=max(z)/2,in the last row we look at the (x−z)-plane at y=max(y)/2. The fields are from left to right: the underlying density field from thesimulation, itsreconstructionwithalognormalprioranditsreconstructionwithaGaussianprior(Wienerfilter). Weplot ln[1+δ]andmarkunphysicalnegativedensitiesintheWienerFilterreconstructioninred. Theresolutionofthereconstruction decreaseswithdistancetotheobserver. Thisisbecause1)thedensityoflinesofsightdecreases2)therearelesssourcesbehind thepointwewanttoreconstruct3)theinformationfromthesesourcesissuppressedbytheshapeoftheintegrationkernel. In all cases, the lognormal reconstruction is superior in capturing the highest values of the density field and avoids unphysically low density contrasts below -1. 10 106 106 106 δ δ δ N body N body N body 105 − 105 − 105 − δˆ δˆ δˆ WF WF WF 104 δˆ 104 δˆ 104 δˆ e ln e ln e ln nc 103 nc 103 nc 103 e e e ur ur ur c 102 c 102 c 102 c c c o o o 101 101 101 100 100 100 10−1 1 0 1 2 3 4 10−1 4 2 0 2 4 6 8 10−1 4 2 0 2 4 6 8 10 − δ − − δ − − δ (a)z=1 (b)z=0.25 (c)z=0 FIG.2. 1-pointPDFsofthesimulateddensitycontrasts(incyan)andtheirWienerfilterandlognormalreconstructions(red anddarkblue)intestsonmockdatacreatedfromdifferentsnapshotsofanN-bodysimulationwithanevensourcedistribution andnegligibleshapenoise(A2). ThePDFofthelognormalreconstructionisslightlymoreskewedthantheWienerfilterPDF, which is close to symmetric in all cases. Both the Wiener filter and the lognormal reconstruction are equally biased in the positionofthepeak,however,themeanvaluesofbothreconstructionsagreewellwiththemeanvalueoftheunderlyingdensity. The mean values of the density contrast are indicated by vertical lines in the same color as their corresponding distribution. We also find that the Wiener filter produces unphysical density contrasts below -1 in all snapshots. 3) A test similar to A2 but with a redshift dependent density field constructed from different snapshots. The physicalsizeoftheboxspans[500h−1Mpc]2 inthe(x y)-planeand4000h−1Mpcintothez-directioncorresponding to a maximum redshift of z = 2.2.3 Results of test A3−are shown in Fig. 3 and Fig. 4. The quantitative comparison betweenthelognormalandthecorrespondingWFrunislistedinthelastrowofTab.Ilabeled“lightcone”. Wefinda similarmatchbetweentheunderlyingandreconstructedfieldsasintestsA2andagainasuperiorityofthelognormal prior over the WFreconstruction. This superiority isstill visible if we smooth the results ofboth reconstructions and the underlying field with a Gaussian kernel with σ =8h−1Mpc (Fig. 3, second row). B. 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.15o. The sources are distributed according to a distribution function of form n(z)=zαexp( (z/z )β), (33) 0 − wherewefitα,β,andz tothepubliclyavailablesourcedistributionoftheCFHTLenSsurvey4. Theresultingsource 0 distributionfunctionisshownintheleftpanelofFig.5. Thefitdoesnotexactlymatchthesurveydistribution,butthe similarity is sufficient for this test. We then draw Ω[armcin2]ρ [gal/arcmin2] source positions from this distribution, s where Ω is the angular opening area of the cone in arcmin2 and ρ the source density of the survey, which we choose s to be ρ = 11gal/arcmin2 corresponding to the source density in CFHTLenS. A scatter plot of the resulting spatial s source distribution is shown in the right panel of Fig. 5. The reconstructions are depicted in Fig. 6. We find that the superiority of the lognormal model remains with a realistic background galaxy distribution. If we restrict our analysis to the region that is covered by the observational cone, we find that the distribution of the reconstructed field values closely follows the underlying density distribution if we apply the lognormal prior. The Gaussian prior fails to capture the skewness and produces negative densities (Fig. 7). 3 The light cone is constructed by merging planes from different snapshots along the z-direction of the box, but the algorithm assumes thatthedistanceandredshiftincreaseinradialdirectionfromtheobserver. Thisleadstoamismatchbetweenthedistancetoa(source) positionassumedbythealgorithmandthedistanceorratherredshiftthatthispositioncorrespondstointhesimulation. Thedistance in the simulation is cosθ times the distance assumed by the algorithm, where θ is the angle between the LOS and the z-axis. This mismatch should only become relevant for relatively large angles θ>π/2. Lines of sight with such angles only cover a minor fraction oftheboxandwethereforeexpectthetesttonotbeseverelyaffectedbythisapproximation. 4 http://www.cfhtlens.org/astronomers/cosmological-data-products