# Lifting weak lensing degeneracies with a field-based likelihood

Natalia Porqueres, Alan Heavens, Daniel Mortlock, Guilhem Lavaux

► To cite this version:

Natalia Porqueres, Alan Heavens, Daniel Mortlock, Guilhem Lavaux. Lifting weak lensing degeneracies with a field-based likelihood. Monthly Notices of the Royal Astronomical Society, 2021, 509 (3), pp.3194-3202. 10.1093/mnras/stab3234. hal-03327445

HAL Id: hal-03327445

https://hal.science/hal-03327445

Submitted on 21 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Lifting weak lensing degeneracies with a field-based likelihood

Natalia Porqueres<sup>1,★</sup> Alan Heavens,<sup>1</sup> Daniel Mortlock<sup>1,2,3</sup> and Guilhem Lavaux<sup>1,4</sup>

Accepted 2021 November 2. Received 2021 October 30; in original form 2021 August 10

## ABSTRACT

We present a field-based approach to the analysis of cosmic shear data to infer jointly cosmological parameters and the dark matter distribution. This forward modelling approach samples the cosmological parameters and the initial matter fluctuations, using a physical gravity model to link the primordial fluctuations to the non-linear matter distribution. Cosmological parameters are sampled and updated consistently through the forward model, varying (1) the initial matter power spectrum, (2) the geometry through the distance-redshift relationship, and (3) the growth of structure and light-cone effects. Our approach extracts more information from the data than methods based on two-point statistics. We find that this field-based approach lifts the strong degeneracy between the cosmological matter density, $\Omega_m$, and the fluctuation amplitude, $\sigma_8$, providing tight constraints on these parameters from weak lensing data alone. In the simulated four-bin tomographic experiment we consider, the field-based likelihood yields marginal uncertainties on $\sigma_8$ and $\Omega_m$ that are, respectively, a factor of 3 and 5 smaller than those from a two-point power spectrum analysis applied to the same underlying data.

**Key words:** gravitational lensing: weak – methods: data analysis – large-scale structure of Universe.

## 1 INTRODUCTION

Weak gravitational lensing by large scale structure modifies the shape of galaxy images through a process known as cosmic shear. It is a powerful probe of cosmology as it is sensitive to the geometry of the Universe and the growth of cosmic structures, and is not dependent on a knowledge of galaxy bias. Cosmic shear data have provided useful constraints on the key cosmological parameters that describe the dark sector of the Universe (Troxel et al. 2018; Hikage et al. 2019; Hamana et al. 2020; Amon et al. 2021; Asgari et al. 2021; Secco et al. 2021), especially the amplitude of matter fluctuations. These cosmological constraints are usually inferred from the analysis of two-point statistics, either the lensing power spectrum (Hikage et al. 2019; Asgari et al. 2021) or the shear two-point correlation functions (Hildebrandt et al. 2017; Troxel et al. 2018; Hamana et al. 2020; Amon et al. 2021; Asgari et al. 2021; Secco et al. 2021). These standard methods based on the two-point summary statistics are sub-optimal for the non-Gaussian shear field and discard phase information. This has led to the development of several alternative inference techniques such as third- and higher-order statistics (Takada & Jain 2002; Pen et al. 2003; Jarvis, Bernstein & Jain 2004; Bernardeau 2005; Kilbinger & Schneider 2005; Semboloni et al. 2011; van Waerbeke et al. 2013; Fu et al. 2014; Jung et al. 2021), lensing peak count statistics (Jain & van Waerbeke 2000; Dietrich & Hartlap 2010; Maturi, Fedeli & Moscardini 2011; Petri et al. 2013; Lin & Kilbinger 2015; Liu et al.

2015; Kacprzak et al. 2016; Peel et al. 2017; Fluri et al. 2018b; Martinet et al. 2018; Shan et al. 2018; Harnois-Déraps et al. 2021) the lensing probability distribution function (Boyle et al. 2021; Martinet et al. 2021), shear clipping (Giblin et al. 2018), and machine learning approaches (Fluri et al. 2018a; Gupta et al. 2018; Ribli et al. 2019; Jeffrey, Alsing & Lanusse 2021). All these methods have the potential to improve the cosmology constraints, but they still rely on summary statistics that require an assumption of their sampling distribution, usually assumed to be Gaussian, which in turn requires a covariance matrix, which is difficult to compute reliably. One way to avoid the use of summary statistics and the problems associated with them is through data assimilation methods, in which the observations are incorporated into a forward model. Some forward modelling approaches have been developed for cosmic shear analysis (Alsing et al. 2016; Alsing, Heavens & Jaffe 2017; Böhm et al. 2017; Fiedorowicz et al. 2021; Porqueres et al. 2021). These methods differ in the quantities they sample and their prior assumptions. We have developed one such method by extending the Bayesian Origin Reconstruction from Galaxies (BORG, Jasche & Wandelt 2013; Lavaux, Jasche & Leclercq 2019) formalism to incorporate weak lensing as BORG-WL (Porqueres et al. 2021). BORG-WL differs from other forward modelling approaches in using a physical model of structure formation to link the primordial fluctuations to the evolved matter distribution, allowing us to sample the initial conditions, which follow simple Gaussian statistics. This data model, more complex and more complete, provides a better link between theory and data than using Gaussian (lognormal) priors for the shear (density) field.

