# hybrid - lenstool: a self-consistent algorithm to model galaxy clusters with strong- and weak-lensing simultaneously

Anna Niemiec, Mathilde Jauzac, Eric Jullo, Marceau Limousin, Keren Sharon, Jean-Paul Kneib, Priyamvada Natarajan, Johan Richard

► To cite this version:

Anna Niemiec, Mathilde Jauzac, Eric Jullo, Marceau Limousin, Keren Sharon, et al. hybrid - lenstool : a self-consistent algorithm to model galaxy clusters with strong- and weak-lensing simultaneously. Monthly Notices of the Royal Astronomical Society, 2020, 493 (3), pp.3331-3340. 10.1093/mn-ras/staa473 . hal-03581360

HAL Id: hal-03581360

https://hal.science/hal-03581360

Submitted on 19 Feb 2022

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>

# hybrid-LENSTOOL: a self-consistent algorithm to model galaxy clusters with strong- and weak-lensing simultaneously

Anna Niemiec<sup>1,★</sup> Mathilde Jauzac<sup>1,2,3,4</sup> Eric Jullo,<sup>5</sup> Marceau Limousin,<sup>5</sup>

Keren Sharon<sup>1</sup>, Jean-Paul Kneib,<sup>6</sup> Priyamvada Natarajan<sup>7</sup> and Johan Richard<sup>8</sup>

Accepted 2020 February 11. Received 2020 February 11; in original form 2019 November 25

## ABSTRACT

We present a new galaxy cluster lens modelling approach, *hybrid-LENSTOOL*, that is implemented in the publicly available modelling software *LENSTOOL*. *hybrid-LENSTOOL* combines a parametric approach to model the core of the cluster, and a non-parametric (free-form) approach to model the outskirts. *hybrid-LENSTOOL* optimizes both strong- and weak-lensing constraints simultaneously (*Joint-Fit*), providing a self-consistent reconstruction of the cluster mass distribution on all scales. In order to demonstrate the capabilities of the new algorithm, we tested it on a simulated cluster. *hybrid-LENSTOOL* yields more accurate reconstructed mass distributions than the former *Sequential-Fit* approach where the parametric and the non-parametric models are optimized successively. Indeed, we show with the simulated cluster that the mass density profile reconstructed with a *Sequential-Fit* deviates from the input by $2-3\sigma$ at all scales while the *Joint-Fit* gives a profile that is within $1-1.5\sigma$ of the true value. This gain in accuracy is consequential for recovering mass distributions exploiting cluster lensing and therefore for all applications of clusters as cosmological probes. Finally we found that the *Joint-Fit* approach yields shallower slope of the inner density profile than the *Sequential-Fit* approach, thus revealing possible biases in previous lensing studies.

**Key words:** gravitational lensing: strong – gravitational lensing: weak – galaxies: clusters: general.

## 1 INTRODUCTION

Gravitational lensing is the bending of the light emitted by a background source as it grazes past the gravitational potential of a foreground object called the lens (for reviews see Kneib & Natarajan 2011; Hoekstra et al. 2013; Bartelmann & Maturi 2017). The lens can be any type of objects with masses ranging from a planet to a massive galaxy cluster. Here, we focus on the lensing of a background galaxy population by a massive, foreground cluster of galaxies. Gravitational lensing is observed in two regimes defined by the intensity of the distortions created by the gravitational potential of the lens: the strong-lensing regime where background galaxies are highly distorted into gravitational arcs and multiple images;

and the weak-lensing regime wherein the distortions of background galaxies are small and need to be treated statistically.

While gravitational lensing generated by a galaxy lens was first observed in 1979 by Walsh, Carswell & Weymann (1979) with the multiply imaged quasar Q0957+561A-B, astronomers had to wait until the late 1980s to confirm the first observation of a gravitational arc in a galaxy cluster (Lynds & Petrosian 1986; Soucail et al. 1987, 1988). Since then, gravitational lensing by galaxy clusters has emerged as a powerful tool to study the Universe. Indeed, gravitational lensing is a unique tool to map the mass distribution of the lenses as it is independent of their dynamical state, thereby providing crucial *in situ* information on the physics of these objects. In the case of cluster lenses, a lot of work has been done from gravitational lensing mass maps and multiwavelength analyses to constrain cluster physics (e.g. Natarajan, Kneib & Smail 2002; Kneib et al. 2003; Clowe, Gonzalez & Markevitch 2004; Bradač

<footnote n="1"><sup>1</sup>Department of Astronomy, University of Michigan, 1085 South University Ave, Ann Arbor, MI 48109, USA</footnote>

<footnote n="2"><sup>2</sup>Centre for Extragalactic Astronomy, Department of Physics, Durham University, Durham DH1 3LE, UK</footnote>

<footnote n="3"><sup>3</sup>Institute for Computational Cosmology, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="4"><sup>4</sup>Astrophysics and Cosmology Research Unit, School of Mathematical Sciences, University of KwaZulu-Natal, Durban 4041, South Africa</footnote>

<footnote n="5"><sup>5</sup>Aix Marseille Univ, CNRS, CNES, LAM, F-13388 Marseille, France</footnote>

<footnote n="6"><sup>6</sup>Laboratoire d'Astrophysique, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, CH-1290 Versoix, Switzerland</footnote>

<footnote n="7"><sup>7</sup>Department of Astronomy, Yale University, New Haven, CT 06511, USA</footnote>

<footnote n="8"><sup>8</sup>Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-Genis-Laval, France</footnote>

<footnote>* E-mail: annaniem@umich.edu</footnote>
et al. 2006; Merten et al. 2011; Jauzac et al. 2012, 2015; Diego et al. 2015; Eckert et al. 2015; Sharon et al. 2015, 2020; Mahler et al. 2018), and dark matter properties (e.g. Natarajan et al. 2002, 2017; Bradač et al. 2008; Harvey et al. 2015; Massey et al. 2015, 2018; Harvey, Kneib & Jauzac 2016; Jauzac et al. 2016a, 2018; Harvey et al. 2017, 2019). Moreover, lensing can give us hints on galaxy evolution (e.g. Natarajan et al. 1998; Limousin et al. 2007; Natarajan et al. 2009; Leauthaud et al. 2012, 2015; Sifón et al. 2015; Li et al. 2016; Niemiec et al. 2017), and on the distant Universe as lenses behave as cosmic telescopes and thus allow us to observe high-redshift galaxies (e.g. D’Aloisio, Natarajan & Shapiro 2014; Atek et al. 2015, 2018; Alavi et al. 2016; Bouwens et al. 2017; Ishigaki et al. 2018; Kawamata et al. 2018), study highly magnified galaxies at intermediate redshifts (e.g. Teplitz et al. 2000; Bayliss et al. 2014; Johnson et al. 2017; Rigby et al. 2018; Bayliss et al. 2020; Chisholm et al. 2019; Rivera-Thorsen et al. 2019; Sharon et al. 2020), and/or lensed transients (e.g. Kelly et al. 2015, 2018; Rodney et al. 2015, 2018; Diego et al. 2016, 2018, 2019; Jauzac et al. 2016b; Treu et al. 2016; Goobar et al. 2017; Smith et al. 2018, 2019). Gravitational lensing can even be used to constrain cosmological parameters, as it is sensitive to the geometry of the Universe itself (Jullo et al. 2010; D’Aloisio & Natarajan 2011; Caminha et al. 2016; Acebron et al. 2017; Suyu et al. 2018; Birrer et al. 2019; Wong et al. 2019). However, to exploit gravitational lensing fully, the mass distribution of the lenses themselves requires to be modelled with high precision and accuracy.

