© The Authors 2024

1. Introduction

Strong gravitational lensing is a formidable tool for the study of galaxy structure and cosmology. Time-delay measurements on strongly lensed variable sources, such as quasars or supernovae, can be used to constrain the Hubble-Lemaître constant, H₀, independently of local distance calibrations (Refsdal 1964; Suyu et al. 2017; Millon et al. 2020; Birrer et al. 2020; Treu et al. 2022). As the precision of time-delay analyses improves (Shajib et al. 2023; Acebron et al. 2023; Kelly et al. 2023), it becomes increasingly important to keep systematic uncertainties under control. The accuracy of a time-delay lensing inference depends critically on how well the mass distribution of the light-deflecting system is characterised, including both the main lens and the mass along the line of sight (Schneider & Sluse 2013; Xu et al. 2016; Sonnenfeld 2018; Kochanek 2020). While great effort has been put into developing methods for constraining both these components (Greene et al. 2013; Shajib et al. 2018; Yıldırım et al. 2023; Wells et al. 2023; Park et al. 2023), it is important to make sure that no potential source of error is left unchecked (Kochanek 2021; Van de Vyvere et al. 2022; Gomer et al. 2022, 2023).

Recently, Baldwin & Schechter (2024) pointed out the existence of a selection effect that might bias inferences of H₀ from time-delay lensing. Baldwin & Schechter (2024) focused on quadruply imaged point sources (quads), which make up the majority of systems used so far for time-delay analyses. In order for a source to be lensed into a quad, it is necessary for it to lie in the region of the source plane enclosed within the tangential caustic of the lens. Baldwin & Schechter (2024) argue that, since galaxies with an intrinsically larger tangential caustic area are more likely to be strong lenses, the selection of quads introduces a bias, as a result of which the estimated caustic area is systematically smaller than the true one. This bias was said to be similar in nature to the Malmquist bias, in which the observed brightnesses of a flux-limited sample of stars or galaxies are systematically overestimated near the detection limit. If present, such a bias would translate into a corresponding bias on H₀.

Here we build upon the work of Baldwin & Schechter (2024) to understand more generally how uncertainties on the estimated caustic area impact the inference of other lens properties. First, we derived an expression for the caustic area bias of a lens survey and applied it to a few example cases. We show that, due to an Eddington bias effect, the sign of the bias at the high mass end is opposite to that due to the Malmquist bias alone. Second, we focus on the lens modelling stage and show how a seemingly innocuous assumption can introduce a positive bias to the caustic area measurement. Along with the examples of biases, we suggest strategies for removing them.

This paper is organised as follows. In Sect. 2 we briefly introduce basic notions about strong lensing. In Sect. 3 we introduce the problem of the caustic area bias and derive an expression for its quantification. In Sect. 4 we show an instance of a caustic area bias introduced when modelling an individual lens and explain how to eliminate it. In Sect. 5 we discuss our results and draw conclusions.

2. Strong lensing basics

A strong lens can be fully described by specifying its dimensionless surface mass density, κ, as a function of position in the sky, θ,

$$\begin{array}{c}\hfill \kappa (\mathit{\theta })\equiv \frac{\mathrm{\Sigma }(\theta )}{{\mathrm{\Sigma }}_{\mathrm{cr}}},\end{array}$$(1)

where Σ is the physical projected mass density and Σ_cr is the critical surface mass density for lensing, determined by the redshift of the lens and of the background source. Given κ, we can compute the deflection angle as a function of position:

$$\begin{array}{c}\hfill \mathit{\alpha }(\mathit{\theta })=\frac{1}{\pi }{\displaystyle {\int }_{{\mathbb{R}}^2}}{d}^2\mathit{\theta }\prime \kappa (\mathit{\theta }\prime )\frac{\mathit{\theta }-\mathit{\theta }\prime }{{|\mathit{\theta }-\mathit{\theta }\prime |}^2}·\end{array}$$(2)

Images of a source at angular position β form at the solutions of the lens equation

$$\begin{array}{c}\hfill \mathit{\beta }=\mathit{\theta }-\mathit{\alpha }(\mathit{\theta }).\end{array}$$(3)