<footnote n="1"><sup>1</sup>Imperial Centre for Inference and Cosmology (ICIC) & Astrophysics group, Imperial College, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK</footnote>

<footnote n="2"><sup>2</sup>Department of Mathematics, Imperial College London, London, SW7 2AZ, UK</footnote>

<footnote n="3"><sup>3</sup>Department of Astronomy, The Oskar Klein Centre, Stockholm University, Albanova, SE-10691 Stockholm, Sweden</footnote>

<footnote n="4"><sup>4</sup>Institut d'Astrophysique de Paris, CNRS & Sorbonne Université, UMR7095, F-75014, Paris, France</footnote>

<footnote>* E-mail: n.porqueres@imperial.ac.uk</footnote>
Here, we present an extended version of BORG-WL in which the cosmological parameters are sampled and updated consistently in the forward model. We vary the primordial matter power spectrum, the geometry of the universe (through the distance-redshift relation) and the growth rate of structures (through the dynamical model and its light-cone). We compared the results of BORG-WL to a standard analysis based on lensing power spectra applied to the same simulated data.

This paper is organised as follows. In Section 2, we describe the data model. In Section 3, we present the inference method used in this work, BORG-WL. Section 4 describes the simulated data employed to test and validate the method. The results are presented in Section 5, which includes the cosmology constraints from BORG-WL and a comparison to the standard analysis based on lensing power spectra. We also discuss the efficiency and convergence of the method and present the results of the validation tests. Section 6 summarises the results.

## 2 DATA MODEL

The effect of weak gravitational lensing on a source can be described by the shear, $\gamma$, which indicates the distortion in the shape of the image, and the convergence field, $\kappa$, which describes the variation in angular size. In the flat-sky approximation, which we assume throughout, the Fourier transforms of these two fields, $\tilde{\gamma}$ and $\tilde{\kappa}$, are related by

$$ \tilde{\gamma}(\ell) = \frac{(\ell_1 + i\ell_2)^2}{\ell^2} \tilde{\kappa}(\ell), \qquad (1) $$

where $\ell = (\ell_1, \ell_2)$ is the wavevector and the shear is written as a complex quantity. The convergence field is connected to the dark matter distribution by integrating the fractional overdensity along the line-of-sight using the lensing weights to give (Kilbinger 2015)

$$ \kappa(\mathbf{v}) = \frac{3H_0^2\Omega_m}{2c^2} \int_0^{r_{\text{lim}}} \frac{rdr}{a(r)} q(r)\delta^f(r\mathbf{v}, r), \qquad (2) $$

where $\mathbf{v}$ is the coordinate on the sky, $r$ is the comoving distance, $r_{\text{lim}}$ is the limiting comoving distance of the galaxy sample, $\delta^f$ is the dark matter overdensity at a scale factor $a$ and

$$ q(r) = \int_{r}^{r_{\text{lim}}} dr' n(r') \frac{r' - r}{r'}, \qquad (3) $$

with $n(r)$ being the redshift distribution of galaxy sources. We assume a spatially flat universe throughout.

These fields can be seen as latent parameters of a Bayesian hierarchical model as represented in Fig. 1. The starting point is the cosmological parameters $\theta$. Given a cosmology, we generate primordial fluctuations, $\delta^{\text{ic}}$, as Gaussian random fields in a rectangular cuboid of voxels, with a covariance matrix corresponding to the initial matter power spectrum. To compute the initial power spectrum, we used the prescription of Eisenstein & Hu (1998), Eisenstein & Hu (1999), which includes baryonic acoustic oscillations. The primordial fluctuations are linked to the evolved large-scale structures by the non-linear gravity model $\delta^f = G(\delta^{\text{ic}})$, which describes the growth of cosmic structures and accounts for light-cone effects. We then use the evolved density field, $\delta^f$, to generate the convergence field $\kappa$. In our discrete implementation and using the Born approximation, the radial line-of-sight integral in equation (2) is approximated by a sum over voxels as

$$ \kappa_{mn}^{b} = \frac{3H_0^2\Omega_m}{2c^2} \sum_{j=0}^{N} \delta_{mnj}^{f} \left[ \sum_{s=j}^{N} \frac{(r_s - r_j) n^b(r_s) \Delta r_s}{r_s} \right] \frac{r_j \Delta r_j}{a_j}, \qquad (4) $$

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Representation of the Bayesian hierarchical model. Blue boxes indicate deterministic functions and green boxes represent probability distributions. $\theta$ are the cosmological parameters, $g$ is a Gaussian field with zero mean and unit variance, $M(g, \theta)$ is the forward model, which consists of a convolution of the Gaussian field with the initial matter power spectrum $P(k)$ to obtain the initial conditions $\delta^{\text{ic}}$, a gravity and structure formation model $G(\delta^{\text{ic}})$ and the data model to obtain the shear $\gamma^b$ for each tomographic bin $b$. $P(d|\gamma_1^b, \gamma_2^b)$ is the likelihood with $d$ being the data.</figcaption></figure>