There are currently two classes of lens modelling algorithms. The first one is based on *parametric* mass models: the total mass distribution of the cluster is decomposed into a finite number of mass components divided into: (i) the large-scale components, representing the mass contribution of the cluster dark matter haloes and gas in the intra-cluster medium (ICM), and (ii) the galaxy-scale components, representing the mass contribution of cluster galaxies. Dark matter haloes are associated with individual cluster galaxies as smaller scale contributions to the overall mass distribution. The mass distribution of each component is then described by an analytical density profile, the most commonly employed ones being the Singular Isothermal Sphere potentials (SIS, see e.g. Binney & Tremaine 1987), the Navarro–Frenk–White potentials (NFW, Navarro, Frenk & White 1996), and the Pseudo Isothermal Elliptical Mass Distribution potentials (PIEMD, Elíasdóttir et al. 2007). Current parametric mass modelling algorithms include LENSTOOL (Jullo et al. 2007), GLAFIC (Oguri 2010), and LTM (Zitrin et al. 2012, 2013).

The second class of algorithms relies on the so-called *free-form* (or non-parametric) models. In this case, the cluster mass distribution is subdivided into a grid of mass 'pixels', and the amplitude (and possibly the shape) of the pixels are optimized so that the overall mass distribution reproduces best the observed lensed image constraints. Free-form reconstruction algorithms include SWUNITED (Bradač et al. 2005, 2009), WSLAP+ (Diego et al. 2005, 2007; Sendra et al. 2014; Diego et al. 2016), GRALE (Liesenborgs, De Rijcke & Dejonghe 2006; Liesenborgs et al. 2009), LENSPERFECT (Coe et al. 2008, 2010), LENSTOOL (for weak-lensing mass reconstruction, see Jauzac et al. 2012; Jullo et al. 2014), SAWLENS (Merten et al. 2009, 2011). We refer the reader to Kneib & Natarajan (2011) for a review on cluster lensing and mass modelling.

The two approaches appear to be complementary to model the different regions of the cluster: in the core, the sparse distribution of the strong-lensing constraints calls for a small number of free parameters, while the geometry of the constraints and the

light distribution can give strong priors on the mass distribution, therefore favouring a parametric modelling approach. In the cluster outskirts, the density of constraints is high, and a more flexible free-form model with many mass 'pixels' would allow a better tracing of the potentially irregular matter distribution, as well as the detection of the presence of (lower density than the cluster) infalling substructures (see for example Meneghetti et al. 2017; Remolina González, Sharon & Mahler 2018, for a comparison between the different types of models). Combining these approaches is the obvious next step, and with this in mind, we have developed a new version of the LENSTOOL software, *hybrid-LENSTOOL*, that combines the parametric modelling approach in the cluster core, with a free-form grid model in the outskirts.

A key challenge when modelling galaxy clusters over an extended spatial scale is the nature of the different lensing constraints depending on the cluster region considered. In cluster cores, where the projected surface mass density is high, gravitational lensing is non-linear. This is the strong-lensing regime. Here background galaxies can be multiply imaged in addition to being extremely distorted. In this case, the positions of the different multiple images of a same background galaxy are used to constrain the projected mass distribution of the lens. In the outskirts of clusters, the surface mass density is lower, images of background galaxies are thus only weakly distorted/sheared. There, gravitational lensing is on average mostly linear, this is the weak-lensing regime. A statistical approach is necessary to infer the projected mass distribution. The combination of the two types of constraints permits self-consistent modelling of the overall cluster mass density distribution.

In this paper, we present the new *hybrid-LENSTOOL* algorithm that we developed in light of the large-scale weak-lensing data that are collected via the Beyond Ultra-deep Frontier Fields And Legacy Observations (BUFFALO) survey (GO-15117, PIs: Steinhardt & Jauzac, Steinhardt et al. 2020). With the availability of these high resolution *Hubble Space Telescope* (HST) observations of the outskirts of clusters, we developed a self-consistent model for the mass distribution of galaxy clusters that includes all scales, using a combination of both strong- and weak-lensing constraints. The goal of the BUFFALO survey is to extend the *HST* coverage of the six Hubble Frontier Field (HFF, Lotz et al. 2018) clusters up to ~3/4 × $R_{vir}$, thus allowing us to extend the HFF high-resolution lens models beyond the cluster cores. Furthermore, we note that the current most up-to-date LENSTOOL strong-lensing model for Abell 370, the first cluster fully observed in BUFFALO, requires a relatively strong external shear component in order to minimize the $\chi^2$ as described in Lagattuta et al. (2019). This strongly motivates the need for improvement of the lens modelling algorithm beyond the parametric version of LENSTOOL. Here, we present this new algorithm, *hybrid-LENSTOOL*, and test the mass reconstruction with a simulated cluster, similar to the BUFFALO clusters.

In Niemiec et al. (in preparation), we apply our newly formulated *hybrid-LENSTOOL* to the extended observations of Abell 370 obtained as part of BUFFALO, and show that the simultaneous modelling of both strong and weak lensing allows us to reduce the external shear component necessary in the parametric strong-lensing only mass model. This demonstrates the power of such an algorithm, and the mandatory need to model clusters consistently incorporating data from all scales comprehensively.

The outline of this paper is as follows: in Section 2, we describe the method used to develop *hybrid-LENSTOOL*, combining the parametric model used in the cluster core and the grid in the outskirts, and how these two modelling frameworks are optimized with both strong- and weak-lensing constraints simultaneously. In
Section 3, we test the algorithm on a simulated cluster, and quantify the improvement in mass modelling compared to the previous versions. We summarize the results and conclude with a discussion in Section 4. Throughout this paper, we use a flat $\Lambda$ cold dark matter cosmology with $\Omega_m = 0.27$ and $h_0 = 0.7$.

## 2 METHOD

### 2.1 Lensing formalism

Gravitational lensing formalism is based on the lens equation:

$$ \beta = \theta - \nabla\psi(\theta), \qquad (1) $$

which relates the angular positions of the image and the source, $\theta$ and $\beta$, respectively, through the gradient of the so-called lensing potential, $\psi$, computed at the image position. The lensing potential is defined as

$$ \psi(\theta) = \frac{2 D_L D_{LS}}{c^2 D_S} \phi(\theta), \qquad (2) $$

where $D_L$, $D_S$, and $D_{LS}$ represent the observer–lens, observer–source, and lens–source angular diameter distances respectively, and $\phi(\theta)$ is the projected Newtonian potential of the lens. The Laplacian of the lensing potential is an indicator of the strength of the lens, and it can be linked using the Poisson equation to the projected surface mass density of the lens, $\Sigma(\theta)$, as

$$ \Delta\psi(\theta) = 2 \frac{\Sigma(\theta)}{\Sigma_{\text{crit}}}, \qquad (3) $$

where the critical surface density of the Universe, $\Sigma_{\text{crit}}$, is defined as

$$ \Sigma_{\text{crit}} = \frac{c^2}{4\pi G} \frac{D_S}{D_L D_{LS}}. \qquad (4) $$

The strength of a lens can thus be determined by comparing its surface mass density with the value of the critical surface density, $\Sigma_{\text{crit}}$, at the corresponding source and lens redshifts. For a cluster at $D_L = 1$ Gpc, and sources at $D_S = 2$ Gpc from the observer, the critical surface density is $\Sigma_{\text{crit}} \sim 0.3 \text{ g cm}^{-2}$. For a typical cluster, the core extends up to 0.5 Mpc with a central mass density $\rho_0 \sim 2 \times 10^{-25} \text{ g cm}^{-3}$ (Bahcall 1977). That gives a surface mass density, $\Sigma_{\text{cluster}} \sim 0.3 \text{ g cm}^{-2} \sim \Sigma_{\text{crit}}$. This quick estimation shows that the centre of clusters, where the density is the highest, can present overcritical strong-lensing regions where multiple images and gravitational arcs can be produced. In this regime, the observed position of multiple images relates to the source position through the lens equation (equation 1), which is degenerated in this case as for a given source position, $\beta$, multiple solutions, $\theta$, can exist.