The magnification of lensed images can be obtained by means of the Jacobian matrix of the lens equation, A ≡ ∂β/∂θ. Lenses can have up to two curves where the A matrix is singular, one for each eigenvalue. These are called critical curves and correspond to points of formally infinite magnification. In our work we exclusively considered lenses with elliptically symmetric mass distributions. For such lenses, the two critical curves are labelled tangential and radial, depending on whether the corresponding magnification is largest in the direction tangential or perpendicular to the isodensity curves. The curve obtained by mapping a critical curve to the source plane with the lens equation is called caustic.

Throughout this work we focus on a particular class of lenses, with elliptically symmetric projected mass distribution and a power-law radial density profile. Their dimensionless surface mass density can be described by

$$\begin{array}{c}\hfill \kappa ({\theta }_{\mathrm{circ}})=\frac{3-\gamma }{2}{\left(\frac{{\theta }_{\mathrm{circ}}}{{\theta }_{\mathrm{Ein}}}\right)}^{1-\gamma },\end{array}$$(4)

where θ_circ is the circularised radius,

$$\begin{array}{c}\hfill {\theta }_{\mathrm{circ}}\equiv \sqrt{q{\theta }_1^2+\frac{{\theta }_2^2}{q}},\end{array}$$(5)

q is the axis ratio, θ_Ein is the Einstein radius, and γ is the corresponding three-dimensional density slope.

Figure 1 shows caustics of elliptical power-law lenses, with different values of q and γ. Starting from the γ = 1.8 case (left-hand panel), the outer, elliptical-shaped, curve is the radial caustic, while the diamond-shaped one is the tangential caustic. The position of the source relative to the caustics determines the number of images that are produced. If the source is inside both caustics, then five images are formed. One of these images is, in all practical cases, extremely de-magnified and hence ignored: the end result is a quad lens. At fixed lens properties, the source-plane area that is mapped to a quad configuration generally consists of the region enclosed within the tangential caustic. The only exception is when the tangential and radial caustic intercept each other, such as in the q = 0.6, γ = 1.8 example. In those cases, only the region enclosed by both caustics gives rise to sets of four images.

Fig. 1. Caustics of power-law lenses with elliptical mass distributions. Each panel corresponds to a value of the density slope, γ. Different colours correspond to different values of the axis ratio, as indicated in the legend. The outer, elliptical-shaped curve is the radial caustic (in the γ = 1.8 case) or the cut (in the γ = 2 case), while the inner diamond-shaped one is the tangential caustic. The major axis of the lens points in the horizontal direction. Caustics were obtained with the software GLAFIC (Oguri 2021).

Power-law ellipsoid lenses with γ ≥ 2 have no radial critical curve, and therefore no radial caustic. In the limit γ → 2, the radial caustic turns into a cut, which is a curve that separates regions of the source plane that are mapped into one or two images. Differently from a caustic, sources that lie exactly on the cut do not give rise to infinitely magnified images. The cut is visible in Fig. 1 for the γ = 2.0 case (middle panel), while is outside of the plotted range in the γ = 2.2 case (right-hand panel). Sources that lie in the region between the tangential caustic and the cut are mapped into two images, while those inside the tangential caustic are mapped into four images.

Figure 1 shows a clear trend between the axis ratio of the lens mass distribution and the area enclosed by the tangential caustic: the higher the ellipticity, the larger the area. Although more difficult to see, the tangential caustic area also decreases with increasing γ, at fixed q. A larger or smaller tangential caustic area generally translates into a higher or lower probability of a lens creating a quad configuration. However, the cross-section for a quad event depends also on the intrinsic brightness of the source and the magnification of the multiple images, because the sole fact that four images are produced does not guarantee that they are detectable. For the sake of simplicity, we ignore magnification throughout this work. This assumption corresponds to the limit of a very bright source, which is a reasonable one in the case of lensed quasars. In Sect. 5 we discuss how the findings of this paper can be generalised to the case of fainter sources.

The mass distribution of a lens can be inferred by fitting a model to the observed lensing data, such as the positions of the multiple images. This is a fundamental step in all time-delay lensing analyses. However, there is a limit to how much information can be extracted in the process. Given a lens with surface mass density distribution κ(θ), a transformation of the kind

