© The Authors 2024

1. Introduction

Strong gravitational lensing (SL) describes the relativistic effect where a massive object along our line of sight distorts the light of a background source into multiple arcs. These objects are called lenses and can be single galaxies or even groups and clusters of galaxies. In the case of lensed quasars, we observe multiple images that are especially prominent as the quasars are extremely bright compared to their host galaxies. Because the light rays interact gravitationally with both baryonic and dark matter (DM), SL allows us to study the local properties of the deflector (Kochanek 1991; Rusin & Ma 2000; Cohn et al. 2001; Suyu et al. 2009; Sonnenfeld et al. 2011; Wagner 2019), as well as the total mass distribution or the DM fraction (Bolton et al. 2006; Barnabe et al. 2011; Gavazzi et al. 2012; Sonnenfeld et al. 2015).

An important part of these studies is to obtain a mass model of the deflector by choosing a mass parametrization and optimizing its parameters to fit the data. For galaxy lenses, the most common parametrization is that of a power-law (PL) profile. The PL class of lens models arises for any ensemble of objects that form a structure via gravity as the dominating interaction (Wagner 2020). SL studies on early-type galaxies (ETGs), such as that of the Sloan Lens ACS (SLACS) Survey (Bolton et al. 2006), have provided important insights. One example is the average slope of the total mass density, which is found to be well described by a PL profile with a small intrinsic scatter around γ ∼ 2 (Auger et al. 2010; Barnabe et al. 2011; Shajib et al. 2021). The three-dimensional mass density is ρ(r)∝r^−γ, where γ is the PL slope.

Insights into other lens-mass parameters of ETGs were gained by Bolton et al. (2008), who modeled 63 SLACS lenses with singular isothermal ellipsoids. The values for the external shear of each lens range from 0 to 0.27 (median of 0.05), and the axis ratios range from 0.51 to 0.97 with a median of 0.79 (Gomer 2020). Independent studies of the lens environment determined similar values for the external shear (Wong et al. 2011). Similarly, the distribution of axis ratios from SLACS is consistent with values from studies of nearby elliptical galaxies (Ryden 1992).

The SLACS sample probes a low-redshift range up to z ∼ 0.5 with a median redshift of z ∼ 0.2 (Gavazzi et al. 2012). Another lens sample is from the Strong Lensing Legacy Survey (SL2S), which spans a range of z = 0.2 − 0.9 with a median redshift of z = 0.5, making it a good high-redshift complement to the SLACS survey, as both surveys probe galaxies with similar sizes and velocity distributions (Jackson et al. 2010). One key result from analyzing the SL2S lenses is that the total mass density slopes of individual galaxies do not appear to evolve over time (Sonnenfeld et al. 2015) while the average slope of the population of lens galaxies becomes steeper over time (e.g., Bolton et al. 2012; Sahu et al. in prep.).

Early-type galaxies are of very high interest in astrophysics as their formation and evolution are still unclear. They are characterized by their old stellar populations, red colors, small amount of cold gas and dust, and lack of spiral arms (Cappellari 2016). The current picture is that ETGs are the result of a hierarchical merging scenario: big structures form through the merger of smaller structures (Toomre & Toomre 1972; Kauffmann et al. 1993; White et al. 1991; Cole et al. 2000). During a merger, the central regions of the galaxies can be disrupted by the interaction of two supermassive black holes (SMBHs), which remove stars from the central regions, leading to the formation of a core (Faber et al. 1997; Milosavljević et al. 2002). Another requirement for core formation is the lack of cold gas, which, if present, would lead to a nuclear star burst that would fill the depleted region in the center. Studies show that the mass of the SMBH in the center correlates with the mass absent from the nucleus (Graham 2005; Kormendy & Bender 2009; Rusli et al. 2013; Dullo & Graham 2014). The cores formed in this way in the center of massive ellipticals typically have sizes of between tens of parsecs and a few hundred parsecs (e.g., Dullo & Graham 2014). Cores larger than 1 kpc are rare. The largest observed cores are of the order of 3 kpc; for example, that observed in the massive brightest cluster galaxy of Abell 2261 (Postman et al. 2012; Nasim et al. 2021).

The core is evident in the surface brightness profile of ETGs, with a shallow inner profile and a steeper outer profile (King et al. 1966; Lauer et al. 1995). Insights into the structural parameters and formation histories of ETGs therefore provide a strong test for the ΛCDM cosmological standard model. Recently, observations by the James Webb Space Telescope in the JADES (Bunker 2019) and CEERS (Finkelstein et al. 2022) surveys revealed tensions with the ΛCDM model. The modeled stellar masses in these galaxies at large redshifts (z > 10) are very high and should be much lower in these young galaxies according to current models (Labbè et al. 2023). The agreement with cosmological simulations is better (McCaffrey et al. 2023), but nevertheless galaxy evolution remains a matter of high interest and studies of ETGs play a crucial role.

Strong gravitational lensing can also be used to detect DM subhalos or low-mass dark galaxies (i.e., galaxies with a high DM content and a low surface brightness). Structure formation through hierarchical merging is incomplete: the inner parts of DM halos can survive and remain as subhalos within their new host. CDM simulations confirm these predictions and show that these subhalos have cusped inner density profiles and an upper mass limit of around 10⁶ M_⊙ − 10⁹ M_⊙ (Dalal & Kochanek 2001; Diemand et al. 2008). This predicted mass range is similar to those of faint, DM-dominated dwarf satellites (Simon & Geha 2007; Vegetti et al. 2012).

There are several approaches to using SL for substructure detection. One approach involves the study of flux-ratio anomalies, where the observed flux ratios of the multiple source images differ from those predicted by the smooth-mass model (Mao & Schneider 1997; Metcalf & Madau 2001; Xu et al. 2014; Nierenberg et al. 2017; Gilman et al. 2019; Hsueh et al. 2020). The flux ratios are sensitive to perturbations in the lensing potential and are therefore an indication of small-scale structure in the halo of the lensing galaxy, such as dark matter substructure.

Another approach to the detection of DM substructure uses a method called gravitational imaging, where a best-fitting smooth-lens model is compared to the data and substructures are detectable through anomalies in the surface brightness distribution of the lensed arc (Falco et al. 1999; Koopmans 2005; Vegetti et al. 2010; Vegetti & Vogelsberger 2014; Nierenberg et al. 2014; Ritondale et al. 2019; Despali et al. 2020; Bayer et al. 2023a,b).