where the index $b$ indicates the tomographic bin, and the sub-indices $m$ and $n$ label the pixel on the sky, which is chosen to be large enough to contain many sources. The index $j$ labels the voxels along the line-of-sight at a comoving distance $r_j$. $N$ is the total number of voxels along the line of sight, which are found using a ray tracer. $\Delta r_j$ is the length of the line-of-sight segment inside the voxel $j$, and $\delta_f$ is the three-dimensional dark matter distribution. The comoving radial distance $r_s$ indicates the distance to the source. The redshift distribution of sources is given by $n^b(z_s)$ for each tomographic bin. Having evaluated $\kappa_{mn}^b$ in this way, we then use equation (1) to obtain the predicted shear field.

In this work, we do not include the important effects of intrinsic alignments or baryon feedback in our forward model, and we leave these to future work, where we will also include the uncertainty in the redshift distributions of the sources. These effects can be included in the inference as associated nuisance parameters.

## 3 METHOD

In this work, we extend the Bayesian hierarchical model presented in Porqueres et al. (2021), BORG-WL, to infer the cosmological parameters jointly with the density field. Here, we briefly describe the method and indicate the changes we have implemented to simultaneously constrain the cosmology and density field. A more detailed description of the general BORG framework can be found in Jasche & Kitaura (2010), Jasche & Wandelt (2013), Lavaux et al. (2019).

The BORG-WL method is a Bayesian framework to infer jointly the cosmology and the linear and non-linear matter distributions from pixelated tomographic shear fields. The underlying idea is to infer the initial conditions and the cosmological parameters from the data, using a full gravitational structure formation model. For that, the BORG framework employs a non-linear gravity model for structure
growth, either based on perturbation theory or N-body simulations. This dynamical model connects the primordial fluctuations to the evolved large-scale structures, allowing us to sample the primordial fluctuations, which follow simple Gaussian statistics. The BORG framework, therefore, uses a Gaussian prior for the initial conditions at an initial scale factor $a \approx 10^{-3}$. This Gaussian prior has zero mean and a covariance matrix corresponding to the initial matter power spectrum. In this extension of BORG-WL for cosmology inference, we vary the initial matter power spectrum according to the value of the cosmological parameters given by the sample, using the transfer function of Eisenstein & Hu (1998), Eisenstein & Hu (1999).

We have implemented a ray tracer to identify the voxels corresponding to each line of sight. With this ray tracer, we drop the distant observer approximation from our previous work and integrate radially from an observation point using the Born approximation. The line-of-sight integration depends on the cosmology through the distance-redshift relation, which we vary when sampling the cosmological parameters. The cosmology sampler, therefore, varies the geometry of the universe consistently with the cosmology.

Another element of the model that is sensitive to cosmology is the structure formation prescription. The dynamics and light-cone effects also vary consistently with the cosmological parameters. BORG-WL, therefore, samples and updates the cosmology consistently throughout the forward model by changing all the quantities that depend on the cosmological parameters. This consistent cosmology sampler is a powerful new addition to the BORG framework. Previously, Ramanah et al. (2019) presented a proof-of-concept to sample cosmological parameters with BORG from galaxy clustering data. However, their cosmology sampler did not update the forward model consistently, and the initial power spectrum was fixed. In contrast, our model is consistent, and its sensitivity to the cosmology is three-fold: Through the geometry of the problem, the growth of structures, and the initial matter power spectrum. This allows us to provide joint constraints on the cosmological parameters and the density field.

Our method proceeds as follows. Given values for the cosmological parameters, we generate a random realisation of the three-dimensional primordial fluctuations as a Gaussian field with a covariance matrix corresponding to the initial matter power spectrum. The structure formation model then evolves these initial conditions to a non-linear realisation of the dark matter distribution, correctly accounting for light-cone effects. By integrating the 3D dark matter field along the lines of sight and applying the lensing data model described in Section 2, BORG-WL predicts the two components of the pixelated tomographic shear fields on the flat sky, $\gamma_{1,mn}^b$, and $\gamma_{2,mn}^b$, where b denotes the tomographic bin and m and n are the sky pixel indices.

The measured shear in each pixel and bin, $\hat{\gamma}_{1,mn}^b$ and $\hat{\gamma}_{2,mn}^b$, will differ from the predicted values due to a combination of the galaxy shape measurement errors and range of intrinsic source shapes of the galaxies in the voxel. This is encoded in a likelihood of the form $P(\hat{\gamma}_{1,mn}^b, \hat{\gamma}_{2,mn}^b | \gamma_{1,mn}^b, \gamma_{2,mn}^b)$. Our numerical implementation does not rely on this likelihood having any specific mathematical form, and in particular can handle shear noise that varies with both sky pixel and tomographic bin. For the simulations used here, we assume that the observations are characterised by shape noise with variance $\sigma_e^2$, and that the associated shear uncertainty is hence given by $\sigma_b = \sigma_e / \sqrt{N_b}$, where $N_b$ is number of sources per pixel in tomographic bin b (assumed here to be independent of sky pixel). For sufficiently high $N_b$ the voxel likelihood can be approximated as being Gaussian with variance $\sigma_b^2$, which leads to a full log-likelihood