$$\begin{array}{cc}& \kappa (\mathit{\theta })\to {\lambda }_{\mathrm{mst}}\kappa (\mathit{\theta })+1-{\lambda }_{\mathrm{mst}},\hfill \end{array}$$(6)

$$\begin{array}{cc}& \mathit{\beta }\to {\lambda }_{\mathrm{mst}}\mathit{\beta }\hfill \end{array}$$(7)

leaves image positions and flux ratios unchanged but varies the time delay between images. The original and the transformed lens are therefore indistinguishable on the basis of the data typically used to constrain lens models, but different values of H₀ are needed to match a given time delay.

Equation (7) correspond to a uniform rescaling of the source plane coordinates. This means that the area enclosed within any caustic is a factor of ${\lambda }_{\text{mst}}^2$ of the original one. In the bright source limit, the probability of creating a quad event is then also rescaled by the same factor. The pair of Eqs. (6) and (7) is referred to as a mass-sheet transformation, and the invariance of lensing observables under it is called mass-sheet degeneracy (Falco et al. 1985). As we discuss in the next section, it plays an important role in the caustic area bias. For a more thorough introduction to strong lensing we refer the reader to Meneghetti (2021) and references within.

3. Population-level bias

We considered a population of foreground galaxies, all at the same redshift, z_g, and a population of bright background point sources, also at a single redshift, z_s > z_g. Each lens has an associated caustic area for the formation of quads, A_quad. The distribution in quad caustic area of the population of foreground galaxies is described by a probability distribution P_g(A_quad). Since the probability of a background source being lensed into a quad by a given foreground galaxy is proportional to A_quad, the distribution in A_quad of the population of quad lenses is

$$\begin{array}{c}\hfill {\mathrm{P}}_{\mathrm{quad}}(A_{\mathrm{quad}})\propto A_{\mathrm{quad}}{\mathrm{P}}_{\mathrm{g}}(A_{\mathrm{quad}}).\end{array}$$(8)

Next we considered observations. Although the caustic area is not a quantity that is usually explicitly fitted for, it is an implicit byproduct of any lens model: given a surface mass density distribution inferred from lensing data, it is always possible to calculate the corresponding value of A_quad. We label the observed A_quad obtained via lens modelling as $A_{\mathrm{quad}}^{(\mathrm{obs})}$. In general, if the lens model is not a perfect representation of the truth, $A_{\mathrm{quad}}^{(\mathrm{obs})}$ will be different from A_quad. Vice versa, any difference between the inferred and the true value of the quad caustic area is a sign that the model is inaccurate. With this definition, we introduced the likelihood of measuring a value $A_{\mathrm{quad}}^{(\mathrm{obs})}$ given the true value, $\mathrm{P}(A_{\mathrm{quad}}^{(\mathrm{obs})}|A_{\mathrm{quad}})$. The form of this distribution depends on the details of how $A_{\mathrm{quad}}^{(\mathrm{obs})}$ is obtained and on the noise properties of the data.

At fixed observed $A_{\mathrm{quad}}^{(\mathrm{obs})}$, the expectation value of the true quad caustic area is simply

$$\begin{array}{c}\hfill \mathrm{E}[A_{\mathrm{quad}}]=\frac{\int d{A}_{\mathrm{quad}}{A}_{\mathrm{quad}}{\mathrm{P}}_{\mathrm{quad}}(A_{\mathrm{quad}})\mathrm{P}(A_{\mathrm{quad}}^{(\mathrm{obs})}|A_{\mathrm{quad}})}{\int d{A}_{\mathrm{quad}}{\mathrm{P}}_{\mathrm{quad}}(A_{\mathrm{quad}})\mathrm{P}(A_{\mathrm{quad}}^{(\mathrm{obs})}|A_{\mathrm{quad}})}·\end{array}$$(9)