Here, we present a modeling analysis of the quadruply lensed quasar HE0230−2130. The configuration of this system is peculiar as there are two lensing galaxies that are possibly part of a galaxy group and may be dynamically interacting to some extent. One of the four quasar images is located between them. For a lensing configuration like this, we would expect five quasar images, and so our goal in this work is to find model parametrizations that can explain a missing quasar image. The lensed quasar PS J0630−1201 (Ostrovski et al. 2018; Lemon et al. 2018) has a very similar configuration to HE0230−2130, but a fifth image is observed. Double galaxy lenses like HE0230−2130 and PS J0630−1201 are uncommon and can cause exotic lens configurations with complex critical and caustic curves (Xivry & Marshall 2009). Understanding these systems might reveal better probes of the underlying lens-mass distribution, or, in the case of high image magnification, may provide ways to resolve even more distant objects than before. Other quasars lensed by binary galaxies are, for example, 2M1310−1714 (Lucey et al. 2018), DES J0408−5354 (Lin et al. 2017; Agnello et al. 2017), and B1608+656 (Myers et al. 1995; Suyu et al. 2009).

We use ground-based imaging data from the Magellan telescope and use the multiple modeled image positions as constraints for the lens-mass distribution. There are HST data available but there is no suitable guide star in the field. A drizzled image of the available unguided and short exposures leads to a noisy image, and therefore the Magellan data are currently the optimal data for our analysis. Approaching this system naively by describing each of the two lensing galaxies with a PL mass distribution results in models that do not fit our observed data: as expected, they always predict one additional, distinctively offset and bright quasar image that is not observed.

To find models that are in agreement with the observations, we analyzed 12 different parametrizations with varying degrees of complexity: singular or cored PL profiles and an additional singular or cored isothermal, spherical mass clump close to the location of the previously predicted fifth image. Each model class is analyzed with and without the presence of external shear. By quantitatively comparing these different classes of models, we can probe possible reasons for the missing fifth image and place constraints on the lens-mass distributions of HE0230-2130.

The outline of the paper is as follows: in Sect. 2 we describe the observations of this lensing system, in particular the imaging data used for the modeling. In Sect. 3, we discuss whether microlensing, dust extinction, or quasar variability could be responsible for the missing image. We then describe our modeling methods and assumptions in Sect. 4. The modeling results are presented and discussed in Sect. 5. In Sect. 6, we summarize our findings.

Throughout the paper, we adopt a flat ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_M = 1 − Ω_Λ = 0.3. Parameter estimates are given by the median of the one-dimensional marginalized posterior probability density functions, and the quoted uncertainties show the 16th and 84th percentiles (corresponding to a 68% credible interval).

2. Observations of HE0230−2130

This quadruply imaged quasar at redshift z_s = 2.162 was discovered by Wisotzki et al. (1999) as part of the Hamburg/ESO survey (Wisotzki et al. 1996). It is located at (right ascension (RA), declination (Dec)) = (38.13829, −21.29046). In 2005, images of this system were taken at the Magellan Clay Telescope with MagIC, which we use here in the modeling. Figure 1 shows a color image, which is a combination of three 400 s exposures in the g-, r-, and i-filter of the Magellan data with a pixel size of 0.069″. The seeing was approximately 0.41″, 0.35″, and 0.35″, respectively.

Fig. 1. Color image of HE0230−2130 created with Magellan imaging data in the g-, r-, and i-band. Image courtesy of Scott Burles.

The brightest lensing galaxy (G1) is located between the multiple lensed quasar images, while the other, fainter galaxy (G2) is slightly offset to the north of image D. The maximum image separation is 2.15″. Faure et al. (2004) found a galaxy overdensity about 40″ south-west of the lens and proposed that G1 and G2 are two members of a physically related galaxy group. This was later confirmed spectroscopically by Eigenbrod et al. (2006), with measured redshifts of z_G1 = 0.523 ± 0.001 and z_G2 = 0.526 ± 0.002. The spectrum of G1 is a good match to that of an ETG. Anguita et al. (2008) obtained similar redshifts for G1 and G2, and found the redshift of a galaxy ∼17″ northwest of image A to be z = 0.518 ± 0.002. To understand whether or not G1 and G2 are potentially dynamically associated with each other we determined whether or not a typical peculiar velocity of cluster galaxies can account for the observed redshift difference of 0.003 (Harrison 1974). We find that we need a peculiar velocity of around 1000 km s⁻¹, which is a typical value for cluster galaxies. Therefore, G1 and G2 seem to be dynamically associated and dynamical interaction is possible to some extent.

3. The missing image

In this section, we discuss some physical reasons that could suppress the fifth image (subsequently referred to as image E) using the predictions from a mass model with two singular PL profiles for the two lens galaxies. Given the predicted proximity of image E to lensing galaxy G2, we need to assess whether microlensing or dust extinction can have a significant effect on the flux of the image and make it unobservable.

3.1. Microlensing

Microlensing is caused by compact objects, such as stars, in the lens galaxy, which can further magnify or demagnify the lensed quasar images. Weisenbach et al. (2021) estimated worst-case microlensing scenarios. Given the median convergence of κ ∼ 0.8 and assuming a stellar mass fraction at the location of image E of ∼0.2, which is realistic for our system, the worst-case magnitude change for this saddle image is Δ₋₋ ∼ 1.5, which corresponds to a magnification ratio of roughly 0.25 between the worst-case demagnified image E and the image E unaffected by microlensing; this means that image E is demagnified through microlensing by a factor of 4 at most.

Figure 2 shows the model-predicted magnifications of image D (red, solid curve), image E (blue, solid curve), and image E in the case of a worst-case microlensing demagnification (blue, dashed curve). The median magnification of image E when unaffected by microlensing is ∼4.8. The median of the distribution of microlensed image E (∼1.2) is lower than the median magnification of image D (∼2.3). Both of these values are above the lowest magnification of ∼0.05 for which an image would be observable in our Magellan image. Therefore, microlensing is an unlikely explanation for the missing image.

Fig. 2. Distributions of predicted image magnifications for the observed image D (red, solid line), the model-predicted image E (blue, solid line), and the demagnified image E in a worst-case microlensing scenario (blue, dashed line). Even if image E were strongly demagnified by microlensing, it would still be observable.

3.2. Dust extinction