of the form

$$ \begin{aligned} \log \mathcal{L} &= \sum_b \sum_{mn} \log [P(\hat{\gamma}_{1,mn}^b, \hat{\gamma}_{2,mn}^b | \gamma_{1,mn}^b, \gamma_{2,mn}^b)] \\ &= -\sum_b \sum_{mn} \frac{(\hat{\gamma}_{1,mn}^b - \gamma_{1,mn}^b)^2 + (\hat{\gamma}_{2,mn}^b - \gamma_{2,mn}^b)^2}{2\sigma_b^2}, \end{aligned} \quad (5) $$

where normalising constants have been ignored and the dependence of the predicted shears on the underlying parameters is left implicit. This simple Gaussian likelihood is sufficient in this first demonstration of this method, but a more realistic form could be used when analysing a real dataset, or the hierarchical model could be extended to sample the shear as well (Schneider et al. 2015).

## 3.1 Sampling scheme

Sampling the density field implies that the amplitude of the density fluctuations in each voxel is a target parameter of the problem. This results in a high-dimensional problem, with typically $10^7$ parameters, that requires advanced statistical techniques to sample from the posterior distribution. We have split the problem in two sampling blocks following a Gibbs sampling scheme, with the density and the cosmological parameters sampled alternately according to

$$ \delta^{\text{ic}} \leftarrow P(\delta^{\text{ic}} | \theta, d), \quad (6) $$

$$ \theta \leftarrow P(\theta | \delta^{\text{ic}}, d). \quad (7) $$

To deal with the high number of parameters in the density field, BORG employs a Hamiltonian Monte Carlo sampler (Neal 2011), which uses the information in the gradients to navigate the parameter space. To sample the cosmology, we use a slice sampler. In both steps of the Gibbs sampling, the cosmology is consistent throughout the forward model.

In a follow-up work, we intend to include systematics modelling in our forward model and add the associated nuisance parameters to the inference. Sampling these few additional parameters will have a minimal extra cost and allow us to propagate their uncertainties automatically.

## 3.2 Setup for this work

In this work, we used Lagrangian perturbation theory (LPT) as our model of gravitational clustering and structure formation. At the resolution we use ($14 \ h^{-1} \ \text{Mpc}$), the LPT model does not show significant deviations from particle-mesh simulations (Tassev, Zaldarriaga & Eisenstein 2013). For future applications of this method at high resolution, we will use a fully non-linear particle-mesh such as the one described in Jasche & Lavaux (2019) to have an accurate description of the matter density at smaller scales.

Here, we focus on constraining three cosmological parameters $\Omega_m$, $h$, and $\sigma_8$, jointly with the density. We consider a flat universe, so $\Omega_\Lambda = 1 - \Omega_m$ can be treated as a derived parameter. We use the following uniform priors for the cosmological parameters: $\Omega_m$; $h$; and $\sigma_8$. To improve the efficiency of the sampler, we perform a coordinate transformation and sample $\theta' = (\Omega_m, h, S_8 \equiv \sigma_8 (\Omega_m/0.3)^{0.5})$ instead of $\theta = (\Omega_m, h, \sigma_8)$. This coordinate transformation reduces the correlation length of the sampler by a factor of 4, making the inference more efficient.
<figure><img src="image_4.png" /><figcaption>Figure 2. Redshift distributions of sources for each tomographic bin in this analysis.</figcaption></figure>

## 4 MOCK DATA

To test our method, we generated mock data using the forward model illustrated in Fig. 1.

To generate the mock data, we assumed a standard $\Lambda$CDM cosmology with the following parameters: $\Omega_m = 0.3175$, $\Omega_\Lambda = 0.6825$, $\Omega_b = 0.049$, $h = 0.677$, $\sigma_8 = 0.8$, $n_s = 0.9624$ with $H_0 = 100 \text{ h km s}^{-1} \text{ Mpc}^{-1}$. Then, we generated Gaussian initial conditions on a Cartesian grid of size $(1 \times 1 \times 4) \text{ h}^{-1}$ Gpc, with $64 \times 64 \times 128$ voxels. We evolved the primordial fluctuations via LPT, accounting for light-cone effects. To obtain the evolved density fields, we used a cloud-in-cell scheme to estimate the density on the Cartesian grid from simulated particles. We generated tomographic shear fields by applying the lensing data model in Section 2 and using the redshift distribution of sources shown in Fig. 2. The angular scale of the voxels for the different tomographic bins is 40, 23, 17, and 15 arcmin. Finally, we added Gaussian pixel noise with a variance that corresponds to 30 galaxies per square arcmin, as expected for the Euclid survey (Laureijs et al. 2011), and with an uncertainty on intrinsic ellipticity $\sigma_\epsilon = 0.3$ (Kilbinger 2015). In this experiment, we set $\sigma_\epsilon$ as the variance of a single component of ellipticity. The sources are then distributed uniformly over the four tomographic bins. The resulting noise level for the farthest tomographic bin is $N = 0.013$, and the data are noise-dominated at all scales.