Any difference between this quantity and $A_{\mathrm{quad}}^{(\mathrm{obs})}$ is a bias on the caustic area. We label this the population-level bias, to distinguish it from a separate effect studied in Sect. 4. If the likelihood is symmetric around the true value, then what determines the sign of the bias is the shape of the distribution in quad caustic area of the lenses, P_quad(A_quad). In the hypothetical, and unrealistic, case of a uniform distribution in the quad caustic area of the foreground galaxies, P_g(A_quad)∼constant, then P_quad(A_quad)∝A_quad and E[A_quad] is larger than the observed value: the quad caustic area is underestimated. This is the scenario implicitly studied by Baldwin & Schechter (2024), who did not consider the distribution in caustic area in their treatment. In this scenario, the root of the bias is uniquely in the presence of an observational scatter and the preferential selection of lenses with larger caustic area, an effect referred to as a Malmquist-type bias. In general, though, P_g(A_quad) is a steeply declining function of A_quad. As a result, the bias can have the opposite sign.

In order to gain intuition on the problem, we carried out the following experiment. We considered a population of singular isothermal ellipsoid (SIE) lenses. These are elliptical power-law mass distributions with a value of the slope γ = 2. We fixed the value of the axis ratio to q = 0.7, for simplicity. With this choice, the quad caustic area is

$$\begin{array}{c}\hfill A_{\mathrm{quad}}^{(\mathrm{SIE})}(q=0.7)=0.066{\left(\frac{{\theta }_{\mathrm{Ein}}}{1″}\right)}^2{\mathrm{arcsec}}^2.\end{array}$$(10)

We then used the following relation, valid for a singular isothermal sphere, to link the Einstein radius to the velocity dispersion, σ_v, of the lens:

$$\begin{array}{c}\hfill {\theta }_{\mathrm{Ein}}=\frac{{\sigma }_v^2}{G{D}_{\mathrm{d}}{\mathrm{\Sigma }}_{\mathrm{cr}}},\end{array}$$(11)

where D_d is the angular diameter distance to the lens. We set z_g = 0.3, z_s = 1.5, and assumed a flat Λ cold dark matter cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_M = 0.3, from which follows

$$\begin{array}{c}\hfill {\theta }_{\mathrm{Ein}}\simeq 0.83″{\left(\frac{{\sigma }_v}{200\mathrm{km}{\mathrm{s}}^{-1}}\right)}^2.\end{array}$$(12)

Finally, we assumed that the values of σ_v of the galaxy population are drawn from a distribution with the following functional form

$$\begin{array}{c}\hfill \mathrm{P}({\sigma }_v)\propto {\left(\frac{{\sigma }_v}{{\sigma }_{\ast }}\right)}^{\alpha -1}exp\{-{\left(\frac{{\sigma }_v}{{\sigma }_{\ast }}\right)}^{\beta }\},\end{array}$$(13)

which is often used to describe the velocity dispersion function of galaxies (Sheth et al. 2003). Based on the measurements of Montero-Dorta et al. (2017), we set α = 6.0, β = 2.4 and σ_* = 130 km s⁻¹. Given Eq. (13), we can compute the probability distribution in quad caustic area of the foreground galaxies, P_g(A_quad). From Eqs. (10) and (11) it follows that $A_{\text{quad}}\propto {\sigma }_v^4$. Therefore,

$$\begin{array}{c}\hfill {\mathrm{P}}_{\mathrm{g}}(A_{\mathrm{quad}})\propto A_{\mathrm{quad}}^{-3/4}\mathrm{P}\left({\sigma }_v(A_{\mathrm{quad}})\right)\end{array}$$(14)

and

$$\begin{array}{c}\hfill {\mathrm{P}}_{\mathrm{quad}}(A_{\mathrm{quad}})\propto A_{\mathrm{quad}}^{1/4}\mathrm{P}\left({\sigma }_v(A_{\mathrm{quad}})\right).\end{array}$$(15)

Defining A_* = A_quad(σ_v = σ_*), we have

$$\begin{array}{c}\hfill {\mathrm{P}}_{\mathrm{quad}}(A_{\mathrm{quad}})\propto A_{\mathrm{quad}}^{(1+\alpha )/4}exp\{-{\left(\frac{A_{\mathrm{quad}}}{A_{\ast }}\right)}^{\beta /4}\}.\end{array}$$(16)