Dust in the lensing galaxy can also reduce the observed fluxes of lensed quasar images. Anguita et al. (2008) found that image D is somewhat affected by dust extinction by fitting an extinction law to the flux-ratio spectrum of image A/image D. To check for dust extinction with our Magellan data, we also modeled the lens and quasar light in the g-band and calculated the color excess compared to the i-band. We find no significant color excess compared to the other three quasar images. As image E has a similar distance to G2, we suspect that it is affected by dust to a similar extent to image D. The effect of dust seems to be small, and so even the combined effect of microlensing and dust is not likely to be enough to completely extinguish image E.

3.3. Intrinsic quasar variability

Another possibility is that the natural variability of the quasar and the time delay between the lensed images have conspired to temporarily reduce the flux of image E below the detection threshold. The magnitudes of the quasar images changed by only about Δm ∼ 0.5 from 2004 to 2017 (Millon et al. 2020). For the image to become unobservable, the magnitudes would need to increase by Δm ∼ 5, and so this is also not a possible explanation for the missing image.

4. Modeling method

Having shown that the missing image cannot be explained by microlensing, dust, or variability, as described in the previous section, we looked for answers in modeling this system with different assumptions on the mass profiles. We model HE0230−2130 with the software Gravitational Lens Efficient Explorer (GLEE, Suyu & Halkola 2010; Suyu et al. 2012) and its Bayesian sampling and optimization methods (simulated annealing and MCMC methods), which were automated by Ertl et al. (2023) and Schuldt et al. (2023) to minimize user interaction time. As these automated methods find the optimal sampling parameters and start new chains completely independently, we used them here to speed up the sampling of multiple different model classes.

4.1. Light and mass parametrization

To obtain accurate image positions and priors for the mass parameters of our models, we started by modeling the observed surface brightness of the quasar images and the two lensing galaxies. We selected the i-band of the Magellan data because the lensing galaxies and quasar images are the brightest in this band.

The point-spread function (PSF) is constructed from the cutout (3.38″×3.38″) of a single star in the field at (RA, Dec) = (38.1548, −21.3166). There are three stars in the MagIC field of this observation, but we obtained the best results by choosing only this particular star for the PSF. As the star is not necessarily located at the center of the brightest pixel, we interpolated and shifted the cutout within fractions of a pixel to center the star. We subsampled the PSF by a factor of 3 because the pixel size of 0.069″ is large relative to the full width at half maximum (FWHM) of the PSF (the FWHM of the PSF covers only a few pixels in the native data resolution). Thus, we avoided numerical inaccuracies caused by an undersampled PSF. The light of the quasar images is modeled by fitting this PSF with three parameters: the x- and y-centroid and the amplitude, which is the flux of the respective image.

We parametrized the light of each lensing galaxy with one Sérsic profile (Sérsic 1963) and the quasar images with a point source. Both were convolved with the PSF to fit to the observed image. The Sérsic profile is parametrized as

$$\begin{array}{c}\hfill I_{\mathrm{S}}(x,y)=A_{\mathrm{S}}exp[-k\{{\left(\frac{\sqrt{{(x-x_{\mathrm{S}})}^2+{\left(\frac{y-y_{\mathrm{S}}}{q_{\text{S}}}\right)}^2}}{r_{\mathrm{eff}}}\right)}^{\frac{1}{n_{\mathrm{s}}}}-1\}],\end{array}$$(1)

where A_S is the amplitude, x_S and y_S are the centroid coordinates, q_S is the axis ratio, and n_S the Sérsic index. The constant k normalizes r_eff such that r_eff is the half-light radius in the direction of the semi-major axis. The light distribution is oriented with a position angle ϕ_S that is measured east of north (where ϕ_S = 0° corresponds to the major axis being aligned with the northern direction).

In general, we parametrize the lens-mass distribution of G1 and G2 with a PL profile. The dimensionless surface mass density (or convergence) of a general PL profile is given by

$$\begin{array}{cc}& {\kappa }_{\mathrm{PL}}(x,y)=\frac{3-\gamma }{1+q}\frac{{\theta }_{\mathrm{E}}^2}{{({\theta }_{\mathrm{E}}^2+r_{\mathrm{c}}^2)}^{\frac{3-\gamma }{2}}-r_c^{3-\gamma }}\hfill \\ \hfill & [{(x-x_{\mathrm{G}})}^2+\frac{{(y-y_{\mathrm{G}})}^2}{q^2}+r_c^2{]}^{\frac{1-\gamma }{2}},\hfill \end{array}$$(2)

where (x_G, y_G) is the lens-mass centroid, q is the axis ratio of the elliptical mass distribution, θ_E is the Einstein radius, r_c is the core radius, and γ is the radial profile slope. The mass distribution is oriented with a position angle of ϕ. For r_c = 0″, this reduces to a singular PL elliptical mass distribution (Barkana 1998). The case with γ = 2 and q = 1 corresponds to an isothermal sphere profile. We refer to isothermal spheres with r_c = 0 and r_c > 0 as a singular isothermal sphere (SIS) and a non-singular isothermal sphere (NIS) profile, respectively. We can account for the tidal gravitational field of surrounding objects by adding an external shear component (Keeton et al. 1997). The shear strength is calculated as ${\gamma }_{\mathrm{ext}}=\sqrt{{\gamma }_{\mathrm{ext},1}^2+{\gamma }_{\mathrm{ext},2}^2}$, with γ_ext, 1 and γ_ext, 2 as the components of the shear matrix. The lens potential for the external shear is parametrized by ${\psi }_{\mathrm{ext}}=\frac{1}{2}{\gamma }_{\mathrm{ext}}[\cos (2{\varphi }_{\mathrm{ext}})(x^2-{\mathit{y}}^2)+2\sin (2{\varphi }_{\mathrm{ext}})x\mathit{y}]$, where ϕ_ext is the shear position angle. The position angle is measured east of north, which means that, for ϕ_ext = 0°, the images are stretched vertically, and for ϕ_ext = 90° the images are stretched horizontally.

4.2. Model classes