## 5 RESULTS

### 5.1 Cosmology constraints and comparison to standard analysis

Here, we present the mock-constraints on cosmology obtained by applying BORG-WL to the simulated data described in Section 4. We also compare our results to one of the common techniques for analysing shear data, i.e. an analysis of the lensing power spectra assuming a Gaussian likelihood with a fixed covariance matrix.

To compare the results of BORG-WL to the standard 2-point statistics approach, we focus on the angular power spectra $C_\ell^{ab}$, defined implicitly by

$$ \langle \kappa^a(\ell) \kappa^{b*}(\ell') \rangle = (2\pi)^2 \delta_D(\ell - \ell') C_\ell^{ab}, \quad (8) $$

where a and b label the tomographic bins. Following the common approach in lensing analyses, we assume that the sampling distribution

<figure><img src="image_5.png" /><figcaption>Figure 3. Comparison of $\Omega_m-\sigma_8$ constraints from our method BORG-WL (blue) and the standard analysis based on lensing tomographic power spectra (orange). The contours correspond to the 68.3% and 95.5% highest posterior density credible regions. Both posteriors are obtained by applying both methods to the same simulated shear data, with 4 tomographic bins and 30 galaxies per square arcmin. The dashed lines indicate the true values of the parameters.</figcaption></figure>

of the $C_\ell^{ab}$ is Gaussian,

$$ \log \mathcal{L} = -\frac{1}{2} [\mathbf{d} - \mathbf{m}(\boldsymbol{\theta})]^\mathrm{T} \mathbf{C}^{-1} [\mathbf{d} - \mathbf{m}(\boldsymbol{\theta})] + \text{constant}, \quad (9) $$

where the data vector $\mathbf{d}$ contains the auto- and cross-power spectra of the tomographic bins. $\mathbf{C}$ is the covariance matrix, which is assumed to be fixed and computed at a fiducial cosmology ($\Omega_m^f, \sigma_8^f, h^f) = (0.3, 0.8, 0.67)$. To compute the covariance matrix, we used the Gaussian approximation in Kilbinger (2015),

$$ \mathrm{Cov}(C_{\ell}, C'_{\ell}) = \frac{1}{f_{\mathrm{sky}} \ell \Delta \ell} (C_{\ell} + N_{\epsilon})^2 \delta_{\ell\ell'}, \quad (10) $$