Figure 2 (left-hand panel) shows the distribution in the quad caustic area of the foreground galaxy population and that of the lens population. The right-hand panel shows the corresponding distributions in Einstein radius. These are qualitatively similar to that of existing samples of galaxy-scale lenses (see e.g. Auger et al. 2009; Sonnenfeld et al. 2013).

Fig. 2. Simulation of foreground galaxies (blue curves) and quad lenses (red curves). Left: distribution of the quad caustic area. Right: distribution of the Einstein radius. All foreground galaxies are SIEs with axis ratio q = 0.7.

The last ingredient needed to compute the bias is the likelihood function, $\mathrm{P}(A_{\mathrm{quad}}^{(\mathrm{obs})}|A_{\mathrm{quad}})$. We assumed a Gaussian form,

$$\begin{array}{c}\hfill \mathrm{P}\left(A_{\mathrm{quad}}^{(\mathrm{obs})}|A_{\mathrm{quad}}\right)=\frac{1}{\sqrt{2\pi }{\sigma }_A}exp\{-\frac{{(A_{\mathrm{quad}}^{(\mathrm{obs})}-A_{\mathrm{quad}})}^2}{2{\sigma }_A^2}\},\end{array}$$(17)

centred on the true value and with dispersion σ_A. We set σ_A to be a fixed fraction s_A of the true value of the quad caustic area,

$$\begin{array}{c}\hfill {\sigma }_A=s_A{A}_{\mathrm{quad}}.\end{array}$$(18)

We explored three different values of s_A: 0.05, 0.10, and 0.20.

Figure 3 shows the relative bias on the quad caustic area, computed using Eq. (9), as a function of the observed value. For small $A_{\mathrm{quad}}^{(\mathrm{obs})}$ the bias is negative: the caustic area is underestimated, as predicted by Baldwin & Schechter (2024). This is a consequence of the fact that, for a sufficiently large value of the low−σ_v slope α, P_quad(A_quad) is an increasing function of A_quad: the observational scatter scrambles the observed values of the caustic area, and there are more lenses with intrinsically larger A_quad that scatter low than the other way around. This effect is known as Eddington bias. However, at large $A_{\mathrm{quad}}^{(\mathrm{obs})}$ the trend is the opposite, because P_quad(A_quad) declines with increasing A_quad. The transition between negative and positive trend occurs at a larger value of $A_{\mathrm{quad}}^{(\mathrm{obs})}$ than the peak in the P_quad(A_quad) distribution simply because there is more volume in the parameter space at large A_quad, giving more weight to larger values in the integrand of Eq. (9).

Fig. 3. Relative bias on the quad caustic area in our simulation, for three different values of sA, as defined in Eq. (18). The vertical dotted line marks the quad caustic area of an SIE lens with q = 0.7 and an Einstein radius of 1.5″, for reference.

Our experiment can also give us an idea of the amplitude of the caustic area bias that we can expect in practice. The key quantity for such an estimate is s_A, which describes the fractional scatter in the caustic area measured from lens modelling with respect to the true value. In the case of time-delay lensing analyses, the dominant source of scatter is the external convergence from line-of-sight matter, κ_ext. To first approximation, the addition of external convergence to a lens model acts like a mass-sheet transformation that re-scales source plane areas by a factor of (1 − κ_ext)². Typical values of κ_ext are ≈0.05 (Wells et al. 2023), which translate to s_A ≈ 0.1. According to our simulation, then, the caustic area bias has a typical amplitude of 1 − 2%. By the same argument, a 2% bias on A_quad is equivalent to a 1% bias on κ_ext. The effect of κ_ext is that of rescaling time delays by 1 − κ_ext; therefore, its bias propagates directly to the inferred value of H₀. A 1% bias on H₀ is smaller than the current statistical uncertainty from time-delay analyses. Moreover, as we discuss in Sect. 5, this is likely an overestimate due to our bright source limit assumption.