After obtaining a lens and quasar light fit, we use only the four modeled quasar image positions as constraints for our models as there is no extended structure (lensed arc) visible in the data. In addition, we assume that G1 and G2 are located at the same redshift. From this lensing configuration, we would expect five images. The straightforward approach of modeling this system by parametrizing G1 and G2 with a PL profile and adding external shear always produces one additional, unobserved quasar image about 0.5″ west of G2 (see Fig. 3). Our goal is therefore to find models that can produce the observed number of quasar images of HE0230−2130. To this end, we set up 12 different model parametrizations, each with varying degree of complexity and physical motivation. G1 and G2 are each described either by a singular (PL) or cored PL (cPL) profile with or without external shear. We also consider models with a singular or cored isothermal, spherical mass clump in the region of the predicted image E (see Fig. 3). We call these profiles “restricted” SIS/NIS (i.e., rSIS and rNIS). Table 1 shows an overview and description of the multiple model classes. In all models, we fix the mass centroids of G1 and G2 to the modeled Sérsic centroid of the galaxies. While the axis ratio of the PLs of G1 and G2 are free to vary, we set a Gaussian prior on the position angle based on the light. The prior on the position of the rSIS and rNIS is Gaussian and centered on a typical position of image E. We choose a small σ = 0.14″ ∼ 2 pixels because the mass clumps tend to move away from the region of predicted image E even with a relatively narrow flat prior. The reason for this is that we do not penalize models that produce more than four images during the sampling. A detailed overview of the priors on the mass and shear parameters in each model is shown in Table 2.

Fig. 3. Modeled image of HE0230−2130. Plotted are the 1σ and 2σ contours of the predicted positions of the unobserved image E in the PL1 + PL2 (red), the PL1+cPL2+γext (green), and the PL1+PL2+rNIS (blue) model class. These predicted fifth quasar images would be above the detection threshold, but are not observed in the imaging data.

4.3. Sampling and analysis

We sampled the parameter space of each model class with MCMC chains, where each accepted point in the chain corresponds to one specific model realization. Each chain consists of 10 000 000 accepted points. As it is difficult to incorporate into the sampling process the fact that there is not a fifth, observable image, we use the four modeled image positions as a constraint for our models and subsequently select those model realizations that predict the correct number of images. Thus, we do not sample a distribution of candidate models, but a distribution where some models might fulfill the criterion of producing four observable images. We account not only for models that directly produce only four images, but also additional images that are too faint and are hidden in the noise of the data. To robustly determine those candidate models in the MCMC chains, we computed the image positions and magnifications for each model realization using the Delaunay triangulation method¹. For this, the image plane is divided into triangles. The image plane triangles that, when mapped to the source plane, contain the mean weighted source position are iteratively refined.

If a model produces only four quasar images, or if all additional images are below the flux limit of our data, we count it as a candidate model. We estimated the flux limit in the following way. Because of surface-brightness conservation, the minimum magnification needed for an image to be observable can be calculated via $|{\mu }_{\min }|=|{\mu }_i|\frac{K_{\min }}{K_i}$, where μ_i is the model-predicted magnification of image i, K_i is the modeled amplitude (flux) of an observed image i (that can be image A, B, C, or D), and K_min is the minimum amplitude for a light source to be above the flux limit. We conservatively chose image i for which the predicted magnification limit μ_min is the lowest. We then compared this technical magnification limit with the predicted magnification of additional images. For predicted images that are very close to the center of G2 (i.e., the innermost 3 pixels ≈0.2″), we obtain a separate magnification limit, which is slightly higher due to the residuals from modeling the light of G2. The priors on the mass and shear parameters in each model are listed in Table 2.

5. Results and discussion

5.1. Light parameters of lens and source

We obtained astrometric positions of the QSO images from the lens and quasar light model by fitting the PSF to the observed quasar image light. Figure 4 shows the results of this light fitting; here we compare the data with our model and present the normalized residuals in a range from −5σ to 5σ. The residuals at the quasar positions are due to the fact that the PSF is constructed from only one star and that no corrections were performed.

Fig. 4. Observed and modeled surface-brightness distribution of HE0230−2130 and normalized residuals in a range of −5σ to 5σ.

Table 3 presents the modeled lens light parameters from fitting one Sérsic profile to both G1 and G2. Throughout the paper, we report positions in the coordinate system of CASTLES², where we take our modeled position of image A as coordinate origin.

G1 and G2 show regular isophotes as both can be fit well with one Sérsic profile each. Given this and the fact that we need a relative velocity of between G1 and G2 of around 1000 km s⁻¹ (see Sect. 1), strong dynamical interaction between the two lensing galaxies is not likely, although we cannot rule it out.

Table 4 presents our modeled image positions and compares them with those measured by Gaia (Gaia Collaboration 2018) and by CASTLES with HST data. In this table, we also report measured fluxes from our models. As a magnitude zero point is missing for the Magellan data, we calibrated one using the magnitudes of stars in the field measured in the Legacy Survey DR10 and also calculated their fluxes in the AB magnitude system. The fluxes from CASTLES are reported in the Vega system.

We find that our modeled image positions agree within 6 mas in the x-direction and 5 mas in the y-direction with Castles and Gaia. The root-mean-square (rms) offset with Gaia is ∼2 mas in both x- and y-directions. We therefore added 2 mas in quadrature to the statistical uncertainties from the MCMC chain to obtain the total uncertainty in the astrometry of the quasars. Given this good agreement, we did not conduct corrections to the PSF to reduce the residuals at the quasar images seen in Fig. 4. We exclude only image D in this comparison, which shows a greater offset to the position for Castles because it is the faintest image and is partially blended with the light of G2.

5.2. Candidate models

We present the final mass and shear parameters of all model classes in Table A.1. We list the distributions of parameters for both the full chain and the candidate model realizations. Table 5 shows an overview of our results. The modeled position of image A is set as the origin of the coordinate system.

Of the 12 model classes we considered, 6 can produce models that match the observed image multiplicity. We find that the simplest models, where G1 and G2 are parametrized by a singular PL profile with or without external shear, consistently produce exactly one additional image, as expected from this lensing configuration. These additional images are located about 0.5″ west of G2, with a magnification that should make it observable. Figure 3 shows the 1σ and 2σ contours of the predicted positions of image E in the PL1 + PL2, the PL1+cPL2+γ_ext, and the PL1+PL2+rNIS model classes.

In the class PL1 + cPL2 + γ_ext, 0.26% of the models produce the correct, observed image multiplicity, while for PL1 + cPL2 (same parameterization for G1 and G2 but without external shear), there are no models that match the observation. Similarly, for the model class cPL1 + cPL2 + γ_ext, 0.26% of the models fit the observations, while there are no candidate models in cPL1 + cPL2.

