title: Cluster Strong Lensing with Hierarchical Inference
authors: Pietro Bergamini Adriano Agnello Gabriel Bartosch Caminha
text:  
  
     
  
     Formalism, functional tests, and public code release
  
  
     
     
     Bergamini, Agnello, Caminha
  
  
     
     
     
     
      
      INAF - Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via Gobetti 93/3, I-40129 Bologna, Italy 
     DARK, Niels Bohr Institute, University of Copenhagen,
      Jagtvej 128, 2200 Copenhagen Ø, Denmark  
  adriano.agnello@nbi.ku.dk 
     Kapteyn Astronomical Institute, University of Groningen, Postbus 800, 9700 AV Groningen, The Netherlands 
     Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany 
  
  
     
  
  
  
  
  
  
   
   
  
  
  
   
  
  
  
   
    
    
     Lensing by galaxy clusters is a versatile probe of cosmology and extra-galactic astrophysics, but the accuracy of some of its predictions is limited by the simplified models adopted to reduce the (otherwise intractable) number of degrees of freedom.
    
     We aim at cluster lensing models where the parameters of all cluster member galaxies are free to vary around some common scaling relations with non-zero scatter and deviate significantly from them if and only if the data require it.
    
     We have devised a Bayesian hierarchical inference framework that enables the determination of all lensing parameters and the scaling relation hyperparameters, including intrinsic scatter, from lensing constraints and (if given) stellar kinematic measurements. We achieve this through BayesLens, a purpose-built wrapper around common parametric lensing codes that can sample the full posterior on parameters and hyperparameters, which we release with this paper.
    
     We have run functional tests of our code against simple mock cluster lensing data-sets with realistic uncertainties. The parameters and hyperparameters are recovered within their 68% credibility ranges and the positions of all the “observed” multiple images are accurately reproduced by the BayeLens best-fit model, without over-fitting.
    
     We have shown that an accurate description of cluster member galaxies is attainable, despite a large number of degrees of freedom, through fast and tractable inference. This extends beyond the state-of-the-art of current cluster lensing models. The precise impact on studies of cosmography, galaxy evolution, and high-redshift galaxy populations can then be quantified on real galaxy clusters. While other sources of systematics exist and may be significant in real clusters, our results show that the contribution of intrinsic scatter in cluster member populations can now be controlled.
  
     
  
     Cluster Strong Lensing with Hierarchical Inference
      
     P. Bergamini <ref> A. Agnello<ref> 
  
     G. B. Caminha <ref>,<ref>  
  
            
      December 30, 2023
  =======================================================================================
  
  
  
  
  
  
  § INTRODUCTION
  
  Galaxy clusters (at z≈0.3-0.9) are ideal “cosmic telescopes” to study background galaxies out to z≈7 and they boost the lensing signal of their own galaxies. When galaxies reside in clusters, their lensing cross-section is enhanced, allowing the study of galaxy populations at z≈0.2 0.4, over a wide mass range (down to M_⋆≈ 10^9.5 M_⊙, e.g., ). If the total mass density distributions of lensing clusters are accurately reconstructed, cosmological parameters can also be inferred <cit.>. Dedicated Hubble Space Telescope (HST) surveys, such as the Cluster Lensing And Supernova survey with Hubble <cit.>, the Hubble/Spitzer Frontier Fields program <cit.>, and the Re-ionization Lensing Cluster Survey <cit.>, have enabled the identification and multi-band characterization of tens, sometimes hundreds, of multiple images of background sources, and hundreds of cluster member galaxies per cluster so as to constrain cluster lens models <cit.>. Subsequent multi-slit spectroscopic follow-up observations, such as the CLASH-VLT survey <cit.> has gathered a wealth of multi-object spectroscopic data on cluster member galaxies and background sources, sampling light-cones around 13 clusters (≈20×20 arcmin^2). 
  Due to their extent and magnifying power, galaxy clusters have also been used to study faint and high-redshift sources <cit.>. Since these studies may have implications on our understanding of star-formation at high-redshift and re-ionization, accurate magnification maps are required.
  
  The advent of the Multi-Unit Spectroscopic Explorer <cit.> has enabled integral-field spectroscopy of galaxy clusters. In particular, in recent years numerous MUSE observations have targeted the cores of several massive clusters  <cit.>.
  
  Every MUSE pointing provides a data-cube with a field of view of 1 arcmin^2 and spatial sampling of 0.2. The spectra cover 4750<λ<9350 Å with a resolution of ∼2.6 Å,  almost constant along the whole wavelength range, and a spectral sampling of 1.25 Å/pix (the spectral resolution is robustly characterized, see ). Integral-field spectroscopy of cluster cores enabled a secure identification of tens of multiple images and their redshifts and complete and pure sets of member galaxies.
  
  Besides the identification of cluster members and multiple images, spectroscopy provides further constraints on mass models through kinematics. The homogeneous spectroscopic coverage of CLASH-VLT means that velocities of individual cluster members can be used to aid the cluster profile reconstruction through dynamical modeling <cit.>. Similarly, integral-field follow-up with MUSE yields the internal stellar kinematics of cluster members and allows to independently constrain their mass profiles.
  
  Lensing clusters have undergone an extensive modeling effort by multiple independent teams. Up to now, the constraints have consisted in the positions of tens of multiple images per cluster <cit.>. Mostly, lens models describe these clusters as a superposition of extended dark matter (DM) halos and more localized over-densities corresponding to cluster member galaxies. So-called “non-parametric”, grid-based models <cit.> add mass over-densities only where they are required by the data. Alternatively, parametric models decompose the cluster potential in different (smooth) components that are expected to follow a physical hypothesis on how mass should be distributed in galaxies and clusters <cit.>. In particular, <cit.> have argued that including cluster member galaxies explicitly, like individual components of cluster mass models, is crucial. This has also emerged from the HFF model comparison project <cit.>, where parametric models routinely outperformed current versions of free-form models, both on semi-analytical and on fully simulated mocks. Towards the era of precision cosmology with galaxy clusters, free-form models may eventually replace parametric models, possibly resolving some of their built-in rigidity <cit.>. However, given the current performance of parametric and free-form models, one may safely adopt the working hypothesis that galaxy clusters can be described as a superposition of relatively simple galaxy/halo components.
  
  The bottle-neck in cluster lens models consists in: ≈0.3 0.6 errors in image-position reconstruction <cit.>, significantly higher than current positional uncertainties from HST data, and the appreciable discrepancies in magnification maps produced by different models on the same clusters <cit.>. One cause of such discrepancies may be the rigidity of the models used to reduce the degrees of freedom associated with cluster members, since galaxies are currently modeled as belonging to zero-scatter scaling relations, typically with hyperparameters that are imposed externally <cit.>.
  The freedom and accuracy in lensing models have gained importance in the era of high-accuracy cosmography <cit.>, and on the reliability of magnification maps for the study of high-redshift galaxy populations <cit.>.
  Up to now, departures of individual galaxies from prescribed scaling relations were determined heuristically, on a galaxy-by-galaxy basis <cit.>. Ideally, all galaxies should be let free to vary around some finite-scatter scaling relations, whose parameters (including scatter) are to be determined directly from the data on each given cluster.
  
  The solution is a hierarchical Bayesian inference formalism, where each galaxy has its own associated parameters, and the parameters of all galaxies are posited to be drawn from common relations with hyperparameters to be determined through lens modeling and (if given) auxiliary kinematic information. In this paper, we illustrate this formalism and its application to lensing models of three simple mock clusters of increasing complexity with HFF/CLASH-VLT data quality, accounting for all observational constraints in a self-consistent manner. This ensures that cluster lens models are as flexible as possible, given the data, and that higher accuracy is reached in predicted image positions. Since only one of the contributions to the inference is from lensing, and in order to ensure a fair comparison with state-of-the-art technology, we build our inference as a modular wrapper of common lens modeling codes. This also ensures that further constraints, such as time delays, flux-ratios, and shapes of background sources, can be easily included in the inference.
  Moreover, we perform different functional tests, comparing a hierarchical approach with state-of-the-art zero-scatter models, to assess which sources of systematics can be resolved by this method, and show that the hierarchical models do not over-fit.
  
  This paper is structured as follows. Our hierarchical inference is detailed in <ref>. <ref> covers some technicalities inherent to our code implementation of the models. In <ref>, we perform functional tests
  on three simple but realistic mock clusters. Results are discussed in <ref>, and we conclude in <ref>.
  
  While a cosmological model is adopted to generate mock clusters and fit them, the main results of this work are general and separate from the choice of cosmology. In fact, if needed, cosmological parameters may be sampled as additional hyperparameters in our inference scheme.
     
  
  
  
  
  § A HIERARCHICAL CLUSTER LENSING MODEL
  
  
  In this section, we generalize the formalism of parametric cluster lensing models to allow for hierarchical inference, in particular scaling relations with non-zero scatter. Given the abundance of symbols, in the following (<ref> and <ref>) all model parameters and hyperparameters are marked by a “hat” symbol.
  
  
  
   §.§ Cluster lensing components
  
  In a parametric cluster lensing model, the total mass distribution is typically divided into a few component classes <cit.>. Here, we indicate each class through its contribution, ϕ^class, to the lensing potential. A first, smooth, cluster-scale component (ϕ^halo) accounts for both the DM content of the cluster and the baryonic intra-cluster gas and light contributions. A second “clumpy” component describes the mass in cluster member halos (ϕ^gal), in DM and baryons. A third component accounts for  the  presence  of  massive  structures  in  the  outer  cluster regions, and possibly additional massive halos along the line-of-sight (ϕ^shear+los). Within this model class, the total cluster potential is then:
  
  
      ϕ_tot = ∑_i=1^N^hϕ_i^halo+ ∑_j=1^N^gal_totϕ_j^gal+ ∑_k=1^N^slϕ^shear+los_k,
  
  
  with N^h, N^gal_tot, N^sl, the number of cluster scale halos, cluster members, and shear plus line-of-sight contributions respectively.
  
  The number of constraints given by the observed multiple image positions is usually not sufficient to fit the parameters of mass profiles used to describe each individual galaxy in a given cluster lens model. For this reason, cluster members are usually parameterized as circular dual pseudo-isothermal mass distributions (circular dPIEs, , ), with negligible core radii, whose parameters are related to galaxy luminosities according to fixed scaling relations. 
  
  The circular dPIE is defined by the 3D mass density <cit.>:
  
  
      ρ(r) =  ρ_0/(1+r^2/r^2_core)(1+r^2/r^2_cut),
  
  
  with:
  
  
      ρ_0 =  σ^2_0/2 π Gr_cut+r_core/r_core^2 r_cut.
  
  
  In these equations, σ_0, r_core, and r_core are the central velocity dispersion, the core radius, and the truncation radius of the dPIE, respectively. For a vanishing r_core, r_cut encircles about half of the total 3D dPIE mass and 60% of the projected mass. In the limit r_cut→∞ and r_core>0, the dPIE coincides with the pseudo isothermal mass distribution defined in <cit.> (PIEMD).
  In lensing models of galaxy clusters, the following scaling relations are commonly adopted for the central velocity dispersions (σ_0^gal), core radii (r_core^gal), and truncation radii (r_cut^gal) of cluster members:
  
  
      σ^gal_0,i= σ^ref_0(  L_i/L_0)^α , r^gal_core,i= r^ref_core(  L_i/L_0)^β_core , r^gal_cut,i= r^ref_cut(  L_i/L_0)^β_cut ,
  
  
  without any intrinsic scatter.
  Under the hypothesis of a power-law scaling M_tot,i/L_i∝ L_i^γ of the total mass of a cluster member, and since the total dPIE mass is M_tot = (πσ_0^2 r_cut)/G, 
  the exponents are therefore related through:
  
  
      β_cut=γ-2α+1 .
  
  
  
  
   §.§ Measured kinematics
  
  
  
  Spectroscopic information, when incorporated inside the lensing models, produces robust and accurate reconstructions of the projected cluster masses. 
  However, strong lensing models still suffer from internal degeneracies between the parameters of their mass components (<ref>). A degeneracy exists, for example, between the two normalizations, σ_0^ref and r_cut^ref in <ref>.
  This degeneracy is easily understandable if we consider that the total mass of a dPIE, with negligible r_core, inside an aperture of radius R is <cit.>:
  
  
      M(R)=πσ_0^2/G(R+r_cut-√(r_cut^2+R^2)),
  
  
  and that multiple images constrain only the total mass within their (projected) distance from the center of a galaxy. Other degeneracies also exist between the cluster-scale DM halo parameters and the clumpy subhalo components since the overall mass may be redistributed differently between the main halo and the subhalos, except where the multiple-image positions are strongly constraining.
  
  Recently, velocity dispersions of cluster members (from fits to spectra with signal-to-noise ratio grater than ten) were used to break, or at least reduce, these internal degeneracies ( and more extensively by ). <cit.> and <cit.> fixed the velocity dispersions of several cluster galaxies to their measured values, thereby breaking their model degeneracies by imposing a strong hypothesis on their mass content. In a more flexible approach, <cit.> used high-quality MUSE spectra of cluster member galaxies to build priors on the parameters of the scaling relations in <ref>. 
  
  In this work, we combine both of the above to fully characterize the subhalo components of clusters. Our purpose-built BayesLens wrapper uses the available measured velocity dispersions of the cluster galaxies, σ_m^gal±δσ_m^gal, to infer the hyperparameters (σ̂^ref, α̂) of the σmag scaling relation. 
  A third hyperparameter, Δ̂σ̂^ref, quantifies the scatter of the measured galaxies around the σmag scaling relation. Gaussian priors, centered on the measured σ_m^gal and with standard deviations equal to five times the measured errors δσ_m^gal, are adopted in the lens model for the dPIE velocity dispersions of galaxies with measured velocity dispersions. For galaxies without a measured velocity dispersion, we assumed Gaussian priors centered on the inferred σmag scaling relation and with a standard deviation equal to Δ̂σ̂^ref. We adopt a uniform prior for the reference cut radius of the r_cutmag scaling relation, and β_cut is obtained from <ref>. Unless otherwise stated, all model hyperparameters are left free to vary to fully explore the model posterior probability as described below. All “hat” symbols are introduced more formally in the following subsection (<ref>).
  
  
  
   §.§ Fitting it all together: the posterior
  
  
  
  In our models, we use the measured velocity dispersions of N_m^gal cluster galaxies, σ_m^gal±δσ_m^gal, together with the positions 𝐱_im of N^im multiple images, from N^fam different sources with positions 𝐱_sou, as observational constraints to the lens model free-parameters. Hereafter, these free-parameters will be marked with a hat symbol. 
  In order to explore the lens models, we sample the total posterior probability function p_tot, expressed as the product of: the odds p_g of drawing each galaxy (independently) from a given scaling relation (with intrinsic scatter); the likelihood p_sr of the scaling relation hyperparameters, given the measured velocity dispersions and luminosities; the likelihood p_im of reproducing the observed image positions, given the cluster model (including the parameters of each cluster member galaxy); a prior p_h on the parameters of the main cluster-scale halo(s). We also include a term p_mg that links the dPIE lensing velocity dispersions of the galaxies with measured kinematics to their stellar velocity dispersions, with large uncertainties. Although not formally needed, this term is used only as a loose regularization to ensure convergence.
  
  In synthesis,
  
      p_tot(σ̂^ref, α̂, Δ̂σ̂^ref,r̂_cut^ref,σ̂_m^gal,σ̂^gal,ϕ̂^h| mag^galσ_m^gal, δσ_m^gal,𝐱_im)∝
      ∝ p_sr(σ̂^ref, α̂, Δ̂σ̂^ref, r̂_cut^ref| mag^galσ_m^gal, δσ_m^gal) ×
      × p_mg(σ̂_m^gal|σ_m^gal, δσ_m^gal)× p_g(σ̂^gal|σ̂^ref, α̂, Δ̂σ̂^ref)×
      × p_h(ϕ̂^h)× p_im( 𝐱_im|σ̂^ref, α̂, Δ̂σ̂^ref,r̂_cut^ref,σ̂_m^gal,σ̂^gal,ϕ̂^h,𝐱_sou).
  
  Each of the five factors on the right-hand-side is discussed below. In the following, the quantities referring to cluster galaxies with measured velocity dispersions are marked with the subscript “m”.
  
  Some of the factors (e.g., the one on hyperparameters) can be interpreted as a posterior on some parameters given some observations, which (by the “Bayes chain rule”) can further be used as a prior for the full hierarchical posterior. To avoid confusion between those posteriors and the final posterior, we alternatively refer to them as “term” or “odds” in the following.
  
  
  
    §.§.§ Odds on the scaling relation hyperparameters, p_sr
  
  
  This factor is responsible for the σmag and r_cutmag scaling relation hyperparameters, given the set of N_m^gal measured cluster galaxies. For the galaxy velocity dispersions, we consider here a scaling relation of the same form of <ref> parameterized by the reference measured velocity σ̂^ref and slope α̂ plus an intrinsic scatter Δ̂σ̂^ref in measured velocity dispersions. Regarding the r_cutmag scaling relation, in the current version of our models we optimize only the reference value r̂_cut^ref while the slope β_cut is determined using <ref> from the inferred α̂ and assuming a fixed mass-to-light scaling for the cluster galaxies. No scatter around this relation is considered.
  
  The term p_sr can be expressed as:
  
  
      p_sr(σ̂^ref, α̂, Δ̂σ̂^ref,r̂_cut^ref| mag^gal, σ^gal_m, δσ^gal_m) ∝
       p_0,sr(σ^gal_m| mag^gal, δσ^gal_m, σ̂^ref, α̂, Δ̂σ̂^ref) ×
      ×ϖ_sr(σ̂^ref, α̂, Δ̂σ̂^ref,r̂_cut^ref),
  
  where the uninformative prior ϖ_sr is defined by:
  
      ln{ϖ_sr(σ̂^ref, α̂, Δ̂σ̂^ref)}=
          
              - ln(Δ̂σ̂^ref),   if σ^ref_min<σ̂^ref<σ^ref_max
         and α_min<α̂<α_max
         and Δσ^ref_min<Δ̂σ̂_m<Δσ^ref_max
         and 
              r^ref_cut,min<r̂_cut^ref<r^ref_cut,max
      
              -∞,                 otherwise,
  
  
  and limits the σmag plus scatter scaling relation hyperparameters and the reference truncation radius, r̂_cut^ref, to lie within suitably chosen boundaries σ^ref_min(max), α_min(max), Δσ^ref_min (max), and r_cut,min(max)^ref.
  The log-likelihood term
  
  
      ln{p_0,sr(σ_m^gal| mag^gal, δσ_m^gal, σ̂^ref, α̂, Δ̂σ̂^ref)} =
      = - 1/2∑_i=1^N_m^gal[ (σ_m,i^gal-σ̂_m,i^sr)^2/(δσ^gal_m,i)^2 + (Δ̂σ̂^ref)^2 + .
      . +ln{2 π[(δσ^gal_m,i)^2 + (Δ̂σ̂^ref)^2 ] }],
  
  
  quantifies the goodness-of-fit of scaling relation hyperparameters to the measured velocity dispersions. In <ref>, we define σ̂_m,i^sr as:
  
  
      σ̂_m,i^sr=σ̂^ref 10^ 0.4 α̂ (mag^ref-mag_i^gal),
  
  where mag^ref corresponds to the reference luminosity L_0 in <ref>.
  
  
  
    §.§.§ Term on measured galaxies, p_mg
  
  This term applies only to the N_m^gal galaxies with measured velocity dispersions. It allows for residual uncertainties (e.g., from mass-anisotropy degeneracy and asphericity) in converting measured velocity dispersions in dPIE σ̂_m^gal, which in turn propagate on the scaling relations that all galaxies are drawn from. This is attained by linking the kinematic sigma to the lensing sigma with some tolerance. If this were not done, one might obtain biased results if the kinematic sigma-mag relation has some intrinsic scatter, but the lensing sigma-mag does not. This term is also needed in order to account for systematic uncertainties in the “measured” dispersions, which in the high signal-to-noise ratio regime can dominate over the statistical uncertainties. Erring on the conservative side, we choose the conversion tolerance to five times the uncertainties in the measured velocity dispersion.
  
  Therefore, p_mg alone consists in Gaussian priors centering the cluster member aperture-average velocity dispersions σ̂_m^gal on their kinematic values σ_m^gal. We choose the standard deviation of Gaussian priors equal to five times the error
  on the kinematic measurements δσ_m^gal. This term is given by:
  
      ln{ p_mg(σ̂_m^gal|σ_m^gal, 5 δσ_m^gal)}=
      = -1/2∑_i=0^N_m^gal{(σ̂_m,i^gal-σ_m,i^gal)^2/(5δσ_m,i^gal)^2+ln[ 2π(5 δσ_m,i^gal)^2 ]}.
  
  
  In other words, this term attributes to galaxies their measured velocity dispersions, with very wide tolerance, unless a deviation from these values produces a significant improvement of the lensing model. We emphasize that this term is, strictly speaking, not necessary, and formally it may prevent the finite-scatter models from reducing to zero-scatter models (if this is needed by the data). However, given the typical δσ_m^gal≳15 km/s uncertainties on measured kinematics, the 5 δσ_m^gal term ensures that this factor only acts as a loose regularization, preventing pathological solutions and aiding the convergence of the models.
  
  
  
    §.§.§ Prior on unmeasured galaxies, p_g
  
  This term is a collection of Gaussian priors on velocity dispersion values for the N^gal=N^gal_tot-N_m^gal galaxies without kinematics measurements. Its form is such that the final posterior prefers lens models in which the unmeasured galaxies lie on the σmag scaling relation inferred by p_sr, unless otherwise required by the lensing data:
  
      ln{p_g(σ̂^gal|σ̂^ref, α̂, Δ̂σ̂^ref)} =
      = - 1/2∑_i=1^N^gal{(σ̂_i^gal-σ̂_i^sr)^2/(Δ̂σ̂^ref)^2 + ln[2 π(Δ̂σ̂^ref)^2 ] }.
  
  The σ̂^sr are computed through <ref>, but now considering galaxies without measured velocity dispersions.
  
  
  
    §.§.§ Prior on halo parameters, p_h
  
  This term consists of flat priors on the smooth cluster-scale halo parameters collectively indicated as ϕ̂^h:
  
  
      ln{ p_h(ϕ̂^h)}=
          
              0,   if ϕ_min^h<ϕ̂^h<ϕ_max^h
      
              -∞,                 otherwise
  
  
  If (as follows) these halos are parameterized as PIEMD profiles, p_h is a prior on the sky coordinates x̂^h and ŷ^h, on the ellipticity ê^h, position angle θ̂^h, core radius r̂_core^h, and central velocity dispersion σ̂_0^h.
  
  
  
    §.§.§ Multiple-image likelihood, p_im
  
  This final term is a likelihood function quantifying the agreement between observed and lens model predicted multiple-image positions. 
  Given N^fam sources with N^im_i multiple images associated to the same i-th source (usually called a “family”), we follow <cit.> and express the likelihood as:
  
  
      p_im( 𝐱_im|σ̂^ref, α̂, Δ̂σ̂^ref,r̂_cut^ref,σ̂_m^gal,σ̂^gal,ϕ̂^h,𝐱_sou) =
      = ∏^N^fam_i=0e^-χ^2_i/2/∏_j Δ x_i,j√(2π),
  
  
  where the χ^2_i associated to the i-th family is:
  
  
      χ^2_i = ∑_j=1^N^im_i𝐱_i,j^obs - 𝐱_i,j^pred^2/Δ x_i,j^2,
  
  
  with 𝐱^obs_i,j the observed positions of the multiple images on the lens plane, Δ x_i,j are the isotropic uncertainties on these positions and 𝐱_i,j^pred are the model predicted positions, given the adopted cosmology and the inferred set of model parameters: σ̂^ref, α̂, Δ̂σ̂^ref,r̂_cut^ref,σ̂_m^gal,σ̂^gal,ϕ̂^h.
  
  
  
  
  
  § TECHNICALITIES
  
  
  
  To sample the complete posterior in <ref>, we use the Affine-Invariant sampling as originally introduced by <cit.>, which is especially suited to our highly-dimensional (and possibly degenerate) parameter space. In particular, to enable full portability and reproducibility, we use the latest  release[https://github.com/dfm/emcee] of  <cit.>.
  
  The first four terms in equation <ref> are directly implemented in our code, while to compute the multiple-images likelihood p_im we exploit the publicly available software LensTool (, , ). 
  The synergy between our code and LensTool requires some technicalities, as described in the following subsections. In any case, our code is fully modular, so that LensTool can be replaced by any other parametric lensing software by changing only a few  lines.
  
  
  
   §.§ Calls to LensTool: input and output files
  
  
  For each parameters combination, corresponding to a given walker position inside the parameters space, BayesLens silently calls LensTool to compute the p_im term of the total posterior. Every LensTool call needs a different input file generated in our code by a specific -function. Another function reads the resulting likelihood, computed by LensTool, from an output file. Since all these output files are saved on the disk using the same name, we create folders with unique (random) names to differentiate each LensTool call. These folders are deleted at the end of every p_im computation.
  
  To sample the posterior p_tot, millions of walker positions are required. Thus, millions of input and output LensTool files are quickly created and deleted on the disk during this process. To avoid the disk wear and the bottleneck represented by the process of writing/reading files, part of our computer RAM is reserved for creating a “RAM-disk” where the input and output files are temporarily saved and deleted.
  
  
  
  
  
   §.§ From measured σ_m^gal to LensTool fiducial σ_LT
  
  
  The dPIE is implemented in LensTool through a fiducial velocity dispersion σ_LT related to the 1D central velocity dispersion σ_0 (in <ref>) by: σ_0=√(3/2)σ_LT. 
  To convert the model predicted aperture-average line-of-sight velocity dispersions σ̂_m^gal and σ̂^gal to their fiducial LensTool values (σ̂_(m),LT^gal), we relate them through
  
  
      σ̂_(m),LT^gal=σ̂_(m)^gal/c_p(R) ,
  
  
  where R is the aperture radius chosen for the cluster member spectral extraction.
  
  Adopting orbital isotropy and a (spherical) galaxy surface brightness profile proportional to the dPIE matter density, the projection coefficient c_p(R) is given by:
  
  
      c_p^2(R) = 6/πr_c+r_t/r_c^2 r_t( √(r_c^2+R^2)-r_c-√(r_t^2+R^2)+r_t)^-1×
      ×∫_0^R∫_R'^∞ R' r_tarctan(r/r_t)-r_carctan(r/r_c)/(1+r^2/r^2_c)(1+r^2/r^2_t)√(r^2-R'^2)/r^2drdR',
  
  
  with r_c=r_core and r_t=r_cut <cit.>.
  
  
  
  § FUNCTIONAL TESTS
  
  
  
  To test the ability of BayesLens in recovering the correct halos and subhalos mass parameters and in predicting the correct multiple-image positions, we develop three mock galaxy clusters, namely Cluster-1, Cluster-2, and Cluster-3, with an increasing degree of complexity. Cluster-1 is a simple and well-controlled cluster mock, while Cluster-2 and 3 have the main scope to analyze the role of low-mass cluster substructures and accuracy in the multiple-image positions. Although none of the models claim to reproduce the complexity of real galaxy clusters, they represent ideal preliminary tests for our code, leaving to future works its application to more complex simulations <cit.> and to real cluster observations. In particular, the mocks allow us to verify the flexibility of our hierarchical approach, their dependence on the lensing or kinematics constraints, and their robustness against over-fitting. Since our simple mocks are developed using LensTool, the results we obtain are directly comparable to model inputs, thus avoiding the possibility of parameter misidentifications in the discussion of the results.
  
  In this section, we describe the main characteristics of the simulated clusters and the settings adopted in the BayesLens optimization. We emphasize that these are functional tests, that is, on the behavior of the method itself in recovering the input parameters of simple but realistic mocks. These are the very first tests that any new method undergoes before it is tested on less-well controlled mocks (e.g., from simulations) to ensure that it can properly reproduce its own input parameters. On fully simulated or real-life clusters, the performance of each separate method may worsen because the true underlying mass distribution is unknown “a priori”, so the only meaningful comparison is among different models on well-controlled mocks.
  
  
  
  
   §.§ Cluster-1: the simplest toy cluster model
  
  
  
  Cluster-1 is based on the best-fit LensTool lens model developed by <cit.> (hereafter ) for the CLASH cluster MACS J1206.2-0847 at redshift z=0.439. With respect to , the cluster-scale component of our simulated cluster is made by a single halo, parametrized through a PIEMD mass profile. The center of the halo has an offset of 0.92 from the brightest cluster galaxy (BCG) reference position. It has a core radius r_core^h=3, while its LensTool fiducial velocity dispersion has a value of σ_LT^h=1000. Its ellipticity e=(a^2-b^2)/(a^2+b^2) is fixed to e^h = 0.7, with a position angle of θ^h = 19.14^∘ counterclockwise from west. 
  
  
  
  
  
  
  
  
  
  The clumpy subhalo component of the mock cluster is composed by 138 cluster galaxies, selected from the  members catalog to have magnitudes in the HST/F160W filter given by m_F160W<22. All these galaxies are parametrized as circular dPIE profiles whose r_core^gal, r_cut^gal, and σ_LT^gal values are determined as follows. For r_core^gal and r_cut^gal we adopt the scaling relations in <ref> 
  with slopes β_core=0.5 and β_cut=0.66. The two normalizations, computed at the BCG luminosity L_0 (m_F160W^ref=17.19), are r_core^ref=0.01 and r_cut^ref=5. We assign line-of-sight stellar velocity dispersions to the cluster member galaxies, averaged within a circular aperture, assuming a 15% Gaussian scatter around the scaling relation in <ref> with α=0.27 and normalization equal to 350 at the BCG luminosity. To determine the LensTool fiducial velocity dispersions, we deproject the aperture-averaged velocity dispersions using <ref> and assuming apertures of radius R=0.8. No shear nor foreground structures are present in our Cluster-1 mock.
  
  
  Given the total mass distribution for the mock cluster, we use LensTool to ray-trace the position of 15 background sources, randomly selected from the  catalog, to their multiple images on the lens plane. The sources are within a redshift range of 1.01≤ z≤6.06 and produce a total of
  85 multiple images. Images are identified using an integer number and a letter in such a way that images with identical numbers in their ID belong to the same family. 26 of the 85 multiple images are excluded from the final sample, either because they are de-magnified or because they are too close to a cluster member to be detectable, in real clusters, due to light contamination. We also exclude one of the 15 families of multiple images (i.e., the 10th family adopting image IDs in Fig. <ref>) since it is constituted by a single lensed image. Therefore, our final mock multiple images catalog consists in 58 multiple images, from 14 background sources, shown in Fig. <ref>. The final simulated cluster, despite a purposely simple mass distribution, contains a number of halos, subhalos, and multiple images comparable to those of some CLASH or HFF clusters. 
  
  In our tests, we consider isotropic statistical errors of 0.2on the multiple-image positions. This is a compromise choice between the sharper accuracies that can be achieved (in principle) with state-of-the-art imaging and the tolerance that the models need in order to efficiently explore the full parameter space, avoiding narrow likelihood peaks around local minima. Choosing a 0.2 tolerance on the image positions also allows us to keep conservative uncertainties on the posteriors in the hyperparameters. In practice, one may apply these models in subsequent iterations where the tolerance is gradually reduced to match the astrometric uncertainties.
  
  We also suppose that the stellar velocity dispersions of 58 luminous cluster galaxies are measured within apertures of radius R=0.8.
  This is comparable to the number of velocity dispersions that <cit.> measured from the MUSE data-cube of MACS J1206.2-0847.
  To associate an error to these simulated measures, we use the following empirical relation derived from the measurements by <cit.>:
  
  
      δσ_m^gal/σ_m^gal = 1.6× 10^-3 m_F160W^3-8.53×10^-2 m_F160W^2
      
          +1.509 m_F160W-8.879
  
  
  
  
   §.§ Cluster-2 and Cluster-3: harder rungs with low-mass substructures and systematics in the astrometry of cluster galaxies and multiple images
  
  
  
  Similarly to Cluster-1, the cluster-scale component of both Cluster-2 and Cluster-3 mocks is described through a single elliptical PIEMD profile. The PIEMD has an off-set of about 1.41 from cluster BCG and it is characterized by an ellipticity e^h=0.72, a position angle of θ^h=19.8^∘, a core radius r_core^h=11, and a LensTool fiducial velocity dispersion σ_LT^h=1080 km s^-1.
  The galaxy-scale component of the new mocks extends over a wider magnitude range from the BCG luminosity (m_F160W^ref=17.19) down to m_F160W=24, and amounts to 258 galaxies extracted from the cluster member catalog by . Cluster galaxies are modeled as circular dPIEs whose r_core^gal, r_cut^gal, and σ_LT^gal are determined as in <ref> but assuming α=0.28, σ^ref=295.5 km s^-1, β_cut=0.64, and r_cut^ref=3.74. As in Cluster-1, a 15% Gaussian scatter around the σ^galm_F160W^gal scaling relation is assumed for the line-of-sight projected, aperture average, stellar velocity dispersion of galaxies. With this choice of dPIE parameter values, the cluster member population in Cluster-2 and Cluster-3 ranges from a total mass of 1.6× 10^12 M_⊙, for the cluster BCG, down to a total mass of 1.1× 10^9 M_⊙, for the less massive galaxy of the cluster.
  In the following, we assume to know the measured stellar velocity dispersion and associated errors of 58 cluster galaxies (see <ref>). These velocities are given as inputs to BayesLens for the lens model optimization.
  
  What differentiates Cluster-2 from Cluster-3 is that in the latter model, we introduce, in addition to cluster galaxies, a large population of low mass subhalos inside the simulated cluster. Thus, Cluster-3 is mainly designed to test the impact of numerous, undetected, substructures on BayesLens performances. To assign a mass to the faint subhalos, we exploit the subhalo mass function fitted by <cit.> and adopted by the MOKA algorithm for simulating the gravitational lensing by galaxy clusters <cit.>. In using this formula, we assume a virial mass of M_vir=1.59 × 10^15 M_⊙ for the host cluster halo. This value corresponds to the total mass of MACS J1206-0847 () contained within a distance from the cluster center of r_200c=2.06 Mpc (362.79 at z=0.439) that is expected to be close to the virial radius of the cluster <cit.>. Note that, a circularly symmetric PIEMD mass distribution with r_core^h=11 and σ_LT^h=1080 km s^-1 has a mass, within a radius of 362.79, exactly equal to 1.59 × 10^15 M_⊙.
  
  Low mass subhalos are spatially distributed according to the function derived by <cit.> fitting the results of cosmological simulations. In particular, we adopt a concentration parameter c_200c=r_200c/r_s=3.5 for the host cluster halo and we assume a=1.944 in 3 of <cit.>. Recent work by <cit.> shows a detailed comparison between subhalos in state-of-the-art N-body and hydrodynamical simulations and observed subhalos from a sizable sample of cluster lens models.
  
  In Cluster-3, we consider all the faint substructures within a distance of 89.65 from the cluster BCG (this is the distance of the most external galaxy considered in the model), and with masses ranging between 1.56 × 10^9 M_⊙, that is about the mass of the least massive cluster member, down to 0.93 × 10^9 M_⊙. Under these assumptions, the low mass subhalo population of the Cluster-3 mock amounts to 910 profiles. 
  
  Each low-mass subhalo in the mock is parameterized as a circular dPIE, and to obtain a value for their σ_LT and r_cut parameters, we adopt the following procedure. Firstly, we compute the total mass for all the simulated cluster member galaxies (M_tot^gal). Then, we fit a M_tot^galm^gal_F160W
  relation between galaxy masses and magnitudes. This relation is adopted to associate a fictitious magnitude to each low-mass subhalo. Finally, we use scaling relations in <ref> to derive velocity dispersions, core radii, and truncation radii for the simulated subhalos from their inferred magnitudes. As for cluster galaxies, also for the low-mass subhalos, we assume 15% Gaussian scatter around the σ^galm_F160W^gal scaling relation.
  
  Using Cluster-2 and Cluster-3 mocks, we map the position of 22 background sources, within the redshift range 1.012≤ z ≤ 5.793 and randomly selected from the  catalog, to their multiple images on the lens plane. Both Cluster-2 and Cluster-3 amount to 69 magnified multiple images, whose positions are used as constraints for the lens model optimizations. As in Cluster-1, we still assume an isotropic statistical error of 0.2 on multiple image positions, but an additional systematic uncertainty on the position of multiple images and cluster member galaxies is considered in the Cluster-2, and Cluster-3 mocks. To do so, we assign random Gaussian displacements to the “true” positions, with a dispersion of 0.01 in both directions. These are comparable to the residual uncertainties from the HST-like quality of imaging data. 
  
  In Fig. <ref>, we plot the total projected mass distribution of Cluster-2 and Cluster-3. Cluster member galaxies with measured velocity dispersions are encircled in green, while we mark the position of the “observed” multiple images using cyan crosses. The spatial distribution of the low-mass subhalo population included in the Cluster-3 mock is shown using magenta data-points on the right panel of the figure.
  
  In Tab. <ref>, we report a summary of the main proprieties of Cluster-1,2,3 mocks, such as: the number of galaxies with and without a measured velocity dispersion (N_m^gal and N^gal), the number of low-mass subhalos (N^DH), the statistical and systematic errors on galaxy and multiple images positions (Δ x^gal)_sys, (Δ x^im)_st, (Δ x^im)_sys, and the number of observed multiple images (N_tot^im). 
  
  
  
  
  
  
  
   §.§ BayesLens on the mock clusters
  
  
  In this subsection, we describe the parameter ranges adopted in BayesLens to probe the mass distribution of the Cluster-1, Cluster-2, and Cluster-3 mocks.
  The cluster-scale component of all three lens models is parametrized using a single PIEMD profile. The x, y coordinates of its center can vary within flat priors, 3 wide, centered on the BCG position, while for the ellipticity ê^h, position angle θ̂^h, and fiducial velocity dispersion σ̂_LT^h, we adopt uniform priors inside the following intervals respectively: 0.0 - 0.9, 5 - 35^∘, and 700 - 1300 . Finally, for the PIEMD core radius r̂_core^h, we assume a uniform prior between 1and 7in the Cluster-1 lens model and between 3and 15for Cluster-2 and Cluster-3.
  The subhalo component of lens models is parametrized using the three scaling relations in <ref> normalized at the BCG luminosity. In particular, we fix r̂_core^ref=0.01 and β̂_core=0.5, while the following flat priors are assumed on the others scaling relation parameters: Cluster-1 (0.20 ≤α̂≤ 0.34, 250≤σ̂^ref≤ 450, 1≤r̂_cut^ref≤ 11); Cluster-2, and Cluster-3 (0.21 ≤α̂≤ 0.35, 150≤σ̂^ref≤ 350, 1≤r̂_cut^ref≤ 9). Note that the slope of the r_cutmag scaling relation, β̂_cut, is not a lens model free-parameter since its value is derived from α̂ through <ref> with γ=0.2. Differently from LensTool, in our code, the σmag cluster member scaling relation refers to measured (line-of-sight) stellar velocity dispersions averaged within apertures of radius R=0.8.
  
  The BayesLens model optimization is performed on the lens plane using the positions of the multiple images inside the simulated catalog and the 58 “measured” stellar velocity dispersions of the cluster galaxies. These velocities are used on the one hand to determine the best σmag scaling relation parameters, and on the other, to derive a Gaussian prior for each measured galaxy. Thanks to these priors, the lens model tends to prefer solutions with cluster members velocities that coincide with the measured values (see <ref>), unless the lensing data require them to deviate. The presence of a non-zero scatter also enables the models to fairly sample the whole parameter space, thereby avoiding the underestimation of systematics.
  
  In running BayesLens on the Cluster-1 mock, we leave all galaxies free to scatter around the fitted σmag scaling relation, while in Cluster-2 and Cluster-3 lens models, only the 58 galaxies with a measured velocity dispersion, and the 22 most luminous unmeasured galaxies are left free to vary. The number of free-to-scatter galaxies is one of the inputs of our code, and in general, a lower number of free galaxies sensibly speed-up the optimization of the lens model.
  
  The total posterior probability distribution (p_tot) of the Cluster-1 lens model is sampled using 298 walkers, while we use 500 walkers for the other two models. 
  The initial walker positions are randomly distributed within the flat priors of the cluster-scale parameters (ϕ^h), but they are initialized inside a Gaussian hyper-sphere, around an estimated maximum of the likelihood (in p_sr,p_mg,p_g), in the subspace of the other free-parameters. 
  In the following, we derive the marginalized posterior distribution of model free-parameters by flattening the final Monte Carlo Markov chain (MCMC) of the walkers after removing a sufficiently large burn-in phase. Moreover, to determine the model-predicted multiple image positions and the r.m.s. displacement, Δ^tot_rms, we use the best-fit BayesLens model. This best-fit model is defined as the set of free-parameter values that maximize the total posterior p_tot.
  
  
  
   
  
  
  § RESULTS AND DISCUSSION
  
  
  
  In this section, we describe the main results obtained from an application of BayesLens on the three simulated clusters, namely Cluster-1, Cluster-2, and Cluster-3.
  
  
  
   §.§ BayesLens results for Cluster-1
  
  
  In Fig. <ref>, we show in black the marginalized BayesLens posterior distributions for the free-parameters of the Cluster-1 lens model. The “true” parameter values assumed in the mock cluster are marked with solid red lines.  
  Panels a and b show the posterior distributions of the scaling relation and cluster-scale halo parameters, respectively. In panel c, we plot the marginalized posterior distributions for the stellar velocity dispersions of the four brightest “measured” cluster galaxies: σ̂_m^BCG, σ̂_m^Gal(1), σ̂_m^Gal(2) and σ̂_m^Gal(3). Similarly, panel d shows the velocity dispersions of the three brightest cluster galaxies without measured kinematics (σ̂^Gal(4), σ̂^Gal(5) and σ̂^Gal(6)) and of the lens galaxy,  σ̂^Gal(7), forming the galaxy scale strong lensing system zoomed in Fig <ref>. In these last two panels, we plot also blue solid lines showing the value of galaxy stellar velocity dispersions, as predicted by the best-fit σ^galm_F160W^gal scaling relation. The “true” and BayesLens optimized values of the lens model free-parameters are reported in Tab. <ref>. The latter correspond to the medians of parameter marginalized posterior distributions, while the 16th and 84th percentiles are quoted as an error.
  Thanks to the “measured” cluster members stellar velocity dispersions and to the observed multiple-image positions, BayesLens recovers, well within the 1σ uncertainty, all of the mock true hyperparameters of the scaling relations and the cluster-scale halo parameters.
  To quantify the accuracy of BayesLens in determining the correct stellar velocity dispersions of the cluster members, we plot in Fig. <ref> the aperture average velocity dispersions within aperture of 0.8radius,  σ^gal, as a function of the galaxy magnitudes, m_F160W. The solid black lines correspond to the scaling relation, with the dashed lines marking the 15% scatter in the mock. The (red) errorbars in the upper panel denote the “measured” velocity dispersions with their uncertainties, and the (green) rectangles mark the 16th-to-84th percentiles from the posterior distributions. All of the galaxies lie on their measured velocity dispersions, and the zero-scatter solution still exists within the sampled parameter space. This shows that the loose prior on “measured” galaxies still allows the models to generalize the zero-scatter ansatz. The bottom panel shows the mock and posterior dispersions for the “unmeasured” galaxies. Most of them are compatible with their true values (red dots) within the uncertainties. Most notably, one galaxy-scale strong lensing system (marked as Gal(7), see Fig. <ref> closeup) departs automatically form the zero-scatter best-fit scaling relation and has a dispersion that is very well constrained by the lensing data. The compatibility of mock and posterior dispersions for all galaxies (within the intrinsic scatter) is displayed in Fig. <ref>, showing that there is no appreciable systematic mismatch trend between the “true” and inferred dispersions.
  
  Fig. <ref> shows the cumulative mass profile of the whole cluster, projected along the line of sight, for the mock and the best-fit (maximum a posteriori) model. Even though the total mass of the mock cluster is redistributed in a slightly different way between the cluster halo and the member galaxies, its value at large radii is accurately recovered by the lens model. At the distances where most of the multiple images form, yielding more constraints to the lens model, the model predicts a correct value for the mass profiles within 1%. This is simply because the projected mass M(<R_E)=πΣ_crR_E^2 within the Einstein radius of a lens does not depend on its density profile.
  
  To relate the goodness of fit of the models to a commonly used benchmark, we examine the r.m.s. offset of “observed” vs inferred positions of the multiple images (Δ_rms^tot). Our best-fit hierarchical model for Cluster-1 can reproduce all images with Δ_rms^tot=0.08, while the lens model total chi-square is χ^2_tot=8.66 (see Tab. <ref>). This non-zero Δ_rms^tot shows that our hierarchical model does not over-fit while also reproducing all observables well within their uncertainties. 
  
  The robustness against over-fitting can be understood in two ways. First, if we re-run the models on a mock with 0.05 statistical uncertainties on the image positions, the maximum-a-posteriori parameters do not change, and the r.m.s. offset is still ≈0.08, while only the confidence intervals in the inferred parameters are shrunk. This is because all image-positions in this first mock are kept on their true locations, and the statistical uncertainties are simply a tolerance parameter that enables a smooth exploration of parameter space. Second, from the posterior we see that only about ten galaxies deviate by more than 1σ from the backbone of the scaling relation; if we re-run a LensTool model with only these galaxies free, and all other galaxies on the zero-scatter scaling relation, the χ^2 improves by ≈33, that is three times the change in degrees of freedom.
  
  
  
  
  
   §.§ BayesLens results for Cluster-2 and Cluster-3
  
  
  A comparison between “true” values of the Cluster-2 and Cluster-3 parameters with BayesLens results is in Tab. <ref>. Note that in this table we are showing the same set of free-parameters for the three mocks, however, their “true” values are different in the Cluster-1 and Cluster-2,3 mocks. 
  
  Despite the increasing complexity of the last two mock clusters, BayesLens recovers well within 1σ all cluster-scale halo and scaling relations parameters. Moreover, it finds similar values for the stellar, aperture-average, velocity dispersions of measured and unmeasured cluster galaxies in the two mocks.
  
  In Fig <ref>, we plot marginalized posterior distributions, obtained from the Cluster-3 MCMC chains, for the same set of model parameters reported in Tab. <ref>. 
  The figure shows that the four measured galaxies BCG, Gal(1), Gal(2), Gal(3), and the two unmeasured galaxies Gal(5), Gal(6) have marginalized posterior distributions that pick very closed to the “true” mock stellar velocity dispersions (red vertical lines). Conversely, we observed large deviations from “true” values for Gal(7) and even more for Gal(4). The reason for these deviations is at the basis of the hierarchical approach implemented in BayesLens, and it is due to the interplay between the p_g and p_im terms of the total posterior in <ref>. In fact, the unmeasured galaxy Gal(4) has marginalized posterior distribution centered on the stellar velocity dispersion value predicted by the best-fit σ^galm_F160W scaling relation (marked with a blue line in Fig <ref>). In the absence of constraints from the position of the observed multiple images, the p_g term encourages a zero-scatter lens model solution with all unmeasured galaxy velocity dispersions that lie on the best-fit scaling relation. 
  Similarly, the displacement of the σ̂^Gal(7) posterior distribution from the scaling relation toward the “true” mock value is possible because the increase of the p_im term (due to a lower r.m.s displacement of the images close to Gal(7)) overcome the decrease of p_g that penalizes scattered solutions.
  
  The best-fit Cluster-2 lens model reproduces the position of multiple images with a Δ_rms^tot=0.09, while the total chi-square is χ^2_tot=12.53. As displayed in Tab. <ref>, the addition of the low-mass subhalo population in the Cluster-3 mock produces a small increase, equal to 1.22, of the lens model total chi-square, however, the final Δ_rms^tot of multiple images remains essentially unchanged.
  
  
  
   §.§ Testing against main halo mismatch
  
  
  
  Different mass models, related through unobservable source-position transformations, may be able to reproduce the same image positions, even though these degeneracies are more severe for systems with less symmetry <cit.>. This occurrence is common to all parametric lensing codes. Here, we check whether BayesLens uses the non-zero scatter to over-fit models that are different from the mocks.
  
  To test whether this happens, we also considered one last simulation whose main halo has a Navarro-Frenk-White (hereafter NFW) density profile, while the BayesLens model still adopt a PIEMD to parameterize the main cluster halo. The cluster-scale NFW halo has a total mass M_200c=1.59×10^15 M_⊙ within a radius r_200c=2.06 Mpc. These values correspond to those measured by <cit.>, using weak lensing techniques, for the cluster MACS J1206.2-0847. Assuming a concentration c_200c=3.5, we derive a scale radius r_s=588.6 kpc. The center position of the NFW profile has an offset of 1.4 from the cluster BCG, an ellipticity of 0.4, and a position angle of 19.8^∘ counterclockwise from west. The subhalo component of the cluster and the position of the lensed sources are assumed the same as the Cluster-2 simulation. The lens model for this cluster is constrained by 69 “observed” multiple images, and systematic isotropic uncertainties of 0.01 are assumed on galaxies and image positions (see Cluster-2 model for more details).
  
  The r.m.s. displacement of the best-fit BayesLens model is 0.40, that is about 0.32 higher than the ≈0.08 of a PIEMD model on a PIEMD mock. Some of the observed multiple images are not produced by the best-fit model, whereas the PIEMD models on PIEMD mocks (Cluster-1,2,3) were all able to produce the whole set of images. We re-run BayesLens while also penalizing models for each multiple image that cannot be produced by the model[We add -100 to the log-likelihood for each multiple image that is not produced by the lens model.]. In this case, the best-fit r.m.s. displacement is just 0.35, and all the images are well reproduced by the model.
  
  The slightly higher r.m.s. displacement, as well as the new issues in producing all multiple images, indicate that a main-halo mismatch may be diagnosed when modeling real clusters. BayesLens does not seem to over-fit the most evident systematics, despite higher freedom in the cluster member galaxies.
  
  
  
  
  
  
  § CONCLUSIONS
  
  
  
  We have shown that, despite the (nominally) large number of degrees of freedom in cluster lensing, the exploration of flexible models is feasible thanks to tractable, hierarchical Bayesian inference. Unlike conventional models, where galaxies are placed on a razor-thin scaling relation, and some are “freed” ad hoc, we populate a scaling relation with non-zero scatter, with hyperparameters that are inferred directly from the lensing constraints and (when available) stellar kinematic data. Our tests on a realistic, albeit simplified, mock cluster (namely Cluster-1) show that the improvement over conventional models may be significant: the velocity dispersions of all the galaxies of our mock cluster are reliably recovered, as are the scaling relation hyperparameters, and the r.m.s. displacement of multiple images between measurements and model predictions can decrease appreciably (in the simplest mocks, from 0.40 to 0.08) unless there are other major sources of systematics. As shown by the tests in <ref>, the impact of the non-zero scatter is mainly in the deflections from the cluster member galaxies and cannot correct for a mismatch between the `true' density profile of the main cluster halo and the one used in the models. This also has some practical implications when modeling real-life clusters: if the r.m.s. displacements do not change appreciably from zero-scatter to hierarchical models, this would indicate that other major sources of systematics must be analyzed (e.g., the density profile of the main halo, or major mass contributions along the line of sight). We also remark that these models do not make any claim on the internal mass content of the galaxies themselves, but just on the fact that they do not necessarily have simple PIEMD profiles, that there is a real and measured scatter in their properties, and that the information from lensing is complementary to that from the measured stellar kinematics (i.e., the hyperparameters are free to differ from the ones from stellar kinematics).
  
  We have also tested the robustness of the hierarchical models against two systematic effects: dark subhalos and small systematic offsets in the measured image and cluster member positions (Cluster-2 and Cluster-3 mocks). Within realistic regimes, compatible with current (HST) or upcoming (ELT) depth and image quality, the model performance on the mock cluster does not degrade appreciably.
  
  A crucial consequence of our hierarchical approach is that it also resolves a “false dichotomy” between well-fitting cluster models and realistic population properties of galaxies. 
  Zero-scatter models around a kinematic prior have realistic galaxy populations but a higher r.m.s. offset than models without any kinematic prior, whose scaling relations are significantly discrepant from the “true” input relation.
  In zero-scatter models, there is typically a trade-off between the r.m.s. and recovering realistic galaxy populations <cit.>.
  Our hierarchical models have realistic galaxy populations and produce a small r.m.s. offset (on our mock cluster). With respect to zero-scatter models, BayesLens automatically identifies the galaxy-scale systems where more freedom from the scaling relation backbone is needed in order to reproduce the positions of multiple images around those galaxies with high accuracy. This makes the hierarchical models robust against over-fitting since, in our mocks, only about ten galaxies deviate from the scaling relation backbone by more than 1σ while the chi-square improves by Δχ^2≈33.
  
  We should emphasize that this is a functional test on mocks, and additional effects may play a role in real-life systems, such as massive substructure, deviations from simple geometry of the main DM halo(s) and cluster members, and additional contributions along the line of sight. 
  However, our hierarchical models can eliminate part of the systematic source of uncertainties and show that internal systematics can be controlled.
  
  The freedom in the parameters of each cluster member galaxy, within the intrinsic scatter of their parent population, has multiple implications. First, the data themselves dictate which galaxy should be “freed” and deviate significantly from a baseline scaling relation. Second, intrinsic scatter is a hyperparameter that is left free in the modeling inference, and this allows for a direct determination of galaxy-population properties, whence accurate studies of galaxy formation and evolution. In particular, accurate mass functions of cluster member galaxies can be compared to predictions from cosmological simulations. Also, the freedom in cluster member parameters allows (in principle) for the quantification of environmental effects on galaxies at different locations across the cluster – which up to now has been possible only on stacked weak-lensing measurements <cit.>.
  
  An important feature of our inference is that BayesLens is fully modular. Changing specific python functions within our code, it is possible, in theory, to allow for calls to any chosen parametric lensing software (besides LensTool, used here as the benchmark), as well as to implement different mass profiles. Moreover, different prescriptions to relate the stellar velocity dispersions and lensing parameters of cluster member galaxies can be used. The currently used scaling relations may be easily replaced with fundamental-plane relations (with free hyperparameters) so as to study the evolution of the fundamental plane across redshift and the environment. Its modular structure also enables the use of additional constraints from, for example, flux ratios or extended-source reconstruction. Samples of magnification maps may be produced simply through an other module, which may enable accurate studies of high-redshift galaxy populations.
  
  This is not central to this first paper, so for the sake of brevity, in Fig. <ref>, we show simply a comparison of magnification maps from the best-fit models of the three mock clusters (namely Cluster-1,2,3). The magnification maps are well reproduced, and as can be expected, the appreciable differences happen only across the critical lines and in high magnification regions (|μ|≳20), where the posterior also predicts a wider variance in the magnification. We remark that these models are already an extension over the state of the art, and the availability of samples from the posterior allows one to evaluate magnification uncertainties at all points.
  
  While our set of hyperparameters is currently limited to those of the scaling relations (intercepts, slopes, scatter widths), this can be in principle extended to cosmological parameters, enabling (fully self-consistent) cosmographic measurements.
  Our code is released publicly[The GitHub link is proprietary until acceptance.] and documented at <https://github.com/pietrobergamini89/BayesLens>.
  
  
  PB and GBC thank Piero Rosati, Claudio Grillo, Massimo Meneghetti and Amata Mercurio who have made this work possible. We thank A. Sonnenfeld and M. Lombardi for useful comments on a previous version of this paper. PB acknowledges financial support from PRIN-MIUR 2017WSCC32ASI and through the agreement ASI-INAF n. 2018-29-HH.0. AA was supported by a grant from VILLUM FONDEN (project number 16599). This project is partially funded by the Danish council for independent research under the project “Fundamentals of Dark Matter Structures”, DFF–6108-00470. GBC acknowledge funding from the European Research Council through the Grant ID 681627-BUILDUP, the Max Planck Society for support through the Max Planck Research Group for S. H. Suyu and the academic support from the German Centre for Cosmological Lensing.
  
  
  
  
  
   
  
  
  
  
  
  
  
  
  