We now go back to considering the general case. We argue that, if the likelihood function is known, then the caustic area bias can be accounted for in a lensing study. The problem of determining the true value of A_quad of a lens is no different from that of determining any noisy property of a galaxy. In both cases, the key is to have an accurate prior on the underlying distribution of the parent population of objects and an accurate knowledge of the likelihood function. For example, assuming a uniform prior on A_quad leads to a bias of the kind shown in Fig. 3, just like assuming a uniform prior on the stellar mass of a galaxy leads to a bias originating from the non-uniformity of the galaxy stellar mass function (i.e. the Eddington bias). One way of removing the bias is by modelling the underlying population of objects. In the lensing case, this can be done by fitting a hierarchical model to a large sample of lenses (see e.g. Sonnenfeld et al. 2015; Shajib et al. 2021). Most importantly, although the lensing selection determines the distribution in the lens properties P_quad, it is perfectly possible to infer the latter without explicitly accounting for the former. As we show in the next section, however, selection effects can play a role at the lens modelling stage.

4. Bias from modelling

In this section we describe a simple scenario in which a seemingly straightforward lens modelling approach leads to a biased estimate of the lens properties. We show that this is a kind of caustic area bias, originating from an improper treatment of the lens selection procedure, and give a prescription for how to remove it.

We took an SIE lens with Einstein radius θ_Ein = 1.0″ and axis ratio q = 0.7. We placed a background point source located within the tangential caustic, hence producing a quad configuration. Then, we fitted a power-law ellipsoid to the four image positions of the lensed source. We assumed the centroid of the lens to be known exactly and we optimised for the following parameters: the Einstein radius θ_Ein, the density slope γ, the axis ratio q, the position angle (PA) of the major axis, and the source position, β. We assumed a uniform prior within generous fixed bounds for all of the lens parameters. We assigned a 0.005″ uncertainty to the observed image positions, then sampled the posterior probability distribution of the model parameters with a Markov chain Monte Carlo (MCMC). We used the software GLAFIC (Oguri 2021) to solve the lens equation for each set of model parameters and EMCEE (Foreman-Mackey et al. 2013) to run the MCMC. The number of constraints (eight) exceeds the number of free parameters (six), and therefore the problem is well posed. The power-law ellipsoid model family used for the fit includes the true values of the lens parameters. Since we did not add noise to the data, one would expect that the fit recovered the input. However, as Fig. 4 shows (blue contours), the posterior probability of the individual lens parameters is biased with respect to the truth. For instance, the median value of the marginalised posterior on the density slope is γ = 2.06, substantially larger than the input value of γ = 2.00.

Fig. 4. Posterior probability distribution of the power-law ellipsoid fit described in Sect. 4. Blue-filled contours represent the posterior probability obtained assuming a fixed uniform prior on the source position, and solid red contours the posterior probability of the model with the prior of Eq. (19). The quantities shown are: the Einstein radius (in arcsec), density slope, position angle of the major axis (in degrees), axis ratio, source position in polar coordinates (distance from the optical axis in arcsec and angle in degrees), and quad caustic area (in square arcsec). The dashed lines indicate the true values of the parameters. The median and 68% credible region of the marginal posterior in each quantity is listed at the top of each column.

The source of the bias is the following. The lens was selected to be a quad, but this piece of information was not used during the fit. To do it, we must assume the following prior on the source position, β:

$$\begin{array}{c}\hfill \mathrm{P}(\mathit{\beta })=\{\begin{array}{cc}\frac{1}{A_{\mathrm{quad}}}\hfill & \mathrm{if}\mathit{\beta }\mathrm{is}\mathrm{within}\mathrm{the}\mathrm{quad}\mathrm{region}\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}.\end{array}$$(19)

In words, the prior of Eq. (19) is constant over the region corresponding to a quad configuration (in most practical cases, the region enclosed within the tangential caustic) and zero outside of it: models with fewer than four images are not considered, because such lenses would not have been selected in the first place. Naively one might think that, since the likelihood of a model with the wrong number of images is zero anyway, there is no practical difference between the original fixed uniform prior and that of Eq. (19). However, the two differ under one crucial aspect: the normalisation. Neglecting it causes models corresponding to larger values of A_quad to be favoured. This can be clearly seen in Fig. 4: the quad caustic area A_quad correlates with all of the other lens model parameters, and the inference of these parameters is biased towards values that correspond to a larger A_quad. Figure 4 also shows the inference obtained with the correct prior of Eq. (19): in this case all of the input values of the parameters are accurately recovered. This bias looks qualitatively similar to one reported by Ruan & Keeton (2023), who performed a similar experiment to the one described here. It is possible that their bias can also be ascribed to the choice of the source position prior.