A similar effect might be produced by placing an SIS or NIS mass distribution close to the saddle point of image E on the lens plane. This mass can be interpreted as DM substructure or an unobservable dark galaxy. We find that by simply placing a SIS or NIS clump without positional prior, they tend to move away from the region close to the saddle point and predominantly end up west of image D. Therefore, these two model classes are unable to suppress image E and we introduce the rSIS and rNIS mass clumps, which have a Gaussian prior on their centroids to keep them closer to the region west of G2.

G1 and G2 are parametrized as singular PL profiles. Of the PL1 + PL2 + rSIS models, 0.03% do not predict a fifth image. This is also the case for 1.7% of the PL1 + PL2 + rNIS models.

In our analysis, we assumed the lensing galaxies to be close to isothermal based on the analysis of SLACS lenses (see Sect. 1 with a Gaussian prior centered on γ = 2.08). We tested the influence of this prior on our results by using a shallower slope of γ = 1.8. We find that the exact candidate fraction is sensitive to the prior on γ but our conclusions remain unchanged: a singular PL model for G1 and G2 cannot explain the observed image multiplicity, while models with cored PLs or a rSIS and rNIS model fit the observations.

Finally, we analyzed the PL1 + PL2 + rSIS and PL1 + PL2 + rNIS models with the presence of external shear. We find that the addition of external shear significantly increases the fraction of candidates to 0.4% in the PL1 + PL2 + rSIS+γ_ext case, and to 12% for the PL1 + PL2 + rNIS+γ_ext class. In the following, we provide a detailed description and interpretation of the results for each model class that produces the observed four images.

5.2.1. PL1 + cPL2 + γ_ext

In this model class, G1 is described by a singular PL profile, while G2 is parametrized as a cored PL profile. Galaxies that are part of clusters often exhibit a core, while a core is not common for isolated field galaxies. An additional external shear component is necessary to reproduce the observed image multiplicity, as the PL1 + cPL2 model class always predicts an additional image. By adding a core to the mass distribution of G2, the saddle point of image E merges with the maximum of G2. In Fig. 5, we plot the Fermat potential of two noncandidate and one candidate model; that is, before and after merging of the saddle with the maximum. The difference between candidate and noncandidate is also evident in the source position with respect to the caustic curves in the source plane: as the radial caustic shrinks and changes in shape, the source crosses the radial caustic into the shaded green region, and the image multiplicity of this model changes³. Furthermore, the Fermat potential is flattened in the vicinity of the vanished image.

Fig. 5. Change of critical curves and caustics between noncandidate models (left and middle columns) and one candidate model (right column). The top and middle rows show the critical curves (blue, solid lines) and caustic curves (red, dashed lines). The bottom row shows the Fermat potential and contours with the color scale in units of arcsec2. As the core radius of G2 increases and the source position crosses the radial caustic into the shaded green region, the time-delay saddle point associated with image E merges with the maximum of G2 and image E disappears. The source position is plotted as a yellow star. The predicted image positions are shown as green crosses in each panel.

Figure B.1 shows the posterior distribution of selected lens-mass and shear parameters. The distributions of candidate models are plotted in red contours, and the distribution of the whole chain (i.e., both candidate and noncandidate models) is plotted in black contours. We find that candidate models have a very round mass distribution of G2 with $q_{\mathrm{G}2,\mathrm{cand}}=0.{97}_{-0.05}^{+0.02}$, while the general distribution of the chain is more uniformly distributed over the whole prior range: $q_{\mathrm{G}2,\mathrm{all}}=0.{77}_{-0.12}^{+0.15}$. Similarly, the general distribution of r_c, G2 is skewed towards a vanishing core radius (i.e., a singular PL profile). In order to produce only four images, a core radius of $r_{\mathrm{c},\mathrm{G}2,\mathrm{cand}}={(0.{14}_{-0.04}^{+0.04})}^{″}$ is needed. This corresponds to a physical core size of $r_{\mathrm{c},\mathrm{G}2,\mathrm{cand}}=(0.{88}_{-0.25}^{+0.25})$ kpc, which is in good agreement with typical observed core sizes in elliptical galaxies (see Sect. 1). The candidate models are all found in the tail of the distribution. The shear strength of candidate models is ${\gamma }_{\mathrm{ext},\mathrm{cand}}=0.{08}_{-0.02}^{+0.02}$. This is a typical value for quadruply lensed quasars as shown by Luhtaru et al. 2021.

5.2.2. cPL1 + cPL2 + γ_ext

Now we parametrize both G1 and G2 with a cored PL profile. As for the PL1 + cPL2 + γ_ext case, we need external shear to reproduce the correct image multiplicity. The mechanism by which image E is eliminated is the same as in the PL1 + cPL2 + γ_ext model class: the saddle point of image E merges with the maximum of G2. We show the posterior distributions of all models compared to the candidates in Fig. B.2. Again, a clear tendency of the candidate models is a high q_G2, cand compared to the more uniform, general distribution. The distribution of candidate core radii are $r_{\mathrm{c},\mathrm{G}1,\mathrm{cand}}=(0.{16}_{-0.06}^{+0.04})$ ″ $\sim (1.{00}_{-0.38}^{+0.25})$ kpc and $r_{\mathrm{c},\mathrm{G}2,\mathrm{cand}}=(0.{13}_{-0.04}^{+0.04})$ ″ ∼ $(0.{82}_{-0.25}^{+0.25})$ kpc, respectively. The external shear has a distribution of ${\gamma }_{\mathrm{ext},\mathrm{cand}}=0.{08}_{-0.02}^{+0.01}$.

5.2.3. PL1 + PL2 + rSIS and PL1 + PL2 + rNIS