In the regions of the cluster where the density is lower, $\Sigma_{\text{cluster}} \ll \Sigma_{\text{crit}}$, i.e. in the weak-lensing regime, both the shape distortion and magnification of the source are very small. These very weak distortions in the subcritical regions require a statistical approach to measure the lensing signal. The mapping from the un lensed to the lensed coordinates can be described by the Jacobian matrix, $A$, also called the amplification matrix:

$$ A_{ij} = \frac{\partial \beta_i}{\partial \theta_j}. \qquad (5) $$

This matrix can be rewritten as

$$ A = \begin{bmatrix} 1 - \kappa & 0 \\ 0 & 1 - \kappa \end{bmatrix} + \begin{bmatrix} -\gamma_1 & -\gamma_2 \\ -\gamma_2 & \gamma_1 \end{bmatrix}, \qquad (6) $$

where we introduce the convergence, $\kappa$, that describes the magnification of the images, and the complex shear, $\gamma = \gamma_1 + i\gamma_2$, that describes the stretching of the images. These parameters are derived from the lensing potential, $\psi(\theta)$ (Kneib & Natarajan 2011):

$$ \kappa(\theta) = \frac{1}{2} (\partial_1^2 + \partial_2^2) \psi(\theta) \qquad (7) $$

and

$$ \begin{align} \gamma_1(\theta) &= \frac{1}{2} (\partial_1^2 - \partial_2^2) \psi(\theta) \notag \\ \gamma_2(\theta) &= \partial_1 \partial_2 \psi(\theta), \tag{8} \end{align} $$

where $\partial_i$ represents the partial derivative with respect to $\theta_i$. We note that equation (3) shows that the convergence is related to the projected surface mass density as $\kappa(\theta) = \Sigma(\theta)/\Sigma_{\text{crit}}$.

### 2.2 Parametric modelling of the cluster core

In the central regions of galaxy clusters, i.e. in the strong-lensing regime, the geometry of multiple image systems and the distribution of cluster galaxies provide information on the priors for the matter distribution. It is therefore more appropriate to use so-called parametric models, which are described by physical quantities that allow a direct interpretation of the results. As described earlier, the matter distribution in this regime with this approach is then typically decomposed into cluster-scale and galaxy-scale haloes Natarajan & Kneib (1997). The way these haloes are modelled in LENSTOOL is described in detail in Jullo et al. (2007). Briefly, each halo is parametrized by its position in the sky (x, y), projected ellipticity, $e$, and angle position, $\theta$. A number of parametric profiles are available to describe the distribution of dark matter within each halo, such as PIEMD (Elíasdóttir et al. 2007), NFW (Navarro et al. 1996), or SIS. Each profile has a different set of parameters to describe the density slope.

Since the number of strong-lensing constraints is small compared to the number of galaxy-scale haloes, the radial profile of each single galaxy-scale subhaloes cannot be individually constrained. Therefore, to decrease the number of free parameters, the mass of each subhalo is coupled to the luminosity of the galaxy it hosts using a global parametric mass-to-light relation. In practice, as initially proposed by Natarajan & Kneib (1997), subhaloes are described with PIEMD profiles, in which free parameters are the core radius, $r_{\text{core}}$, the cut-off radius, $r_{\text{cut}}$, and the velocity dispersion $\sigma_0$. These parameters are in turn related to the galaxy luminosity, $L$:

$$ \left\{ \begin{array}{l} \sigma_0 = \sigma_0^* \left(\frac{L}{L^*}\right)^{1/4}, \\ r_{\text{core}} = r_{\text{core}}^* \left(\frac{L}{L^*}\right)^{1/2}, \\ r_{\text{cut}} = r_{\text{cut}}^* \left(\frac{L}{L^*}\right)^{\alpha} \end{array} \right. \qquad (9) $$

where $L^*$ is the typical luminosity of a galaxy at the cluster redshift, and $r_{\text{cut}}^*, r_{\text{core}}^*$, and $\sigma_0^*$ are its PIEMD parameters. These are the free parameters used to describe the mass of galaxy-scale subhaloes. The total mass of one subhalo can then be written as

$$ M = (\pi/G) (\sigma_0^*)^2 r_{\text{cut}}^* (L/L^*)^{1/2+\alpha}. \qquad (10) $$

We fix $\alpha = 1/2$, following e.g. Jullo & Kneib (2009) and Richard et al. (2010). In the rest of this paper we denote by $\Theta$ the vector containing the set of free parameters corresponding to the parametric part of the model.
## 2.3 Grid modelling

In the cluster outskirts a more flexible approach is necessary to account for the potentially irregular shape of the cluster, and to allow for effective substructure detection. This can be best achieved using a non-parametric model where the mass distribution is reconstructed using a grid of ‘mass pixels’. As described in Jullo & Kneib (2009) and Jullo et al. (2014), the projected density (or convergence) field is decomposed into a grid of radial basis functions (RBFs). More precisely, a grid covering the field to be modelled is set up, and an RBF is then fixed at each node, described by a truncated isothermal mass distribution (TIMD, the circularized version of a PIEMD).

The true convergence field, $κ(θ)$, is therefore approximated as

$$κ(θ) = \frac{1}{\sum_{i} v_i^2} \sum_{i} f(||θ_i - θ||, s_i, t_i), \qquad (11)$$

where the RBF on grid node $θ_i$ is defined as

$$f(R, s, t) = \frac{1}{2G} \frac{t}{t-s} \left( \frac{1}{\sqrt{s^2 + R^2}} - \frac{1}{\sqrt{t^2 + R^2}} \right). \qquad (12)$$

For a TIMD profile, the weight of the RBF, $v^2$, is the velocity dispersion at the centre of the gravitational potential,<fnref n="1" /> and the RBF parameters, $s$ and $t$, represent the core and cut radii of the profile, respectively. The core radius, $s$, is typically fixed to the distance between two nodes of the grid, and the cut radius to $s = 3t$.

The shear field can now be approximated by the RBFs as

$$\gamma_1(\theta) = \sum_i v_i^2 \Gamma_1(||\theta_i - \theta||, s_i, t_i) \\ \gamma_2(\theta) = \sum_i v_i^2 \Gamma_2(||\theta_i - \theta||, s_i, t_i), \qquad (13)$$

where analytical expressions for $\Gamma_1$ and $\Gamma_2$ can also be derived (see equation A8 in Elíasdóttir et al. 2007).

We denote by $w$ the vector containing all the RBF weights, $v_i^2$, which corresponds to the set of free parameters for the grid model. We set as a prior that the RBF weights, $v_i^2$, are positive, and following Jullo et al. (2014), that they are also described by a Poisson probability distribution function (PDF):

$$P(v_i^2) = \exp(-v_i^2/q)/q, \qquad (14)$$

where $q$ is a nuisance parameter, described by the PDF:

$$P(q) = q_0^2 q \exp(-q/q_0), \qquad (15)$$

where the parameter $q_0$ is fixed to 10 following Jullo et al. (2014).

## 2.4 Likelihood definitions

In the combined mass modelling, the two types of modelling frameworks described above (parametric + grid) are optimized jointly. To perform the optimization, two different types of constraints are used: the strong-lensing regions are constrained by the positions of the multiple images of a same background source, while the weak-lensing regions are constrained by the shapes of the distorted images of the background sources producing the shear field.

In order to consolidate these two approaches to produce a combined mass model, we take the following approach. The free parameters of the parametric part of the model (cluster-scale and galaxy-scale haloes) are arranged in a vector $\Theta$, and the grid model

composed of N RBFs with weight $v_i^2$, are ordered in a vector $w = [v_1^2, ..., v_N^2]$. This allows us to combine and derive the total likelihood describing our model, which is written as

$$\mathcal{L}(\Theta, w) = \mathcal{L}_{\text{SL}}(\Theta, w) \times \mathcal{L}_{\text{WL}}(\Theta, w). \qquad (16)$$

We describe further in this section how we compute the strong- and weak-lensing likelihoods, $\mathcal{L}_{\text{SL}}$ and $\mathcal{L}_{\text{WL}}$, respectively.

### 2.4.1 Computing the strong-lensing likelihood

We consider a set of $M_{\text{SL}}$ background sources strongly lensed so that each source $i$ has $n_i$ multiple images. Considering that the noise in the image position measurement of different images is uncorrelated, the noise covariance matrix is diagonal, and the likelihood can be written as

$$\mathcal{L}_{\text{SL}} = \prod_{i=1}^{M_{\text{SL}}} \frac{1}{\prod_j \sigma_{ij} \sqrt{2\pi}} \exp^{-\frac{\chi_j^2}{2}}, \qquad (17)$$

where $\sigma_{ij}$ is the error on the position of image $j$ of the source $i$. The contribution of a multiple image system $i$ to the total $\chi^2$ can then be expressed as

$$\chi_i^2 = \sum_{j=1}^{n_i} \frac{\|x_{\text{obs}}^j - x^j(\Theta, w)\|^2}{\sigma_{ij}^2}, \qquad (18)$$

where $x_{\text{obs}}^j$ is the measured position of the multiple image $j$, and $x^j(\Theta, w)$ is the position of the image $j$ predicted by the model, in which the free parameters are $\Theta$ for the parametric part, and $w$ for the grid.

In the case of the combined model, we compute the $\chi^2$ in the source plane (for a discussion on the pros and cons of computing the $\chi^2$ in the image or source plane, see Jullo et al. 2007), which gives for one system:

$$\chi_{S,i}^2 = \sum_{j=1}^{n_i} \frac{\|x_S^j(\Theta, w) - \langle x_S^j(\Theta, w) \rangle \|^2}{\mu_j^{-2} \sigma_{ij}^2}, \qquad (19)$$

where $x_S^j(\Theta, w)$ is the position of the source galaxy corresponding to the image $j$ projected to the source plane by the lens equation, $\langle x_S^j(\Theta; w) \rangle$ is the barycentre of the positions of the source corresponding to all the images in system $i$, and $\mu_j$ is the magnification for image $j$. The source position $x_S^j(\Theta, w)$ can be calculated from the measured position of the image, $x_{\text{obs}}^j$, by linearly adding the deflection angle of the parametric model and the RBFs, as

$$x_S^j(\Theta, w) = x_{\text{obs}}^j - \alpha(x_{\text{obs}}^j, \Theta) - \sum_i v_i^2 A(||x_{\text{obs}}^j - x_i||, s_i, t_i), \qquad (20)$$

where $\alpha(x_{\text{obs}}^j, \Theta)$ is the deflection angle produced at the observed image position by the mass distribution included in the parametric model, and $v_i^2 A(||x_{\text{obs}}^j - x_i||, s_i, t_i)$ is the deflection angle produced at the image location by the RBF located at position $x_i$ (see Elíasdóttir et al. 2007, for an analytical expression of $A(r, s, t)$).

### 2.4.2 Computing the weak-lensing likelihood

We now consider a set of $\lambda$ background galaxies, each with a measured ellipticity, $e^i = [e_1^i, e_2^i]$, ordered in a vector of size $2\lambda$, $e = [e_1, e_2]^\top = [e_1^i, ..., e_1^\lambda, e_2^i, ..., e_2^\lambda]^\top$. We denote $\gamma^i = [\gamma_1^i, \gamma_2^i]$, the shear produced on the image of the galaxy $i$ by the grid model. The full shear vector of size $2\lambda$, $\gamma = [\gamma_1, \gamma_2]^\top =$

<footnote n="1">¹The velocity distribution is usually noted $σ_0^2$, but we choose the notation $v^2$ to avoid confusion with the variance $σ$ that appears later in the text.</footnote>
$[\gamma_1^1, ..., \gamma_1^\lambda, \gamma_2^1, ..., \gamma_2^\lambda]$, computed on the locations of the $\lambda$ background galaxies is then

$$ \gamma = M_{\gamma w} w, \qquad (21) $$

where $w$ are the weights of the $N$ RBFs, and $M_{\gamma w} = [\Delta_1, \Delta_2]^\top$ is a $2\lambda \times N$ matrix. Its elements represent the contributions of each unweighted RBF $j$ to the shear of image $i$:

$$ \begin{aligned} \Delta_1^{(i,j)} &= \frac{D_{LSi}}{D_{OSi}} \Gamma_1^i (||\theta_j - \theta_i||, s_j, t_j), \\ \Delta_2^{(i,j)} &= \frac{D_{LSi}}{D_{OSi}} \Gamma_2^i (||\theta_j - \theta_i||, s_j, t_j). \end{aligned} \qquad (22) $$

In the case of the combined model, part of the lens mass distribution is described by the parametric model, and the contribution to the total shear of this mass component must also be taken into account. We denote this shear component with $\gamma'(\Theta)$, and we remind the reader that $\Theta$ comprises the free parameters of the parametric part of the model.

In the linear weak-lensing approximation:

$$ e = 2\gamma + 2\gamma' + n, \qquad (23) $$

where $n$ represents the intrinsic shapes of the galaxies in the source plane, described by a Gaussian distribution with mean 0 and standard deviation $\sigma_n$. We note that the factor 2 comes from the complex ellipticity definition used in LENSTOOL, where its amplitude is expressed using the axial ratio, $r$, as $|e| = (1 - r^2)(1 + r^2)^{-1}$ (corresponding to the $\chi$ notation in Bartelmann & Schneider 2001).

As the intrinsic galaxy ellipticity distribution is considered to be well described by a Gaussian, the weak-lensing likelihood can be expressed as

$$ \mathcal{L}_{\text{WL}} = \frac{1}{Z_L} \exp^{-\frac{\chi_{\text{WL}}^2}{2}}, \qquad (24) $$

where

$$ \chi_{\text{WL}}^2 = (e - 2M_{\gamma w}w - 2\gamma'(\Theta))^\top N_{ee}^{-1} (e - 2M_{\gamma w}w - 2\gamma'(\Theta)), \qquad (25) $$

where $N_{ee} = < ee^\top >$ is the covariance matrix of the measured ellipticities. Following Jullo et al. (2014), the matrix is considered diagonal, and its diagonal elements are expressed as

$$ N_{ee}^{(i,j)} = \sigma_m^2 + \sigma_{\text{int}}^2, \qquad (26) $$

where $\sigma_m$ is the measurement uncertainty, and $\sigma_{\text{int}}$, the galaxy shape noise, which is defined as the scatter in the galaxy intrinsic shape distribution.

The normalization factor is written as $Z_L = \sqrt{(2\pi)^{2\lambda} \det N_{ee}}$.

## 2.5 Implementing the MCMC sampling

The PDFs of the lens model-free parameters are sampled using the Monte Carlo Markov chain (MCMC) algorithm BAYESYS (Skilling 2004) implemented in LENSTOOL, as described in detail in Jullo et al. (2007). In short, the free parameters are sampled with 10 Markov Chains that explore the parameter space following a variant of the Metropolis-Hastings algorithm (Metropolis et al. 1953; Hastings 1970). The selective annealing variant used in BAYESYS ensures a progressive convergence of the chains from the prior to the posterior distribution without being trapped in any local minima.

As described in Jauzac et al. (2012) and Jullo et al. (2014), hybrid-LENSTOOL also includes the BAYESYS extension MASSINF. This extension is useful when some of the free parameters have

a linear contribution to the mass model. It allows us to find their values at each step of the MCMC through a Gibbs sampling, thus drastically decreasing the convergence time compared to a pure BAYESYS sampling (for details see Jauzac et al. 2012). In our case, these linear parameters are the weights of the RBFs (see equation 13).

Throughout the sampling, BAYESYS + MASSINF does not assume that all the RBFs are necessary to reconstruct the mass distribution of the lens, but will rather use a number $n$ of them at each step. This effective number of RBFs is described with a geometric probability distribution:

$$ P(n) = (1 - c)c^{n-1}, \quad \text{with } c = \frac{\alpha}{\alpha + 1}. \qquad (27) $$

Following Jullo et al. (2014), the parameter $\alpha$ is fixed at 2 per cent of the total number of RBFs.

We implement a block-wise sampling in the BAYESYS algorithm, and sample the parametric model and the grid model parameters alternately. In a regular Metropolis-Hasting sampling algorithm, at each step of the MCMC new values for all free parameters are drawn, and then accepted or rejected based on the values of the likelihood function. Given the large number of free parameters, it can become quite difficult, and onerous computing-wise, to find a new acceptable set of parameters. Instead, we here use component-wise sampling, i.e. we split the ensemble of parameters into two blocks, $\Theta$ and $w$, and update the blocks alternately. A more formal description of component-wise sampling in general can be found for instance in Johnson, Jones & Neath (2009).

# 3 TESTS OF hybrid-LENSTOOL ON A SIMULATED CLUSTER

## 3.1 The simulation

We test hybrid-LENSTOOL on a simulated cluster which replicates the mass distribution of the cluster Abell 2744, as described in Jauzac et al. (2016a). This simulated cluster is composed of

(i) two central large-scale potentials in the cluster core – these potentials are modelled with PIEMDs, which parameters are summarized in Table 1; the total mass of the PIEMDs is given by $M_{\text{tot}} = 2\pi\rho_0 \frac{r_{\text{core}}^2 r_{\text{cut}}^2}{r_{\text{core}} + r_{\text{cut}}}$ (Elíasdóttir et al. 2007), where $\rho_0 = \frac{(1.46\sigma_0)^2}{2\pi Gr_{\text{core}}^2}$ (Jullo et al. 2007). This gives $M_{C1} = 1.03 \times 10^{13} M_{\odot}$, $M_{C2} = 8.80 \times 10^{12} M_{\odot}$;

(ii) six large-scale potentials to model surrounding substructures, located within 1 Mpc of the cluster centre. The substructures are also modelled with PIEMDs, with masses $M_N = 2.41 \times 10^{12} M_{\odot}$, $M_{NW} = 2.85 \times 10^{12} M_{\odot}$, $M_{S1} = 7.82 \times 10^{11} M_{\odot}$, $M_{S2} = 6.83 \times 10^{11} M_{\odot}$, $M_{S3} = 1.99 \times 10^{12} M_{\odot}$ and $M_{S4} = 8.98 \times 10^{11} M_{\odot}$;

(iii) a catalogue of 246 galaxy-scale potentials, with parameters $r_{\text{cut}}^* = 14$ kpc, $r_{\text{core}}^* = 0.15$ kpc, and $\sigma_0^* = 155$ km s$^{-1}$. We note that positions and shapes of these potentials correspond to the true galaxy distribution in Abell 2744 as measured in Mahler et al. (2018), and covers only the core region of the simulated cluster.

We use LENSTOOL to compute the deflection and shear maps corresponding to this mass distribution, and create the strong- and weak-lensing constraints by tracing back the positions and shapes of the background sources from the source to the image planes. We thus obtain

(i) 15 multiple image systems in the strong-lensing region covering the redshift range $1.5 < z < 5$, which roughly corresponds to
<table><caption>**Table 1.** Parameters of the large-scale potentials for the simulated cluster used in this analysis. Potentials are modelled with PIEMDs: input (left), best fit from the *Sequential-Fit* (middle) and the *Joint-Fit* (right). The simulated cluster is located at redshift $z = 0.308$.</caption><thead><tr><th rowspan="2"></th><th colspan="8">Input</th><th colspan="2">Sequential-Fit</th><th colspan="2">Joint-Fit</th></tr><tr><th>C1</th><th>C2</th><th>N</th><th>NW</th><th>S1</th><th>S2</th><th>S3</th><th>S4</th><th>C1</th><th>C2</th><th>C1</th><th>C2</th></tr></thead><tbody><tr><td>x (arcsec)</td><td>-2.1</td><td>-17.7</td><td>29.8</td><td>103.0</td><td>-55.5</td><td>-39.4</td><td>139.0</td><td>191.1</td><td>-3.1 ± 2.4</td><td>-16.6 ± 1.5</td><td>-0.2 ± 0.9</td><td>-17.3 ± 0.7</td></tr><tr><td>y (arcsec)</td><td>1.4</td><td>-15.7</td><td>153.3</td><td>84.3</td><td>91.9</td><td>155.7</td><td>95.0</td><td>110.1</td><td>1.8 ± 2.9</td><td>-15.4 ± 1.0</td><td>-0.1 ± 1.1</td><td>-15.7 ± 0.6</td></tr><tr><td>e</td><td>0.82</td><td>0.51</td><td>0.40</td><td>0.60</td><td>0.20</td><td>0.40</td><td>0.30</td><td>0.60</td><td>0.64 ± 0.11</td><td>0.38 ± 0.08</td><td>0.63 ± 0.06</td><td>0.41 ± 0.05</td></tr><tr><td>θ (deg)</td><td>90</td><td>45</td><td>85</td><td>30</td><td>140</td><td>110</td><td>40</td><td>0</td><td>77 ± 8</td><td>39 ± 9</td><td>94 ± 4</td><td>54 ± 3</td></tr><tr><td>r<sub>core</sub> (arcsec)</td><td>18.8</td><td>10.7</td><td>8.4</td><td>8.3</td><td>7.1</td><td>6.2</td><td>6.5</td><td>6.7</td><td>23.3 ± 1.59</td><td>13.8 ± 2.5</td><td>18.4 ± 1.3</td><td>12 ± 1</td></tr><tr><td>r<sub>cut</sub> (arcsec)</td><td>221</td><td>221</td><td>110</td><td>221</td><td>221</td><td>221</td><td>221</td><td>221</td><td>[221]</td><td>[221]</td><td>[221]</td><td>[221]</td></tr><tr><td>σ<sub>0</sub> (km s<sup>-1</sup>)</td><td>607</td><td>743</td><td>439</td><td>480</td><td>272</td><td>292</td><td>454</td><td>300</td><td>759 ± 111</td><td>754 ± 102</td><td>641 ± 38</td><td>762 ± 38</td></tr></tbody></table>

the redshift range of multiple images with measured spectroscopic redshifts in A2744;

(ii) a catalogue of weakly lensed background sources with a density of 45 sources arcmin⁻<fnref n="2" />, and covering the redshift range 0.5 < z < 1.5. This source density and redshift range correspond to the expected depth of the BUFFALO survey. We draw the two components of the intrinsic ellipticities of the galaxies in a Gaussian PDF of width $σ_{int} = 0.27$, while the size is considered constant. The sources are uniformly distributed in the source plane, and we remove the sources located in the strong-lensing region. To estimate the impact of the background source density, we perform a second mass reconstruction, with a weak-lensing source density $N_s = 100$ sources arcmin⁻<fnref n="2" />, which corresponds to a typical source density that can be obtained with deep HST observations (see for example Jauzac et al. 2015).

The projected mass distribution of the simulated cluster is presented in the left-hand panel of Fig. 1, along with the large-scale potentials as white ellipses in the top panel, with substructures named as in Jauzac et al. (2016a). Positions of the grid potentials are shown as white circles in the bottom panel, where the size of the circles is set to the potential core radii of these clumps.

## 3.2 LENSTOOL mass reconstructions

We perform two mass reconstructions of this simulated cluster: (1) the *Sequential-Fit* – we model the strong- and weak-lensing regions successively by first optimizing the parametric model in the cluster core with strong-lensing constraints; then we fix this part to its best-fitting values, and optimize the grid model with weak-lensing constraints; (2) the *Joint-Fit* – we simultaneously optimize the parametric+grid models including both strong- and weak-lensing constraints, following the method presented in this paper.

We note that for both models presented in this section the parameters describing the galaxy-scale potentials are fixed to the input values. To test the implications of such choice, we perform a fit with these parameters set as free, and the *Joint-Fit* does not show any improvement on the recovered values compared to the *Sequential-Fit*.

### 3.2.1 Sequential-Fit

**Strong-lensing modelling.** We first reconstruct the cluster mass distribution following the method described in Jauzac et al. (2015), i.e by modelling successively strong- and weak-lensing regions. Therefore we start by modelling the core of the cluster using strong-lensing constraints only, as if we had no knowledge on the presence of substructures in the outskirts. We reconstruct the mass

distribution by optimizing the parameters describing the potentials C1 and C2. As is often done in similar analyses wherein the modelling is focused on the central region of the cluster, we keep the cut radius fixed to $r_{cut} = 1000$ kpc.

To be able to compare between the two methods we perform this strong-lensing optimization in the source plane. The best-fit values of the model parameters are presented in the middle part of Table 1.

**Weak-lensing component modelling.** As a second step, we add a set of 353 RBFs to model the outskirts of the cluster. They are located at the nodes of a multiscale grid. The resolution of the grid traces the mass distribution of the simulated cluster, with higher resolution in the densest region (for more details on the multiscale grid, see for instance Jullo & Kneib 2009; Jauzac et al. 2012). The resulting grid potentials have core radii varying between $s = 11$ and 87 arcsec (i.e. 50–400 kpc). We remove the RBFs covering the cluster core, and model this central region with the best-fitting values of the parametric model described in the previous paragraph (see Table 1). The grid of RBFs is created in a LENSTOOL input format using our set of publicly available scripts.<fnref n="2" />

We present the projected mass distribution maps resulting from these two successive models in the middle panel of Fig. 1, and refer to it as the *Sequential-Fit*. The presented mass maps are the mean maps computed as the average of 1000 MCMC realizations.

### 3.2.2 Joint-Fit

We now reconstruct the cluster mass distribution by optimizing jointly the parametric and the grid models, with both strong- and weak-lensing constraints. We use the same priors on the free parameters as for the *Sequential-Fit*, and the same multiscale grid. The best-fitting parameters of the large-scale potentials are given in the right-hand panel of Table 1, and the resulting mass map is presented in the right-hand panel of Fig. 1. We refer to this mass reconstruction as the *Joint-Fit*.

## 3.3 Convergence diagnostics

As we modified the MCMC sampling algorithm for *hybrid-LENSTOOL*, we provide in this section a few diagnostics on the convergence of the chains. We use the 10 chains from the *Joint-Fit* with $N_s = 45$ sources arcmin⁻<fnref n="2" />, where each chain was evolved for 100 steps after the burn-in phase.

We first verify the chains have converged by inspecting them visually. We then measure the serial correlation of the chains, and plot the autocorrelation functions for each of them, i.e the amount

<footnote n="2">²https://github.com/AnnaNiemiec/grid_lenstool</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Projected mass maps. Left column: Input simulation, with the large-scale potentials shown as white ellipses (top panel), and the potentials of the multiscale grid shown as white circles, which sizes are set to the core radii, *s*, of the potential (bottom panel). Middle column: *Sequential-Fit*, with source density $N_s = 45$ sources arcmin$^{-2}$ (top) and $N_s = 100$ sources arcmin$^{-2}$ (bottom). Right column: *Joint-Fit*, with source density $N_s = 45$ sources arcmin$^{-2}$ (top) and $N_s = 100$ sources arcmin$^{-2}$ (bottom). The *Sequential-* and *Joint-Fit* mass maps are means over a 1000 MCMC samples.</figcaption></figure>

of autocorrelation between the terms of the chain as a function of the lag, and verify that they rapidly decrease, and become consistent with 0 starting at lag ~5. Finally, we perform a Geweke diagnostic on each chain to check that they have reached a stationary state.

## 3.4 Comparison of the Sequential-Fit and the Joint-Fit mass reconstructions

To compare the reconstructed mass maps obtained with the different models, we compute the normalized residual maps, defined as $(M_{\text{model}} - M_{\text{input}})/M_{\text{input}}$. We present the resulting maps in Fig. 2 for the *Sequential* (left-hand panel) and *Joint-Fit* (right-hand panel), and for source densities $N_s = 45$ sources arcmin$^{-2}$ (top), and $N_s = 100$ sources arcmin$^{-2}$ (bottom). The two models appear to well reproduce the simulated cluster, both in the core and the outskirts. As expected, mass reconstructions with a higher source density are closer to the input, i.e. they are better at detecting lower mass substructures such as S4 and S2, and trace better the true shape of substructures. For the two source densities, the *Joint-Fit* mass reconstructions have a lower overall bias than the *Sequential-Fit* reconstructions.

Indeed, the best-fitting parameters presented in Table 1 indicate that the combined modelling allows us to decrease the bias in the parametric model. When the core of the cluster is first modelled alone as in the *Sequential-Fit*, the mass in the centre can often

be overestimated to compensate for the mass missing in the outskirts. This is reflected in the best-fitting values of the parameters $r_{\text{core}}$ and/or $\sigma_0$ for both C1 and C2 clumps. This gives $M_{\text{C1}} = (2.00 \pm 1.13) \times 10^{13} \text{ M}_\odot$, and $M_{\text{C2}} = (1.17 \pm 1.16) \times 10^{13} \text{ M}_\odot$ for the *Sequential-Fit*, and $M_{\text{C1}} = (1.13 \pm 0.45) \times 10^{13} \text{ M}_\odot$, and $M_{\text{C2}} = (1.04 \pm 0.45) \times 10^{12} \text{ M}_\odot$ for the *Joint-Fit*.

To quantify the model deviations from the input simulation, we measure the projected density profiles for each of the models. For each of the 1000 model realizations corresponding to the MCMC samples, we compute a radial density profile by azimuthally averaging the reconstructed mass maps, taking the cluster centre close to the centre of the mass clump, C1. We then average the 1000 profiles, and compute error bars by taking the standard deviation over all 1000 measurements. Figs 3 and 4 show the density profiles for source densities $N = 45$, and 100 sources arcmin$^{-2}$, respectively. In both figures, the top panel shows the density profiles, while the bottom panel shows the relative deviations of models as lines and relative errors as shaded areas.

For the two background source number densities, the *Sequential-Fit* deviates from the simulation by ~2–3σ on all scales and consistently predicts a higher mass density. As explained before, this is due to the overestimation of the parameters which define the amplitude and size of the central clumps, as this impacts the mass distribution over all scales. Fitting the distribution in the core and in the outskirts self-consistently allows us to avoid this systematic
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Normalized residual maps for the *Sequential-Fit* (left) and *Joint-Fit* (right), and for source density $N_s = 45$ sources arcmin$^{-2}$ (top) and $N_s = 100$ sources arcmin$^{-2}$ (bottom).</figcaption></figure>

effect. Indeed, Figs 3 and 4 show that the *Joint-Fit* gives a density profile within $1-1.5\sigma$ of the true value over all scales. Finally, the mass density profiles confirm that increasing the source number density allows a better recovery of the shapes of substructures, and produces a smoother mass distribution.

We also compute the integrated mass profiles, and derive the total mass within 1 Mpc, for each model. The mass of the simulated cluster is $M(< 1 Mpc) = 9.02 \times 10^{14} M_{\odot}$. As can be inferred from the density profiles, the *Sequential-Fit* overestimates the total enclosed mass, i.e. $M(< 1 Mpc) = (11.52 \pm 0.68) \times 10^{14} M_{\odot}$ [$(12.43 \pm 1.45) \times 10^{14} M_{\odot}$] for $N = 45$ (100) sources arcmin$^{-2}$. The *Joint-Fit* shows a reduced bias on the total mass estimate, and gives values consistent with the true mass: $M(< 1 Mpc) = (9.43 \pm 0.92) \times 10^{14} M_{\odot}$ [$(9.70 \pm 0.82) \times 10^{14} M_{\odot}$] for $N = 45$ (100) sources arcmin$^{-2}$.

A commonly used indicator to quantify the goodness of fit for a strong-lensing model is the root mean square, noted as rms, i.e. the distance between the observed positions of multiple images and the ones predicted by the best-fitting model. We compute the best-fit rms values for the different models, and find that the combined modelling with $N_s = 45$ sources arcmin$^{-2}$ indeed decreases the rms to 0.57 arcsec, compared to 0.63 arcsec for the strong-lensing only model of the cluster core. The rms value is further slightly decreased to 0.54 arcsec for the *Joint-Fit* with $N_s = 100$, which points that improving the modelling of substructures impacts the quality of fit in the cluster core, in agreement with conclusions from Acebron et al. (2017) analysis. We note that the decrease in rms value can only be considered as a significant estimator of the goodness of fit when comparing between different models of the same cluster, as is the case with our analysis.

The last point we would like to address is the speed of the algorithm. The *Joint-Fit* mass reconstruction performed with *hybrid*-

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top panel: Projected mass density profiles of the input simulation (black line), the *Sequential-Fit* (solid blue line), and the *Joint-Fit* (dashed orange line) for the source density $N_s = 45$ sources arcmin$^{-2}$. Bottom panel: Relative deviation for the *Sequential-Fit* (solid blue line) and *Joint-Fit* (orange dashed line). The shaded regions represent the relative errors for both models. The radial range shown by the grey shaded area corresponds to the strong-lensing region of the simulated cluster.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Same as Fig. 3 for source density $N_s = 100$ sources arcmin$^{-2}$.</figcaption></figure>

LENSTOOL is significantly slower than the one performed with the regular version of LENSTOOL. For instance, on a standard laptop, for the models presented here with the source density $N_s = 45$, the *Joint-Fit* can be obtained in 15 mn, and the *Sequential-Fit* in 5 mn (3 mn for the parametric model and 2 mn for the grid). While this result could be expected and is not particularly problematic in this case, this time difference could become more of an issue when modelling very complex systems. In addition, keeping the parameters of the galaxy-scale potentials free further slows down the modelling, and in a larger proportion in the case of the *Joint-Fit* compared to the *Sequential-Fit*. Indeed at each step of the optimization, it requires to compute the shear produced by each of the cluster galaxies (which can be hundreds) at the position of each weakly lensed galaxies
(which can be thousands). This can be made faster by implementing
parallelization in the computation of the weak-lensing likelihood,
as already done for the strong-lensing one.

4 SUMMARY AND CONCLUSIONS

In this paper, we present *hybrid*-LENSTOOL, a new method im-
plemented in the publicly available lens reconstruction algorithm
LENSTOOL. It combines strong- and weak-lensing constraints to self-
consistently reconstruct the cluster mass distribution at all scales.
*hybrid*-LENSTOOL combines a parametric model in the cluster core
with a free-form grid model in the outskirts. It takes advantage
of the complementary strengths of the two types of modelling to
recover the shape and amplitude of the mass distribution with high
precision.

We tested this new method on a simulated cluster composed of a bi-modal mass distribution in the core and six massive substructures in the outskirts. We found that the *Joint-Fit* modelling recovers well the shape and position of the substructures, and gives a more accurate reconstructed mass density profile for the cluster compared to a *Sequential-Fit*, where the core and the outskirts of the cluster are modelled separately (method that was used in the past). In addition, the *Joint-Fit* performs better at predicting the position of multiple images in the cluster core, reducing the rms from 0.63 to 0.57 arcsec.

After demonstrating the power of this new algorithm with simu-
lated data in this method paper, we will present the first application of
*hybrid*-LENSTOOL to real observations in a forthcoming paper
(Niemiec et al., in preparation). As mentioned before, the ongoing
BUFFALO survey (GO-15117, PIs: Steinhardt & Jauzac, Steinhardt
et al. 2020) is extending the *HST* coverage of the six HFF clusters
(PI: Lotz, Lotz et al. 2018), and will complement the strong-
lensing constraints in the core of these clusters with high resolution
weak-lensing data. Combining these datasets, i.e. high resolution
constraints both in the strong- and weak-lensing regions of the
clusters, with *hybrid*-LENSTOOL will produce high-precision models
of the mass distribution of the HFF clusters up to ~ 3/4 R<sub>vir</sub>.
Higher fidelity mass distributions for cluster lenses are important for
utilizing the full potential of clusters to probe dark matter properties
and cluster physics, study the distant Universe that they magnify,
and be used as cosmological probes, as observed lensing effects are
not impacted by the dynamical complexity of clusters.

ACKNOWLEDGEMENTS

AN would like to thank Jonathan Vacher for the helpful discussions. Support for BUFFALO surveys Program number GO-15117 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. MJ is supported by the United Kingdom Research and Innovation (UKRI) Future Leaders Fellowship ‘Using Cosmic Beasts to uncover the Nature of Dark Matter’ [grant number MR/S017216/1]. This project was also supported by the Science and Technology Facilities Council [grant number ST/L00075X/1]. AN, PN, and KS acknowledge support from the BUFFALO programme. ML, EJ, and AN acknowledges CNRS and CNES for support.

REFERENCES

Acebron A., Jullo E., Limousin M., Tilquin A., Giocoli C., Jauzac M., Mahler G., Richard J., 2017, MNRAS, 470, 1809
Alavi A. et al., 2016, ApJ, 832, 56

Atek H. et al., 2015, ApJ, 800, 18
Atek H., Richard J., Kneib J.-P., Schaerer D., 2018, MNRAS, 479, 5184
Bahcall N. A., 1977, ARA&A, 15, 505
Bartelmann M., Maturi M., 2017, Scholarpedia, 12, 32440
Bartelmann M., Schneider P., 2001, Phys. Rep., 340, 291
Bayliss M. B. et al., 2020, Nat. Astron., 4, 159
Bayliss M. B., Rigby J. R., Sharon K., Wuyts E., Florian M., Gladders M. D., Johnson T., Oguri M., 2014, ApJ, 790, 144
Binney J., Tremaine S., 1987, Galactic Dynamics. Princeton University Press
Birrer S. et al., 2019, MNRAS, 484, 4726
Bouwens R. J., Illingworth G. D., Oesch P. A., Atek H., Lam D., Stefanon M., 2017, ApJ, 843, 41
Bradač M. et al., 2006, ApJ, 652, 937
Bradač M. et al., 2009, ApJ, 706, 1201
Bradač M., Schneider P., Lombardi M., Erben T., 2005, A&A, 437, 39
Bradač M., Allen S. W., Treu T., Ebeling H., Massey R., Morris R. G., von der Linden A., Applegate D., 2008, ApJ, 687, 959
Caminha G. B. et al., 2016, A&A, 587, A80
Chisholm J., Rigby J. R., Bayliss M., Berg D. A., Dahle H., Gladders M., Sharon K., 2019, ApJ, 882, 182
Clowe D., Gonzalez A., Markevitch M., 2004, ApJ, 604, 596
Coe D., Fuselier E., Benítez N., Broadhurst T., Frye B., Ford H., 2008, ApJ, 681, 814
Coe D., Benítez N., Broadhurst T., Moustakas L. A., 2010, ApJ, 723, 1678
D'Aloisio A., Natarajan P., 2011, MNRAS, 414, 1628
D'Aloisio A., Natarajan P., Shapiro P. R., 2014, MNRAS, 445, 3581
Diego J. M. et al., 2016, MNRAS, 456, 356
Diego J. M. et al., 2018, ApJ, 857, 25
Diego J. M., Protopapas P., Sandvik H. B., Tegmark M., 2005, MNRAS, 360, 477
Diego J. M., Tegmark M., Protopapas P., Sandvik H. B., 2007, MNRAS, 375, 958
Diego J. M., Broadhurst T., Molnar S. M., Lam D., Lim J., 2015, MNRAS, 447, 3130
Diego J. M., Hannuksela O. A., Kelly P. L., Pagano G., Broadhurst T., Kim K., Li T. G. F., Smoot G. F., 2019, A&A, 627, A130
Eckert D. et al., 2015, Nature, 528, 105
Elíasdóttir Á. et al., 2007, preprint (arXiv:0710.5636)
Goobar A. et al., 2017, Science, 356, 291
Harvey D., Massey R., Kitching T., Taylor A., Tittley E., 2015, Science, 347, 1462
Harvey D., Kneib J. P., Jauzac M., 2016, MNRAS, 458, 660
Harvey D., Robertson A., Massey R., Kneib J.-P., 2017, MNRAS, 464, 3991
Harvey D., Robertson A., Massey R., McCarthy I. G., 2019, MNRAS, 488, 1572
Hastings W., 1970, Biometrika, 57, 97
Hoekstra H., Bartelmann M., Dahle H., Israel H., Limousin M., Meneghetti M., 2013, Space Sci. Rev., 177, 75
Ishigaki M., Kawamata R., Ouchi M., Oguri M., Shimasaku K., Ono Y., 2018, ApJ, 854, 73
Jauzac M. et al., 2012, MNRAS, 426, 3369
Jauzac M. et al., 2015, MNRAS, 446, 4132
Jauzac M. et al., 2016a, MNRAS, 463, 3876
Jauzac M. et al., 2016b, MNRAS, 457, 2029
Jauzac M. et al., 2018, MNRAS, 481, 2901
Johnson A. A., Jones G. L., Neath R. C., 2013, Statistical Science, 28, 360
Johnson T. L. et al., 2017, ApJ, 843, L21
Jullo E., Kneib J.-P., 2009, MNRAS, 395, 1319
Jullo E., Kneib J.-P., Limousin M., Elíasdóttir Á., Marshall P. J., Verdugo T., 2007, New J. Phys., 9, 447
Jullo E., Natarajan P., Kneib J. P., D'Aloisio A., Limousin M., Richard J., Schimd C., 2010, Science, 329, 924
Jullo E., Pires S., Jauzac M., Kneib J.-P., 2014, MNRAS, 437, 3969
Kawamata R., Ishigaki M., Shimasaku K., Oguri M., Ouchi M., Tanigawa S., 2018, ApJ, 855, 4
Kelly P. L. et al., 2015, Science, 347, 1123
Kelly P. L. et al., 2018, Nat. Astron., 2, 334

Kneib J.-P. et al., 2003, ApJ, 598, 804

Kneib J.-P., Natarajan P., 2011, A&AR, 19, 47

Lagattuta D. J. et al., 2019, MNRAS, 485, 3738

Leauthaud A. et al., 2012, ApJ, 744, 159

Leauthaud A. et al., 2015, MNRAS, 446, 1874

Li R. et al., 2016, MNRAS, 458, 2573

Liesenborgs J., De Rijcke S., Dejonghe H., 2006, MNRAS, 367, 1209

Liesenborgs J., de Rijcke S., Dejonghe H., Bekaert P., 2009, MNRAS, 397, 341

Limousin M., Kneib J. P., Bardeau S., Natarajan P., Czoske O., Smail I., Ebeling H., Smith G. P., 2007, A&A, 461, 881

Lotz J. M., et al., 2017, ApJ, 837, 97

Lynds R., Petrosian V., 1986, BAAS, 18, 1014

Mahler G. et al., 2018, MNRAS, 473, 663

Massey R. et al., 2015, MNRAS, 449, 3393

Massey R. et al., 2018, MNRAS, 477, 669

Meneghetti M. et al., 2017, MNRAS, 472, 3177

Merten J. et al., 2011, MNRAS, 417, 333

Merten J., Cacciato M., Meneghetti M., Mignone C., Bartelmann M., 2009, A&A, 500, 681

Metropolis N., Rosenbluth A. W., Rosenbluth M. N., Teller A. H., Teller E., 1953, J. Chem. Phys., 21, 1087

Natarajan P. et al., 2017, MNRAS, 468, 1962

Natarajan P., Kneib J.-P., 1997, MNRAS, 287, 833

Natarajan P., Kneib J.-P., Smail I., Ellis R. S., 1998, ApJ, 499, 600

Natarajan P., Loeb A., Kneib J.-P., Smail I., 2002, ApJ, 580, L17

Natarajan P., Kneib J.-P., Smail I., 2002, ApJ, 580, L11

Natarajan P., Kneib J.-P., Smail I., Treu T., Ellis R., Moran S., Limousin M., Czoske O., 2009, ApJ, 693, 970

Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563

Niemiec A. et al., 2017, MNRAS, 471, 1153

Oguri M., 2010, PASJ, 62, 1017

Remolina González J. D., Sharon K., Mahler G., 2018, ApJ, 863, 60

Richard J. et al., 2010, MNRAS, 404, 325

Rigby J. R. et al., 2018, AJ, 155, 104

Rivera-Thorsen T. E. et al., 2019, Science, 366, 738

Rodney S. A. et al., 2015, ApJ, 811, 70

Rodney S. A. et al., 2018, Nat. Astron., 2, 324

Sendra I., Diego J. M., Broadhurst T., Lazkoz R., 2014, MNRAS, 437, 2642

Sharon K. et al., 2015, ApJ, 814, 21

Sharon K. et al., 2020, ApJS, 247, 12

Sifón C. et al., 2015, MNRAS, 454, 3938

Skilling J., 2004, http://www.inference.org.uk/bayesys/

Smith G. P. et al., 2019, MNRAS, 485, 5180

Smith G. P., Jauzac M., Veitch J., Farr W. M., Massey R., Richard J., 2018, MNRAS, 475, 3823

Soucail G., Fort B., Mellier Y., Picat J. P., 1987, A&A, 172, L14

Soucail G., Mellier Y., Fort B., Mathez G., Cailloux M., 1988, A&A, 191, L19

Steinhardt C. L. et al., 2020, preprint (arXiv:2001.09999)

Suyu S. H., Chang T.-C., Courbin F., Okumura T., 2018, Space Sci. Rev., 214, 91

Teplitz H. I. et al., 2000, ApJ, 533, L65

Treu T. et al., 2016, ApJ, 817, 60

Walsh D., Carswell R. F., Weymann R. J., 1979, Nature, 279, 381

Wong K. C. et al., 2019, preprint (arXiv:1907.04869)

Zitrin A. et al., 2013, ApJ, 762, L30

Zitrin A., Broadhurst T., Bartelmann M., Rephaeli Y., Oguri M., Benítez N., Hao J., Umetsu K., 2012, MNRAS, 423, 2308

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