Rescaling the source position prior on the basis of the caustic area is not a standard practice in lens modelling. This raises the question of why such a large bias on the density slope of power-law models has not been reported before. The answer is in the quality of the data used for the fit. In our experiment, we only fitted image positions: they can allow us to measure the Einstein radius of a lens but have little constraining power on the density profile. Instead, state-of-the-art lens models are obtained by fitting the full surface brightness distribution of an extended source (in the case of a lensed quasar, the host galaxy). Therefore, we believe that the amplitude of bias is a function of the constraining power of the lensing data. We verified this conjecture by repeating the experiment while fitting flux ratios, in addition to the image positions: we did not detect any bias in that case. This is in agreement with results from tests of power-law ellipsoid fits to simulated lenses with smooth extended sources (Suyu 2012). Tests on more complex simulations, such as those of the time-delay lens modelling challenge (Ding et al. 2021), have revealed biases on the inferred lens parameters, though they are likely the result of incorrect modelling of the source structure or the point-spread function (Galan et al. 2024).

5. Discussion and conclusions

In Sect. 3 we argue that when the likelihood of the lensing data is known, the caustic area bias can be accounted for by modelling the population of lenses. In Sect. 4 we show that although the standard way of setting the source position prior can lead to a biased inference of the parameters of a lens, such a bias can be removed by accounting for the lens selection criteria. We also point out how this lens modelling bias disappears in the limit in which the model is well constrained by the data. Although this work deals exclusively with quad lenses, all findings can be generalised to lenses of any configuration.

The entire analysis was carried out in the limit of a very bright source. Nevertheless, we can easily drop that assumption. The first step is to introduce the source flux parameter, f_s, and its corresponding population distribution, P_s(f_s). Second, we introduce the strong lensing cross-section, as defined by (Sonnenfeld et al. 2023):

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{SL}}={\displaystyle {\int }_{{\mathbb{R}}^2}}d\mathit{\beta }{\mathrm{P}}_{\det }({\mathit{\psi }}_{\mathrm{g}},\beta ,f_{\mathrm{s}}|\mathrm{quad}).\end{array}$$(20)

In this equation, P_det is the probability of detecting a strong lens given the parameters describing the lens galaxy (ψ_g), the background source flux and position, and the selection criterion, which in our case is the detection of four images. Next, to repeat the analysis from Sect. 3, we need to marginalise σ_SL over the source flux, to obtain an effective cross-section averaged over the background source population:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{SL}}^{(\mathrm{eff})}\equiv {\displaystyle \int \mathrm{d}}{f}_{\mathrm{s}}{\sigma }_{\mathrm{SL}}{\mathrm{P}}_{\mathrm{s}}(f_{\mathrm{s}}).\end{array}$$(21)

At this point we can simply replace A_quad with ${\sigma }_{\mathrm{SL}}^{(\mathrm{eff})}$. Similarly, we can generalise the prior of Eq. (19) as follows:

$$\begin{array}{c}\hfill \mathrm{P}(\mathit{\beta },f_{\mathrm{s}})=\frac{1}{{\sigma }_{\mathrm{SL}}^{(\mathrm{eff})}}{\mathrm{P}}_{\det }({\mathit{\psi }}_{\mathrm{g}},\mathit{\beta },f_{\mathrm{s}}){\mathrm{P}}_{\mathrm{s}}(f_{\mathrm{s}}).\end{array}$$(22)

This formalism can allow us in principle to describe scenarios with arbitrarily complex selection criteria, which are encoded in P_det. The difficult part is obtaining an accurate description of these criteria on real lenses. Even assuming that P_det is known exactly, implementing a prior of the kind of Eq. (22) in a lens modelling code might not be easy: it would require, for each possible set of lens model parameters (ψ_g), the lensing cross-section to be computed for all possible values of the source flux.