The physical motivation behind the model classes with an SIS or NIS profile is to simulate a dark mass distribution in the plane of the lensing galaxies; for instance, DM substructure or a dwarf galaxy. The PL1 + PL2 + SIS and PL1 + PL2 + NIS model classes cannot reproduce the correct number of images as the centroids of the mass clumps are moving away from the critical region where the saddle of image E can be merged with the maximum of G2. We assume a Gaussian prior on the centroid of the SIS/NIS clump to keep it close to the critical region shown in Fig. 3. As mentioned earlier, we refer to these as restricted SIS and NIS profiles, or rSIS and rNIS. While the PL1 + PL2 + SIS and PL1 + PL2 + NIS model classes do not predict models with only four images, 0.03% of PL1 + PL2 + rSIS models produce no fifth image. For the PL1 + PL2 + rNIS model class, this percentage is 1.7%. We can see in Figs. B.3 and B.4, that, with the presence of an additional mass component close to G2, the Einstein radius θ_E, G2 decreases. The PL1 + PL2 + rSIS candidate models predict an especially low Einstein radius of G2 with ${\theta }_{\mathrm{E},\mathrm{G}2,\mathrm{cand}}={(0.{19}_{-0.03}^{+0.02})}^{″}$, while it is higher for the candidates in the rNIS model: ${\theta }_{\mathrm{E},\mathrm{G}2,\mathrm{cand}}={(0.{33}_{-0.14}^{+0.11})}^{″}$. The Einstein radius of the rSIS in candidate models is ${\theta }_{\mathrm{E},\mathrm{rSIS},\mathrm{cand}}={(0.{18}_{-0.05}^{+0.05})}^{″}$, which is of a similar size to θ_{E, G2, cand}. The PL1 + PL2 + rNIS candidates have an rNIS Einstein radius of ${\theta }_{\mathrm{E},\mathrm{rNIS},\mathrm{cand}}={(0.{13}_{-0.09}^{+0.12})}^{″}$, which is smaller than θ_{E, G2, cand} in this model class. It is also smaller than the inferred core radius of $r_{\mathrm{c},\mathrm{rNIS},\mathrm{cand}}={(0.{34}_{-0.13}^{+0.19})}^{″}$.

While subhalos are modeled with a Navarro-Frenk-White (NFW; Navarro et al. 1996, 1997) profile in many studies, we chose to model the subhalo with an isothermal sphere profile to keep the number of free parameters as small as possible. Other studies use a Pseude-Jaffe (PJ; Dalal & Kochanek 2001; Vegetti et al. 2014; Despali et al. 2018) profile, which is an isothermal profile with a truncation radius. We conducted a quick model run and find that both the NFW and PJ models are also able to suppress image E with a similar candidate fraction in the MCMC chains. The success of a model with an additional mass clump is therefore independent of the specific choice of mass profile. Furthermore, we do not have any observational constraints on the position of the subhalo along the line of sight, which was shown to be degenerate with the mass of the subhalo by Vegetti & Vogelsberger (2014) and Wagner (2018). The lensing effect more robustly depends on the projected mass of the subhalo. Minor et al. (2017) showed that the projected mass can be robustly inferred even if the mass density profile of the subhalo is not known. We discuss the projected masses of the rNIS clumps in Sect. 5.5.

5.2.4. PL1 + PL2 + rSIS + γ_ext and PL1 + PL2 + rNIS + γ_ext

By adding external shear to the PL1 + PL2 + rSIS and PL1 + PL2 + rNIS model classes, the fraction of candidate models increases substantially. The effect of adding external shear to our models can be seen in the corner plots in Figs. B.4 and B.6. The physical meaning of external shear has recently been challenged, as shear values inferred from strong lensing are inconsistent with independent measurements from, for example, weak lensing (Etherington et al. 2023). The fact that the inclusion of external shear in our models generally increases the candidate fraction may therefore mean that there is some complexity in this lensing system that cannot be accounted for with just two (cored) PL profiles for G1 and G2. This is not surprising as G1 and G2 are possibly embedded in a common DM halo. In particular, 12% of the models of the PL1 + PL2 + rNIS + γ_ext class produce the correct number of images. The 1σ and 2σ contours of the PL1 + PL2 + rNIS centroids are plotted on the modeled image in Fig. 6. This shows that the mass clumps in models with the correct image multiplicity are located about 0.5″ northwest of G2. The posterior distributions in Fig. B.6 show that the PL1 + PL2 + rNIS + γ_ext candidates have an rNIS Einstein radius of ${\theta }_{\mathrm{E},\mathrm{rNIS}+{\gamma }_{\mathrm{ext}},\mathrm{cand}}={(0.{24}_{-0.16}^{+0.17})}^{″}$ and a core radius of $r_{\mathrm{c},\mathrm{rNIS}+{\gamma }_{\mathrm{ext}},\mathrm{cand}}={(0.{24}_{-0.12}^{+0.16})}^{″}$. The external shear strength is ${\gamma }_{\mathrm{ext},\mathrm{cand}}=0.{05}_{-0.01}^{+0.02}$. As in the PL1 + PL2 + rSIS model, the Einstein radius θ_{E, rNIS + γext, cand} is of similar size to θ_{E, G2, cand}, which implies that the mass clump should be sufficiently massive to be luminous.

Fig. 6. 1σ and 2σ contours of the mass centroid of the candidate rNIS clumps in the PL1 + PL2 + rNIS model class. These candidate models predict only four quasar images, or additional images that are too faint to be observed (i.e., below the flux limit of the data).

5.3. Model comparison

As there are multiple model classes that are able to produce the correct number of observed images, we investigated which of the model classes is more likely by comparing their likelihoods. Rather than using the minimum ${\chi }_{\text{im}}^2$ for the comparison, we calculated the Bayesian information criterion (BIC) to weigh the model according to both their complexity and goodness of fit. The BIC is defined as

$$\begin{array}{c}\hfill \text{BIC}=kln(N_{\mathrm{data}})-2ln(\widehat{\mathcal{L}}),\end{array}$$(3)

where k is the number of free parameters, N_data is the number of data points, and $\widehat{\mathcal{L}}$ is the maximum likelihood of the candidate models. The relative probability of a model M with BIC compared to the most likely model M^* with BIC^* is

$$\begin{array}{c}\hfill p(M)/p(M^{\ast })=\text{exp}(-\frac{\text{BIC}-{\text{BIC}}^{\ast }}{2})\times \frac{f_{\mathrm{cand},M}}{f_{\mathrm{cand},M^{\ast }}},\end{array}$$(4)

where f_cand is the weighted fraction of candidate models in the MCMC chain. In terms of free parameters k, we have the free lens galaxy model parameters, plus two for the source position. Because we calculate the flux limit that is specific for each model class and which we use to select candidate models, we count the source intensity as an additional free parameter and the image amplitude K_i (see Sect.4.3) as an additional data point to the eight data points from the four image positions. Additionally, we count the Gaussian priors as both free parameter and data point. Table 6 lists the BIC values and relative probabilities of our six candidate models.

The most likely model class is PL1 + PL2 + rNIS, which is slightly favored over the model class PL1 + cPL2 + γ_ext, which has a relative probability of 94%. Although cPL1 + cPL2 + γ_ext has a lower ${\chi }_{\text{im}}^2$ and a similar fraction of candidates compared to PL1 + cPL2 + γ_ext, the latter is simpler in terms of model complexity (and has fewer free parameters) and is therefore favored by the BIC. In addition to this model ranking based on Bayesian statistics, we checked the physical properties, such as the mass-to-light ratios of the lens galaxies and the mass of the mass clumps in the PL1 + PL2+ rNIS model class in the following two subsections in order to show which models are physically reasonable.