where $f_{\mathrm{sky}}$ is the fraction of the sky observed by the survey, $\Delta\ell$ is the width of the $\ell$-bin, $\delta_{\ell\ell'}$ is a Kronecker delta, which makes the covariance diagonal, and $N_\epsilon$ is the noise power given by $\sigma_\epsilon^2/2n_g$, with $\sigma_\epsilon$ being the intrinsic dispersion of galaxy ellipticities and $n_g$ being the mean number density of source galaxies on the sky. For this analysis, we use broad uniform priors for $\Omega_m, \sigma_8$, and $h$ on the intervals $[0.05, 1.45]$, $[0.3, 1.45]$, and $[0.4, 0.9]$, sufficiently broad that the posterior is unaffected. To compute the model prediction $\mathbf{m}(\boldsymbol{\theta})$, we use the Core Cosmology Library (Chisari et al. 2019).

Fig. 3 shows the 68.3% and 95.5% credible regions for $\Omega_m$ and $\sigma_8$ obtained with BORG-WL. For comparison, we added the constraints from the standard power spectra analysis on the same underlying data, showing that BORG-WL lifts the weak lensing degeneracy and yields much tighter constraints on these parameters than the standard analysis based on two-point summary statistics. We note that the size of the posterior distribution is significantly smaller than the adopted size for the prior distributions for $\Omega_m$ and $\sigma_8$. We conclude that the constraints on these parameters should not be affected by the specific choice of prior. Fig. 3 is the main result of this paper.

Fig. 4 shows the full corner plot including $h$ and the marginal distributions. We computed the BORG-WL posteriors from 10 000 samples after the burn-in phase, which is determined by examining the trace plots of the parameters in Fig. 5. These results show that BORG-WL recovers the correct value of $\Omega_m$ and $\sigma_8$. Since the lensing data is not informative about the value of $h$ (Hall 2021), this parameter is not well constrained, and the marginal distribution corresponds
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Posterior distribution of the sampled cosmological parameters from the BORG-WL analysis (blue) and the standard analysis based on the lensing power spectra (orange). The dashed lines indicate the ground truth.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Trace plot for the cosmology sampler. The burn-in phase can be identified as the drift of the parameters in the first ~2500 samples, in particular $S_8$. The dashed lines indicate the true value of the parameters.</figcaption></figure>

closely to the prior. For comparison, we added the constraints from the standard analysis based on the lensing power spectra. The corresponding posterior distributions are evaluated in grids of 50 × 50 points. The comparison of the marginal distributions shows that BORG-WL provides far more precise cosmology constraints than the standard method, improving the error bars of $\Omega_m$ by a factor 5 and $\sigma_8$ by a factor 3 (see Table 1). The ratio of the constraints

on $S_8$ from $C_\ell$ versus BORG-WL is 1.2. The improvement on the constraints of $S_8$ is small compared to $\Omega_m$ and $\sigma_8$ because $S_8$ is a combination of cosmological parameters found to be optimally constrained by the standard analysis with power spectra.

By building this Bayesian field-based approach, we can extract information from the data beyond the two-point statistics and lift the $\Omega_m-\sigma_8$ degeneracy that results from using the standard approach in
<table><caption>**Table 1.** Mean of the posterior, and 68.3% credible intervals for $\Omega_m$ and $\sigma_8$ from tomographic power spectra $C_l^{ab}$ and BORG-WL.</caption><thead><tr><th></th><th>Power spectra</th><th>BORG-WL</th></tr></thead><tbody><tr><td>$\Omega_m$</td><td>$0.3^{+0.5}_{-0.2}$</td><td>$0.34^{+0.07}_{-0.05}$</td></tr><tr><td>$\sigma_8$</td><td>$0.6^{+0.4}_{-0.3}$</td><td>$0.79^{+0.11}_{-0.11}$</td></tr><tr><td>$S_8$</td><td>$0.80^{+0.15}_{-0.10}$</td><td>$0.83^{+0.10}_{-0.10}$</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Auto-correlation of the cosmological parameters as a function of the sample number in the Markov chain. The correlation length of the sampler can be estimated as the point when the correlation drops below 0.1.</figcaption></figure>

weak lensing analysis to date. This demonstrates that the high-order statistics, which are overlooked by the standard methods, contain significant cosmological information.

The spectacular improvement is in line with the results of Leclercq & Heavens (2021), who showed for a simplified lognormal model that data assimilation approaches can lift degeneracies that are present in analyses based on two-point summary statistics, even for nearly Gaussian fields.

## 5.2 Efficiency of the cosmology sampler

By construction, subsequent samples in Markov chains are correlated. The statistical efficiency of an MCMC is determined by the number of independent samples that can be drawn from a chain of a certain length. To estimate the efficiency of the sampler, we measure the correlation length of the cosmological parameters along the Markov chain. We typically define the correlation length as the lag at which the auto-correlation drops below 0.1 (Gregory 2010) Fig. 6 presents the results of this test, showing that the correlation length for the parameters $\Omega_m$ and $\sigma_8$ are 600 and 200 samples, respectively. The parameter $h$ shows a very short correlation length because the data is not constraining.

## 5.3 Convergence tests

To assess whether the Markov chain has fully explored the target posterior distribution or more samples are required, we consider the convergence of the Markov chain. A common approach to assess the convergence is based on comparing multiple chains with different starting points. One way to compare the chains is through a Gelman–Rubin diagnostic (Gelman & Rubin 1992), where the convergence is assessed by comparing the intra-chains and inter-chains variances

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Contours of 68.3% and 95.5% highest posterior density credible regions for two chains with different starting points. As can be noted here, both chains converge qualitatively to the same distribution.</figcaption></figure>

for each model parameter. If the difference between these variances is large, the chains have not converged.

In this work, we run two chains with starting points ($\Omega_m, h, \sigma_8$) = (0.31, 0.677, 0.8) and ($\Omega_m, h, \sigma_8$) = (0.4, 0.7, 0.55), respectively. The initial conditions are also randomised. After 3000 samples (discarding the burn-in phase), the Gelman–Rubin diagnostic is $R < 1.03$ for all the cosmological parameters. This indicates that the chains have converged reasonably well. Fig. 7 presents the $\Omega_m-\sigma_8$ constraints from both chains, showing that they converge to the same results.

## 5.4 Inferred density fields

Jointly with the cosmology results presented in the previous sections, BORG-WL also infers the 3D matter density field. Since we have shown that the method correctly recovers the underlying cosmology, here, we focus on testing and validating the field results. To do so, we compare the inferred density fields to the truth and test whether the inferred cosmology and densities can explain the data. Although the tests presented here are the same as in our previous work Porqueres et al. (2021), there is a significant difference: Here, we sample the cosmology and the density field jointly, while in our previous work, we used a fixed cosmology. Therefore, it is important to validate this extension of the method by running a similar set of consistency tests to confirm that BORG-WL recovers the density fields when the cosmology is sampled.

Since BORG-WL uses MCMC to explore the parameter space, we have access to the full posterior distribution, and we can draw posterior samples of the primordial matter fluctuations and the matter distribution. We use these samples to compute the mean and uncertainty of the inferred fields. Fig. 8 shows the projection on the sky of the true density fields and the mean and variance of the posterior samples. Since the lensing data are not very informative in the radial direction, we projected the three-dimensional fields on the sky, weighting the planes of the density field with the distribution of sources. A visual comparison between the true and inferred fields shows that the method correctly recovers the structures from the truth. As expected, the estimated mean of the fields shows some degree of smoothing, a feature familiar from Wiener filtering (Gregory 2010). This smoothing is an effect of averaging several samples, but the individual posterior samples have the correct amplitude and power
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Projection of the initial (upper panels) and evolved (bottom panels) density fields. The first column shows the true fields that we used to generate the data, the second column shows the ensemble mean, and the third column shows the standard deviation of the fields. The mean and standard deviation are estimated from 500 MCMC samples.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Posterior initial matter power spectrum, measured from the posterior samples of the three-dimensional initial conditions. The input power spectrum is shown as the dashed curve.</figcaption></figure>

spectrum. Fig. 9 shows the matter power spectrum of the initial conditions. The power spectrum is directly measured from the three-dimensional samples and agrees with the truth, indicating that the posterior samples have the correct clustering properties.

To test whether the inferred cosmology and densities provide accurate explanations of the data, we do a posterior predictive test (Gelman et al. 2004). This consists of applying the data model to the

inferred quantities and checking whether we recover the structures in the data. Fig. 10 shows the results of this test for one tomographic bin, showing that the posterior-predicted shear field has the features of the true shear, and, therefore, the model can explain the data. At the resolution of this work, we are in the weakly non-Gaussian regime. We can expect a further improvement of the cosmology constraints when going to a higher resolution, but at smaller scales the baryon effects will be important. In a follow-up work, we will incorporate baryon feedback into the data model and increase the resolution of the analysis.

## 6 SUMMARY AND CONCLUSIONS

We have presented a field-based approach to infer jointly cosmological parameters and the dark matter distribution, extending our earlier BORG-WL implementation (Porqueres et al. 2021), by sampling cosmological parameters and upgrading the line-of-sight integrations. The forward model consists of uniform priors for the cosmological parameters, a Gaussian prior for the primordial fluctuations, and a physical gravity model for the non-linear dynamics that link the initial conditions and the evolved density field. The cosmological parameters are sampled and updated consistently throughout the forward model by varying the initial matter power spectrum, the geometry of the lines of sight (through the distance-redshift relation) and the growth of cosmic structures (through the dynamical model and light-cone). This allows us to constrain the cosmology and the matter distribution simultaneously.

We demonstrated the application of this method on simulated shear data with four tomographic bins and the number of sources per square-arcmin expected for *Euclid*. We sampled $\Omega_m$, $h$, and $\sigma_8$ and compared the constraints from our field-based approach to the standard analysis based on the power spectra, with the same simulated data. The results showed that BORG-WL lifts the $\Omega_m-\sigma_8$ degeneracy
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Posterior predicted maps for one tomographic bin. The first column shows the noisy shear mock data, the second column shows the true shear, and convergence fields, and the third and fourth columns are the mean and standard deviation computed from 500 posterior-predicted fields.</figcaption></figure>

and provides far more precise constraints from the same underlying data: The uncertainties in $\Omega_m$ and $\sigma_8$ are decreased enormously, by a factor of 5 and 3, respectively. This demonstrates the advantage of using data assimilation techniques and applying the likelihood to the shear fields themselves, rather than to summary two-point statistics.

In future work, we will extend our forward model to include systematic effects such as intrinsic alignments, baryon feedback, and uncertain redshift distributions, and apply this method to real data to constrain jointly cosmological parameters and the underlying matter distribution.

## ACKNOWLEDGEMENTS

This work was supported by STFC through Imperial College Astro-physics Consolidated grant ST/5000372/1. GL acknowledges financial support from the ANR BIG4, under reference ANR-16-CE23-0002. This work was carried out within the Aquila Consortium.<fnref n="1" />

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

## REFERENCES

* Alsing J., Heavens A., Jaffe A. H., Kiessling A., Wandelt B., Hoffmann T., 2016, MNRAS, 455, 4452

* Alsing J., Heavens A., Jaffe A. H., 2017, MNRAS, 466, 3272

* Amon A. et al., 2021, preprint (arXiv:2105.13543)

* Asgari M. et al., 2021, A&A, 645, A104

* Bernardeau F., 2005, A&A, 441, 873

* Böhm V., Hilbert S., Greiner M., Enßlin T. A., 2017, Phys. Rev. D, 96, 123510

* Boyle A., Uhlemann C., Friedrich O., Barthelemy A., Codis S., Bernardeau F., Giocoli C., Baldi M., 2021, MNRAS, 505, 2886

* Chisari N. E. et al., 2019, ApJS, 242, 2

* Dietrich J. P., Hartlap J., 2010, MNRAS, 402, 1049

* Eisenstein D. J., Hu W., 1998, ApJ, 496, 605

* Eisenstein D. J., Hu W., 1999, ApJ, 511, 5

* Fiedorowicz P., Rozo E., Boruah S. S., Chang C., Gatti M., 2021, preprint (arXiv:2105.14699)

* Fluri J., Kacprzak T., Refregier A., Amara A., Lucchi A., Hofmann T., 2018a, Phys. Rev. D, 98, 123518

* Fluri J., Kacprzak T., Sgier R., Refregier A., Amara A., 2018b, J. Cosmology Astropart. Phys., 2018, 051

* Fu L. et al., 2014, MNRAS, 441, 2725

* Gelman A., Rubin D. B., 1992, Stat. Sci., 7, 457

* Gelman A., Carlin J. B., Stern H. S., Rubin D. B., 2004, Bayesian Data Analysis, 2nd edn. Chapman and Hall/CRC, New York

* Giblin B. et al., 2018, MNRAS, 480, 5529

* Gregory P., 2010, Bayesian Logical Data Analysis for the Physical Sciences. Cambridge University Press, Cambridge

* Gupta A., Matilla J. M. Z., Hsu D., Haiman Z., 2018, Phys. Rev. D, 97, 103515

* Hall A., 2021, MNRAS, 505, 4935

* Hamana T. et al., 2020, PASJ, 72, 16

* Harnois-Dérap s J., Martinet N., Castro T., Dolag K., Giblin B., Heymans C., Hildebrandt H., Xia Q., 2021, MNRAS, 506, 1623

* Hikage C. et al., 2019, PASJ, 71, 43

* Hildebrandt H. et al., 2017, MNRAS, 465, 1454

* Jain B., van Waerbeke L., 2000, ApJ, 530, L1

* Jarvis M., Bernstein G., Jain B., 2004, MNRAS, 352, 338

* Jasche J., Kitaura F. S., 2010, MNRAS, 407, 29

* Jasche J., Lavaux G., 2019, A&A, 625, A64

* Jasche J., Wandelt B. D., 2013, MNRAS, 432, 894

* Jeffrey N., Alsing J., Lanusse F., 2021, MNRAS, 501, 954

<footnote n="1">¹https://aquila-consortium.org</footnote>
Jung G., Namikawa T., Liguori M., Munshi D., Heavens A., 2021, J. Cosmology Astropart. Phys., 2021, 055

Kacprzak T. et al., 2016, MNRAS, 463, 3653

Kilbinger M., 2015, Rep. Prog. Phys., 78, 086901

Kilbinger M., Schneider P., 2005, A&A, 442, 69

Laureijs R. et al., 2011, preprint (arXiv:1110.3193)

Lavaux G., Jasche J., Leclercq F., 2019, MNRAS, preprint (arXiv:1909.06396)

Leclercq F., Heavens A., 2021, MNRAS, 506, L85

Lin C.-A., Kilbinger M., 2015, A&A, 583, A70

Liu J., Petri A., Haiman Z., Hui L., Kratochvil J. M., May M., 2015, Phys. Rev. D, 91, 063507

Martinet N. et al., 2018, MNRAS, 474, 712

Martinet N., Harnois-Déraps J., Jullo E., Schneider P., 2021, A&A, 646, A62

Maturi M., Fedeli C., Moscardini L., 2011, MNRAS, 416, 2527

Neal R., 2011, MCMC Using Hamiltonian Dynamics. Chapman and Hall/CRC, New York, p. 113

Peel A., Lin C.-A., Lanusse F., Leonard A., Starck J.-L., Kilbinger M., 2017, A&A, 599, A79

Pen U.-L., Zhang T., van Waerbeke L., Mellier Y., Zhang P., Dubinski J., 2003, ApJ, 592, 664

Petri A., Haiman Z., Hui L., May M., Kratochvil J. M., 2013, Phys. Rev. D, 88, 123002

Porqueres N., Heavens A., Mortlock D., Lavaux G., 2021, MNRAS, 502, 3035

Ramanah D. K., Lavaux G., Jasche J., Wandelt B. D., 2019, A&A, 621, A69

Ribl D., Pataki B. Á., Zorrilla Matilla J. M., Hsu D., Haiman Z., Csabai I., 2019, MNRAS, 490, 1843

Schneider M. D., Hogg D. W., Marshall P. J., Dawson W. A., Meyers J., Bard D. J., Lang D., 2015, ApJ, 807, 87

Secco L. F. et al., 2021, preprint (arXiv:2105.13544)

Semboloni E., Schrabback T., van Waerbeke L., Vafaei S., Hartlap J., Hilbert S., 2011, MNRAS, 410, 143

Shan H. et al., 2018, MNRAS, 474, 1116

Takada M., Jain B., 2002, MNRAS, 337, 875

Tassev S., Zaldarriaga M., Eisenstein D. J., 2013, J. Cosmology Astropart. Phys., 2013, 036

Troxel M. A. et al., 2018, Phys. Rev. D, 98, 043528

van Waerbeke L. et al., 2013, MNRAS, 433, 3373

This paper has been typeset from a *TEX/LATEX* file prepared by the author.