One important consequence of dropping the bright source assumption is that the amplitude of the caustic area biases decreases. This was already pointed out by Baldwin & Schechter (2024) and follows from the fact that, in all practical cases, increasing the caustic area at a fixed Einstein radius causes the magnification to decrease. As a result, the variation in ${\sigma }_{\mathrm{SL}}^{(\mathrm{eff})}$ with lens model parameters is milder than in the bright source case, which means that the likelihood $\mathrm{P}(A_{\mathrm{quad}}^{(\mathrm{obs})}|A_{\mathrm{quad}})$ becomes narrower.

In light of our findings, a time-delay cosmography study that aims to remove all sources of caustic area bias must do the following: First, model individual lenses while ensuring that the lensing data used for this purpose have a high signal-to-noise ratio. The actual signal-to-noise ratio threshold needs to be determined via dedicated simulations and depends on the level of accuracy desired. The end product of the individual modelling phase will be lens models that are constrained to a high precision but are still potentially subject to the mass-sheet degeneracy, from line-of-sight structure or from a mismatch between the model family and the true density profile of the lenses (such as the power-law models used by Birrer et al. 2020). Second, obtain measurements of an ensemble of observables that correlates with the mass-sheet transformation parameter, λ_mst, at a fixed Einstein radius in a well-understood way (again, in relation to the desired level of accuracy). Such observables could be the number of galaxy counts around the main lens, which can be related to the external convergence, or the velocity dispersion profiles of the lenses, which constrain their internal structure. Third, fit a hierarchical model describing how the lens parameters, ψ_g, which now include λ_mst, are distributed among the population of lenses. State-of-the-art time-delay lensing analyses are already going in this direction: Birrer et al. (2020) included a hierarchical model of the lens population in their study of seven time-delay lenses, although only two of the lens model parameters were modelled hierarchically (λ_mst and the orbital anisotropy parameter). Future lens population models should aim to capture the distribution of all parameters related to the mass density profile of the lenses. By doing so, the Eddington bias on all quantities that correlate with the time delays, including but not limited to the caustic area, can be taken care of. In any case, given its very small amplitude, removing the caustic area bias is not a priority for current time-delay lensing analyses.

Acknowledgments

I thank Paul Schechter for his constructive feedback on this manuscript. This work was supported by the National Science Foundation of China, with an Excellent Young Scholars Fund Overseas grant.

References

All Figures

Fig. 1. Caustics of power-law lenses with elliptical mass distributions. Each panel corresponds to a value of the density slope, γ. Different colours correspond to different values of the axis ratio, as indicated in the legend. The outer, elliptical-shaped curve is the radial caustic (in the γ = 1.8 case) or the cut (in the γ = 2 case), while the inner diamond-shaped one is the tangential caustic. The major axis of the lens points in the horizontal direction. Caustics were obtained with the software GLAFIC (Oguri 2021).

In the text

	Fig. 2. Simulation of foreground galaxies (blue curves) and quad lenses (red curves). Left: distribution of the quad caustic area. Right: distribution of the Einstein radius. All foreground galaxies are SIEs with axis ratio q = 0.7.
In the text

	Fig. 3. Relative bias on the quad caustic area in our simulation, for three different values of sA, as defined in Eq. (18). The vertical dotted line marks the quad caustic area of an SIE lens with q = 0.7 and an Einstein radius of 1.5″, for reference.
In the text

Fig. 4. Posterior probability distribution of the power-law ellipsoid fit described in Sect. 4. Blue-filled contours represent the posterior probability obtained assuming a fixed uniform prior on the source position, and solid red contours the posterior probability of the model with the prior of Eq. (19). The quantities shown are: the Einstein radius (in arcsec), density slope, position angle of the major axis (in degrees), axis ratio, source position in polar coordinates (distance from the optical axis in arcsec and angle in degrees), and quad caustic area (in square arcsec). The dashed lines indicate the true values of the parameters. The median and 68% credible region of the marginal posterior in each quantity is listed at the top of each column.

In the text