5.4. Plausibility of inferred parameter values

To check the plausibility of the model parameter values of our two favored model classes, we estimate the mass-to-light ratio $\frac{M}{L}$ of G1 and G2 for both classes. We obtain the luminosity L by numerically integrating the modeled Sérsic light profile (see Eq. (1)) within the effective radius r_eff using a circular aperture. From the integrated fluxes, we obtain apparent magnitudes of m_{G1, AB, i} = 20.39 and m_{G2, AB, i} = 20.79 in the observed i-band. The two lensing galaxies are located at redshifts of 0.523 and 0.526, respectively, which approximately corresponds to rest-frame g-band. We conducted a K correction (Hogg et al. 2002) to transform the apparent magnitudes in the i-band to apparent magnitudes in the rest-frame g-band via m_AB, g = m_AB, i − K_gi with K_gi = −0.64. We used the galaxy spectra from Eigenbrod et al. (2006) and fitted an ETG spectral template. The rest-frame apparent magnitudes are therefore m_G1, AB = 21.03 and m_G2, AB = 21.43, which convert to L_G1 = 3.83 × 10¹⁰ L_⊙ and L_G2 = 2.69 × 10¹⁰ L_⊙. We infer the mass by integrating the κ_PL profile of the respective model class within the same effective radius as the light. From a distribution of candidate models, we find $M_{\mathrm{G}1,\mathrm{cPL}2}=8.1_{-0.7}^{+0.9}\times {10}^{10}{M}_{\odot }$ and $M_{\mathrm{G}2,\mathrm{cPL}2}=3.5_{-0.5}^{+0.6}\times {10}^{10}{M}_{\odot }$ for the PL1+cPL2+γ_ext model class, and $M_{\mathrm{G}1,\mathrm{rNIS}}=7.5_{-0.5}^{+0.6}\times {10}^{10}{M}_{\odot }$ and $M_{\mathrm{G}2,\mathrm{rNIS}}=2.5_{-0.7}^{+2.0}\times {10}^{10}{M}_{\odot }$ for the PL1+PL2+rNIS model class. For G1, we therefore find mass-to-light ratios of $(\frac{M}{L}{)}_{\mathrm{G}1,\mathrm{cPL}2}=2.1_{-0.2}^{+0.2}$ and $(\frac{M}{L}{)}_{\mathrm{G}1,\mathrm{rNIS}}=2.0_{-0.1}^{+0.2}$. For G2, we find $(\frac{M}{L}{)}_{\mathrm{G}2,\mathrm{cPL}2}=1.3_{-0.2}^{+0.2}$ and $(\frac{M}{L}{)}_{\mathrm{G}2,\mathrm{rNIS}}=0.{94}_{-0.2}^{+0.7}$. Galaxies with a similar magnitude are expected to have ${(\frac{M}{L})}_g\sim 1-5$ (Kauffmann et al. 2003; Bell et al. 2003), and therefore our inferred masses are plausible. We note that the mass of G2 is degenerate with the mass of the NIS subhalo, which means that M_G2, rNIS, and therefore the mass-to-light ratio, might be underestimated.

5.5. Detection of DM substructure?

The model class PL1+PL2+rNIS, which is preferred by the BIC, implies the presence of a mass clump about 0.5″ to 1″ northwest of G2. We subsequently investigated the nature of the predicted rNIS clumps in the PL1 + PL2 + rNIS model class by comparing their mass with the expected mass range for subhalos from simulations. For a total mass estimate of the rNIS clump, we calculated its effective subhalo lensing mass (Minor et al. 2017), which is robust to changes in the mass density profile and is a good estimate of the total mass of the subhalo. The effective subhalo lensing mass is the mass enclosed within a subhalo perturbation radius, scaled by the slope of the host galaxy, which we assume to be G2.

We show the distribution of M_rNIS of all the candidate models in the chain on a logarithmic scale in Fig. 7. The distribution spans a mass range of ∼8 × 10⁴ − 6 × 10¹⁶ M_⊙. Approximately 2% of all candidate models predict 10⁶ M_⊙ ≤ M_rNIS ≤ 10⁹ M_⊙, which is the range of predicted masses for DM substructure or dark dwarf-galaxies in simulations (Dalal & Kochanek 2001; Simon & Geha 2007; Diemand et al. 2008; Vegetti et al. 2012). Furthermore, ∼14% of all candidate models predict M_rNIS ≤ 10¹⁰ M_⊙. For halos with mass above 10⁹ M_⊙, we would expect them to host a galaxy and therefore be luminous. The large fraction of rNIS models with M_rNIS > 10¹⁰ results from the combination of a small Einstein radius of G2 and a large distance of the subhalo to G2. This leads to a bigger subhalo perturbation radius and therefore to a higher effective subhalo lensing mass. We conclude that the presence of unobserved dark matter structure close to the lens is unlikely to be responsible for the missing image because the effective subhalo masses are mostly too large compared to the expected masses for subhalos from simulations. It is possible that the rNIS tries to fit to the group halo of the galaxy group that G1 and G2 are part of, but the lack of high-resolution imaging data prevents us from further investigating this point at present. Our preferred explanation is therefore the presence of a core in G2, which changes the landscape of critical curves and caustics, and ultimately eliminates the fifth quasar image.

Fig. 7. Logarithmic distribution of the masses of rNIS clumps (in solar masses) in the candidate models of the PL1 + PL2 + rNIS + γext model class.

6. Summary

We modeled and analyzed the strongly lensed quasar HE0230−2130 using ground-based imaging data from the Magellan telescope. This lensing system shows a peculiar configuration with two lensing galaxies at similar redshift and four quasar images. The lensing galaxies are possibly part of a larger galaxy group and it is possible that the two lensing galaxies are dynamically associated or are embedded in an overall DM halo. As we expect five quasar images for such a lensing configuration, our aim is to find models that produce only four images. After we obtained the image positions by modeling the lens and quasar light in the observed image, we use them to constrain the mass parameters of our models. The straightforward approach of modeling both lens galaxies with a singular PL profile agrees with our expectations and always results in one additional predicted image, which is not observable in the data. We test a dozen different model parametrizations where we add core radii to the PL profiles describing G1 and G2, and/or we add a singular or cored isothermal sphere profile at the same redshift as G1 and G2, mimicking a dark mass distribution. We find that introducing a core to G2, or indeed to G1 and G2, can explain the missing image, as can a cored or singular isothermal sphere close to the predicted position of the missing image. In all of these cases, the saddle of the missing image merges with the maximum of G2, such that only four observable quasar images remain. All model classes consistently predict an external shear strength between 0.06 and 0.08. Adding external shear to the models generally increases the candidate fraction. This could mean that this system is potentially not dynamically relaxed, as the external shear parameter has recently been shown to compensate for missing model complexity. We compared the likelihoods of each candidate model class with the BIC, taking into account the fraction of candidate models. The most likely model class is cPL1 + cPL2 + rNIS, where 1.7% of all models produce the correct number of images, but the median mass of these NIS clumps is ∼2.9 × 10¹⁰ M_⊙, which is one order of magnitude larger than the upper mass limit of the predicted subhalos from simulations. Only about 2% of the candidate models have clumps in the predicted mass range of 10⁶ M_⊙ ≤ M_rNIS ≤ 10⁹ M_⊙ from simulations. Most of the clumps are too massive to be interpreted as undetected dark matter substructure. We therefore favor the explanation that a core in the mass distribution of G2 is responsible for the missing image.

We conclude that we find models that can describe this peculiar lensing system. Our favored explanation is the presence of a core in G2 in combination with external shear. Follow-up high-resolution and high-signal-to-noise-ratio imaging observations, which could display the quasar host galaxy as arcs or give us a better view of the region around G2, would help to further constrain our models and shed light on the unusual nature of HE0230−2130.

Acknowledgments

We thank R. Levinson for reduction of the imaging data, and S. Burles for the color image. We are grateful to the referee C. D. Fassnacht for constructive and detailed comments that improved the presentation and clarity of our work. S.E. and S.H.S. thank the Max Planck Society for support through the Max Planck Research Group and the Max Planck Fellowship for S.H.S. This work is supported in part by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. S.S. acknowledges financial support through grants PRIN-MIUR 2017WSCC32 and 2020SKSTHZ. This work made use of NumPy (Oliphant 2015), SciPy (Jones, et al. 2001), astropy (Astropy Collaboration 2018), and matplotlib (Hunter 2007).

References

Appendix A: Final model parameter values of all candidate model classes

Appendix B: Corner plots for candidate model classes

In this Appendix, we present constraints on a selection of mass and external shear parameters for each candidate model class. The distributions of candidates are plotted in red over the general distribution of the complete chain, which is colored in black.

Fig. B.1. PL1 + cPL2 + γext: posterior distributions for several lens-mass parameters of the power-law profile (see Sect. 4) and the external shear strength. We show the distributions for candidate models (i.e., models that predict four observable quasar images) in red, and the distribution of the whole MCMC chain is plotted in black. The two contours show the 1σ and 2σ credible regions. The one-dimensional histograms show the marginalized posterior distribution for the selected mass parameters, and the vertical lines mark the 1σ confidence intervals.

Fig. B.2. Posterior distributions for the cPL1 + cPL2 + γext model class.

Fig. B.3. Posterior distributions for the PL1 + PL2 + rSIS model class.

Fig. B.4. Posterior distributions for the PL1 + PL2 + rNIS model class.

Fig. B.5. Posterior distributions for the PL1 + PL2 + rSIS + γext model class.

Fig. B.6. Posterior distributions for the PL1 + PL2 + rNIS + γext model class.

All Tables

All Figures

	Fig. 1. Color image of HE0230−2130 created with Magellan imaging data in the g-, r-, and i-band. Image courtesy of Scott Burles.
In the text

	Fig. 2. Distributions of predicted image magnifications for the observed image D (red, solid line), the model-predicted image E (blue, solid line), and the demagnified image E in a worst-case microlensing scenario (blue, dashed line). Even if image E were strongly demagnified by microlensing, it would still be observable.
In the text

	Fig. 3. Modeled image of HE0230−2130. Plotted are the 1σ and 2σ contours of the predicted positions of the unobserved image E in the PL1 + PL2 (red), the PL1+cPL2+γext (green), and the PL1+PL2+rNIS (blue) model class. These predicted fifth quasar images would be above the detection threshold, but are not observed in the imaging data.
In the text

	Fig. 4. Observed and modeled surface-brightness distribution of HE0230−2130 and normalized residuals in a range of −5σ to 5σ.
In the text

Fig. 5. Change of critical curves and caustics between noncandidate models (left and middle columns) and one candidate model (right column). The top and middle rows show the critical curves (blue, solid lines) and caustic curves (red, dashed lines). The bottom row shows the Fermat potential and contours with the color scale in units of arcsec2. As the core radius of G2 increases and the source position crosses the radial caustic into the shaded green region, the time-delay saddle point associated with image E merges with the maximum of G2 and image E disappears. The source position is plotted as a yellow star. The predicted image positions are shown as green crosses in each panel.

In the text

	Fig. 6. 1σ and 2σ contours of the mass centroid of the candidate rNIS clumps in the PL1 + PL2 + rNIS model class. These candidate models predict only four quasar images, or additional images that are too faint to be observed (i.e., below the flux limit of the data).
In the text

	Fig. 7. Logarithmic distribution of the masses of rNIS clumps (in solar masses) in the candidate models of the PL1 + PL2 + rNIS + γext model class.
In the text

Fig. B.1. PL1 + cPL2 + γext: posterior distributions for several lens-mass parameters of the power-law profile (see Sect. 4) and the external shear strength. We show the distributions for candidate models (i.e., models that predict four observable quasar images) in red, and the distribution of the whole MCMC chain is plotted in black. The two contours show the 1σ and 2σ credible regions. The one-dimensional histograms show the marginalized posterior distribution for the selected mass parameters, and the vertical lines mark the 1σ confidence intervals.

In the text

	Fig. B.2. Posterior distributions for the cPL1 + cPL2 + γext model class.
In the text

	Fig. B.3. Posterior distributions for the PL1 + PL2 + rSIS model class.
In the text

	Fig. B.4. Posterior distributions for the PL1 + PL2 + rNIS model class.
In the text

	Fig. B.5. Posterior distributions for the PL1 + PL2 + rSIS + γext model class.
In the text

	Fig. B.6. Posterior distributions for the PL1 + PL2 + rNIS + γext model class.
In the text