A Model-Independent Characterisation of Strong Gravitational Lensing by Observables

Abstract

When light from a distant source object, like a galaxy or a supernova, travels towards us, it is deflected by massive objects that lie in its path. When the mass density of the deflecting object exceeds a certain threshold, multiple, highly distorted images of the source are observed. This strong gravitational lensing effect has so far been treated as a model-fitting problem. Using the observed multiple images as constraints yields a self-consistent model of the deflecting mass density and the source object. As several models meet the constraints equally well, we develop a lens characterisation that separates data-based information from model assumptions. The observed multiple images allow us to determine local properties of the deflecting mass distribution on any mass scale from one simple set of equations. Their solution is unique and free of model-dependent degeneracies. The reconstruction of source objects can be performed completely model-independently, enabling us to study galaxy evolution without a lens-model bias. Our approach reduces the lens and source description to its data-based evidence that all models agree upon, simplifies an automated treatment of large datasets, and allows for an extrapolation to a global description resembling model-based descriptions.

1. Introduction

Forty years ago, in 1979, two images of the quasar QSO 0957+561 were observed, [1], which marked the discovery of strong gravitational lenses. The first description of the image configuration in form of a lens model followed in [2]. In 1986, luminous arcs, highly magnified images of background galaxies, were discovered behind the galaxy clusters Abell 370 and Cl2244-02, [3,4]. The first lens models for Abell 370 were set up in [5,6]. Since then, observations of multiple images have been routinely used to probe the deflecting mass-density distribution, in particular in order to infer its dark matter content and dark matter properties—see, e.g., [7,8] for galaxy-scale lens surveys, or [9,10] for galaxy-cluster-scale lens surveys. Recent benchmarks with simulations of cluster-scale lenses show that a lot of different approaches accurately and precisely reconstruct the mass-density distribution in the vicinity of extended multiple images up to a few percent, [11]. For galaxy-scale lenses, an analogous project is currently being pursued, [12]. Furthermore, strong gravitational lensing of time varying objects has been used as a cosmological probe to determine the Hubble-Lemaître constant, $H_0$, [13,14,15]. Attempts to infer the cosmic spatial curvature, the cosmic matter density parameter, and dark energy properties of the current cosmological concordance model and its potential extensions are also being pursued with galaxy-scale and cluster-scale lenses, [16,17,18]. Even after forty years of strong gravitational lensing studies, it still is a subject of high research interest because unprecedented observations, like supernova (SN) Refsdal, [19], or the recently discovered fast radio bursts (FRBs), [20], have continuously contributed to develop lens reconstruction methods further and thereby broaden the range of strong-lensing applications.

In our research, as being developed in [21,22,23,24,25,26,27], we pursue the question which properties of a strong gravitational lens can be directly inferred from the observables that describe a multiple-image configuration. Contrary to most approaches, we do not aim for a global mass-density reconstruction of the entire lens that usually employs the observables as constraints in a model-fitting problem. For galaxy-scale lenses, such global lens reconstructions using parametric models or free-form methods can be found in [28,29,30,31,32,33], for instance. Galaxy-cluster mass densities can be reconstructed with global lens reconstructions as, for instance, detailed in [34,35,36,37]. References to further approaches and a comparison between the approaches can be found in [11]. Some of the approaches originally developed for galaxy-scale lenses can also be employed and extended to reconstruct cluster-scale lenses and vice versa, see [38] for an example. Due to the high symmetry of the lenses on galaxy scale, fully automatic characterisations are easier to implement than for cluster-scale lenses. The latter usually require at least a visual inspection of the results, or manual fine-tuning at intermediate stages. In addition, aiming at a global characterisation, the sparsely distributed multiple images are not sufficient to fully constrain the mass density within the cluster, [39]. As a consequence, different approaches employ different additional information and assumptions to reconstruct the mass-density distribution, which complicates the comparison of the results.

Instead of finding the global deflecting mass density that can cause the multiple images observed, we separate data-based information from model-based assumptions to find those lens properties that all approaches agree upon. We employ the most general equations of the standard gravitational lensing formalism to directly determine local lens properties from the observables. These local properties are unique and calculated in the same way for all lenses from symmetric galaxy-scale lenses to amorphous galaxy-cluster-scale lenses. Thus, the approach is highly efficient and robust to handle large datasets with a minimum amount of manual intervention.

The cosmic distances between us, the lens, and the source are usually determined according to the cosmological standard model and, as such, are not model-independent, see e.g., [40]. To make them independent of any assumption about the nature and the overall abundance of dark matter and dark energy in the universe, we set up data-based cosmic distances from the most recent Pantheon sample of SNe type Ia, [41]. This covers SNe up to redshift 2.3 such that data-based distances are available for a lot of observed lensing configurations. In this way, our approach is independent of a specific lens model and of a specific parametrisation of the cosmological background model, which is only assumed to be spatially homogeneous and isotropic and should fulfil the Einstein field equations. To set up the data-based cosmic distances, we employ the analytic orthonormal basis introduced in [42] for the sake of run-time efficiency and numerical stability. Other approaches to develop distance measures based on the standardisable candles are, for instance, [43,44,45,46].

In this work, we summarise the status quo of our approach, show a comparison to different global lens reconstruction ansatzes, based on [47,48], and detail its contribution to determine the invariance transformations of the most general equations of the gravitational lensing formalism, especially the time-delay equation which is used to infer $H_0$. Based on these results, we outline future developments of lens reconstructions that could be enabled and tested by the recently discovered FRBs, as further described in [49].

The next sections are organised as follows: In Section 2, we introduce the theory of the model-independent lens characterisation for the leading-order configurations of multiple images. We list the local lens characteristics that can be determined by the different observables in each multiple-image configuration. Section 3 shows examples for these different cases on galaxy and on galaxy-cluster scale and the results obtained by our approach. We also estimate to which precision $H_0$ can be determined from the time-delay equation employing the Pantheon-data-based angular diameter distances. After these proof-of-concept examples, we summarise the assets and applications of our approach in Section 4 and conclude the review with a diagrammatic summary of our method.

2. Materials and Methods

In this section, we introduce the theoretical concepts to determine local properties of a strong gravitational lens without assuming a specific lens model and without assuming a specific parametrisation of a spatially homogeneous and isotropic Friedmann background cosmology. Most cases of strong gravitational lensing are well described by an effective theory of light deflection by a mass distribution in the static, thin-screen, weak-field limit in a spatially homogeneous and isotropic background cosmology, as summarised in [50,51]. This means that the deflecting mass distribution is assumed to be located in one (or several) plane(s) orthogonal to the line of sight, the lens plane(s). The mass density in these lens planes can be described by a Newtonian gravitational potential on top of the cosmic background metric which is usually based on a Friedmann–Lemaître–Robertson–Walker metric (FLRW metric). Depending on the mass-density distribution, as investigated in [52], this Newtonian potential can deflect light rays emitted from a source, lying in a source plane, into multiple images. It is assumed that the distances on the sky between the angular positions of the source and the multiple images are small compared to the distances that the light rays travel along the line of sight. For the light deflection to be static, all motions of the observer, the source, and the lens with respect to each other are considered negligible during the observing time. If the source object has an intrinsic time-varying intensity, the deflected light rays arrive at different times at the observer because they cross environments of different mass densities and travel distances of varying length from the source to the observer. Figure 1 depicts an example of the strong gravitational lensing effect in which one background source is mapped to two images by a deflecting mass density in one lens plane.

In principle, there are other ways to describe the deflection of light rays in our universe which are detailed in [27]. For instance, the approach in which the background metric accounts for the entire deflecting mass-density distribution is an equivalent ansatz. However, determining the geodesics in such an inhomogeneous universe is more difficult than treating gravitational lensing in an effective theory based on lens planes in a homogeneous and isotropic background cosmology. Another option partitions the mass density along the line of sight into a main deflecting object and small-scale perturbations scattered along the light path in an FLRW background metric without projecting the masses to lens planes.

In order to avoid degeneracies within the description, we should choose the background metric, small-scale weakly deflecting masses, and deflecting masses that cause multiple images in such a way that the observables can directly constrain all parts of the chosen partition. We detail our choices of an FLRW background cosmology and a single-lens-plane Newtonian deflection potential in Section 2.1 and Section 2.2.

Given these prerequisites, singularity theory has been established as the mathematical foundation to characterise the generic, nonlinear mappings between multiple-image configurations and their sources, [53]. Originally developed by Vladimir Arnold to describe the set of singular points arising on manifolds that undergo (sudden) changes, e.g., projections to lower dimensions, singularity theory is the ideal tool to describe gravitational lensing in the lens-plane formalism introduced above. It thus provides a framework to solve the under-constrained inverse problem of reconstructing the gravitational lens and the source from an observed multiple-image configuration.

2.1. Observables in Multiple Images

As the (time-varying) brightness profile of multiple images caused by the strong gravitational lensing effect depends on the (time-varying) brightness profile of the source object and the deflecting mass distribution, one may expect a large variance in the appearances of multiple images and an equally high variance in the arrival time patterns of time-varying multiple images. However, the current measurement precision and resolution allow us to categorise the observations into three types of multiple images sharing the same observable features as summarised in

Table 1. Types of multiple images and quantitative observables that can be extracted from each image. Each type can occur in different multiple-image configurations shown in Table 2. Image credits: NASA/ESA/HST.

Point Image	Unresolved Image	Resolved Image
angular position on the sky	centre of light, quadrupole	vectors of distances between features

. The number of multiple-image configurations is also limited and most of the occurring cases can be categorised into the three configurations shown in

Table 2. Special multiple-image configurations of a common source. The figure shows an observed example, the schematic underneath the observables for unresolved images. Point-like and resolved images can also form these configurations. The bottom rows list the quantitative observables for each case. For arcs, we assume the centre of brightness of the lens is known. Image credits: NASA/ESA/HST.

Fold	Cusp	Giant Arcs
rel. distances ${\delta }_{ij}$ rel. orientations ${\phi }_i$ (time delays)	rel. distances ${\delta }_{ij}$ rel. orientations ${\tilde{\phi }}_i$ (time delays)	radial distances $r_i$, arc lengths $r_i{\phi }_i$ (time delays)

. Here, we focus on the observational properties of these configurations. In Section 2.2, we give an explanation for the limited amount of possible configurations based on singularity theory. Furthermore, general rules about the order in which light from these multiple-image configurations reaches the observer can also be derived from singularity theory, [53].

Any point source or object far below the resolution limit of the telescope will be mapped to multiple image points by the gravitational lensing effect with the angular position on the sky as the quantitative observable. Examples for such cases are multiple images of SNe, FRBs, or most quasars (QSOs), i.e., objects of short duration or with a time-varying intensity, see Table 1 (left).

We denote images above the resolution limit of the data acquisition as unresolved images if these extended images do not show intrinsic features in the brightness profile above noise level, as shown in Table 1 (centre). The brightness profile of these images can be decomposed into multipoles. We obtain the centre of light of the brightness profiles of the images and the quadrupole around the centre of brightness as observables. Higher-order moments may be extractable, as we are currently investigating. From the quadrupole, we obtain the vectors of the semi-minor and semi-major axes with respect to the centre of light. The ellipticity of an unresolved image is calculated from the axis ratio between the semi-minor and the semi-major axis. The orientation of the image is defined as the direction of the semi-major axis in the coordinate system of the telescope observation. Decomposing the brightness profiles into other sets of basis functions that rely on specific profile models may lead to biases when inferring lens properties, as found e.g., in [54,55] for weakly-distorted galaxy brightness profiles. Therefore, any set of basis functions that is based on non-parametric, statistical properties of the brightness profiles is preferred, see e.g., [56] for a moment-based decomposition of weakly-distorted galaxy brightness profiles, or [57,58] for another suitable basis set.

In the optimal case, the source object has a brightness profile with intrinsic features and the lensing effect is so strong that these features can be identified in all images like in Table 1 (right). As quantitative observables in each image, we obtain the vectors of the angular distances between the brightness features. Our approach requires at least two such non-parallel vectors per multiple image. This implies that any unresolved image can also be treated as a resolved image by means of the semi-major and semi-minor axes of the quadrupole. If the source is a galaxy, the brightness features can be star-forming regions. Thus far, three cases of resolved multiple images of background galaxies caused by a galaxy cluster as lens have been observed: five multiple images in the galaxy cluster CL0024+17 as found by [59], a giant arc and its counter image in the galaxy cluster RCS2 032727-132623 as discovered by [60], and three images of a spiral galaxy in the galaxy cluster MACS J1149.5+2223 as first detailed in [61]. In the radio bands, some multiple images of QSOs show resolved features like core-jet structures, see e.g., [62] or [63].

Table 2 shows the most common multiple-image configurations, i.e., the relative positions and orientations of multiple images with respect to each other. The configurations are shown for unresolved images. However, they are equally observed for point-like and resolved images. As we detail in Section 2.2, these configurations are derived from a leading order approximation to the general nonlinear lens mapping. In addition to these configurations, central images close to the lens centre can occur. They remain unobserved in most cases because they are demagnified by the lens and are usually outshone or occluded by a bright galaxy-scale lens or by the brightest cluster galaxy in the centre of a galaxy-cluster-scale lens.

The first case in Table 2 (left) is a fold configuration of two images that are mirror images of each other. This is often expressed as images A and B having opposite parity (handedness). The second case in Table 2 (centre) is a cusp configuration of three multiple images arranged in a parabolic shape. To leading order, this configuration can be decomposed into the superposition of two fold configurations formed by images A and B and A and C, as further detailed in Section 2.4. The last case in Table 2 (right) occurs for (almost) axisymmetric lenses if the position of the background source falls close to the point in the source plane to which the symmetry centre of the lens is back-projected. For a perfect alignment of the source position and the symmetry centre of the lens along the observer’s line of sight, the source is mapped into an Einstein ring. Small deviations from this alignment or from a perfect axisymmetric lens lead to two giant arcs, a larger one and a smaller one. The latter is called the counter-image. These arcs surround the lens as shown in Table 2 (right). Sometimes, the larger arc can be split further into a cusp configuration, if the deviations from an axisymmetric mass-density distribution are so strong that the ellipticity of the deflecting mass profile becomes relevant.

Measuring the redshift of the images and the lens, the distances between the observer, the source, and the lens can be determined, employing the distance measures of the background metric. For a spatially flat FLRW metric, the angular diameter distances between two redshifts $z_1$ and $z_2$ with $z_2>z_1$ are given as

$$D_{\mathrm{A}}(z_1,z_2)={\displaystyle \frac{c}{H_0}}{\displaystyle \frac{z_1+1}{z_2+1}}\underset{z_1}{\overset{z_2}{\int }}{\displaystyle \frac{\mathrm{d}z}{E\left(z\right)}},$$

(1)

in which $E\left(z\right)$ is the expansion function of the universe. For the standard cold-dark-matter model with a cosmological constant $\Lambda$ ($\Lambda$CDM) based on a spatially flat FLRW metric, it reads

$$E\left(z\right)=\sqrt{{\Omega }_r{(1+z)}^4+{\Omega }_m{(1+z)}^3+{\Omega }_{\Lambda }}.$$

(2)

The ${\Omega }_i$, $i=r,m,\Lambda$, denote the dimensionless contributions of radiation, matter, and the cosmological constant to the total energy density of the universe, normalised such that $E\left(0\right)=1$. The term for the spatial curvature with ${\Omega }_K$ is discarded, as [64] found that the curvature is compatible with zero. Alternatively, any data-based expansion function can be inserted, as, for instance set up for our approach in [26]. Section 2.5 details these data-based distance measures. The distances between the observer, the source, and the lens are only relevant in the time-delay equation, i.e., for cases involving images of time-varying intensity, as detailed in Section 2.4.

2.2. The Standard Gravitational Lensing Formalism

Based on the physical prerequisites of Section 2 and the detailed derivation in [27], we set up the standard gravitational lensing formalism for one lens plane as follows. First, the light propagation along the line of sight is divided into the light propagation through the flat FLRW metric and the light deflection at the lens plane. Relative to the unperturbed light ray that propagates along a null-geodesic of the background metric, let the light be subject to a time delay $t_{\mathrm{s}}$ due to the deflection in the lens plane. Due to this deflection, the light also travels the additional time $t_{\mathrm{g}}$ in the flat FLRW metric relative to the unperturbed ray. Minimising the total travel time with respect to an unperturbed light ray, $t=t_{\mathrm{g}}+t_{\mathrm{s}}$, we obtain

$$t={\displaystyle \frac{(1+z_{\mathrm{l}})}{c}}{\displaystyle \frac{D_{\mathrm{A}}(0,z_{\mathrm{s}})D_{\mathrm{A}}(0,z_{\mathrm{l}})}{D_{\mathrm{A}}(z_{\mathrm{l}},z_{\mathrm{s}})}}\left({\displaystyle \frac{1}{2}}\mathit{\alpha }{\left(\mathit{x}\right)}^2-\psi \left(\mathit{x}\right)\right),$$

(3)

for a light ray starting at the source redshift $z_{\mathrm{s}}$ and being deflected by a two-dimensional gravitational potential $\psi \left(\mathit{x}\right)$ at redshift $z_{\mathrm{l}}$, as shown in Figure 1. $\mathit{\alpha }\left(\mathit{x}\right)\in {\mathbb{R}}^2$ denotes the deflection angle caused by the gravitational potential in the lens plane. It is given by the difference of the angular position of a multiple image in the lens plane, $\mathit{x}=(x_1,x_2)\in {\mathbb{R}}^2$, and the angular position of the source in the source plane, $\mathit{y}=(y_1,y_2)\in {\mathbb{R}}^2$,

$$\mathit{\alpha }\left(\mathit{x}\right)=\mathit{x}-\mathit{y}.$$

(4)

In addition, the deflection angle $\mathit{\alpha }\left(\mathit{x}\right)$ is given as the gradient of the projected deflection potential $\psi \left(\mathit{x}\right)$ with respect to $\mathit{x}$, $\mathit{\alpha }\left(\mathit{x}\right)=\nabla \psi \left(\mathit{x}\right)$ such that Equation (4) is often stated as

$$\mathit{y}=\mathit{x}-\nabla \psi \left(\mathit{x}\right).$$

(5)

For the sake of convenience, we abbreviate

$$\Gamma \equiv {\displaystyle \frac{(1+z_{\mathrm{l}})}{c}}{\displaystyle \frac{D_{\mathrm{A}}(0,z_{\mathrm{s}})D_{\mathrm{A}}(0,z_{\mathrm{l}})}{D_{\mathrm{A}}(z_{\mathrm{l}},z_{\mathrm{s}})}},\varphi (\mathit{x},\psi )\equiv {\displaystyle \frac{1}{2}}{\left(\mathit{x}-\mathit{y}\right)}^2-\psi \left(\mathit{x}\right)={\displaystyle \frac{1}{2}}{\left(\mathit{\alpha }\left(\mathit{x}\right)\right)}^2-\psi \left(\mathit{x}\right).$$

(6)

$\varphi (\mathit{x},\psi )$ is called the Fermat potential. From the viewpoint of singularity theory, $\varphi (\mathit{x},\psi )$ is a function of $\mathit{x}$ parametrised by $\psi \left(\mathit{x}\right)$. Its critical points, $\nabla \varphi (\mathit{x},\psi )=0$, are the positions of the multiple images given by Equation (5). The second derivative, the Hessian matrix of $\varphi (\mathit{x},\psi )$, also called distortion matrix,

$$A\left(\mathit{x}\right)=\left(\begin{array}{cc}{\varphi }_{11}\left(\mathit{x}\right)& {\varphi }_{12}\left(\mathit{x}\right)\\ {\varphi }_{12}\left(\mathit{x}\right)& {\varphi }_{22}\left(\mathit{x}\right)\end{array}\right),{\varphi }_{\alpha \beta }\left(\mathit{x}\right)\equiv {\displaystyle \frac{{\partial }^2\varphi (\mathit{x},\psi )}{\partial x_{\alpha }\partial x_{\beta }}},\alpha ,\beta =1,2,$$

(7)

is a two-dimensional mapping from vectors in the lens plane around $\mathit{x}$ to vectors in the source plane around $\mathit{y}$. Hence, it describes the distortion of light bundles in the vicinity of $\mathit{x}$ caused by the gravitational lensing effect. $A\left(\mathit{x}\right)$ is often parametrised by the scaled, projected mass density called convergence $\kappa \left(\mathit{x}\right)$ and the shear $\gamma \left(\mathit{x}\right)=({\gamma }_1\left(\mathit{x}\right),{\gamma }_2\left(\mathit{x}\right))$, such that $A\left(\mathit{x}\right)$ is decomposed into a magnifying trace part and a distorting symmetric part

$$A\left(\mathit{x}\right)=\left(\begin{array}{cc}1-\kappa \left(\mathit{x}\right)& 0\\ 0& 1-\kappa \left(\mathit{x}\right)\end{array}\right)-\left(\begin{array}{cc}{\gamma }_1\left(\mathit{x}\right)& {\gamma }_2\left(\mathit{x}\right)\\ {\gamma }_2\left(\mathit{x}\right)& -{\gamma }_1\left(\mathit{x}\right)\end{array}\right)=(1-\kappa \left(\mathit{x}\right))\left(\begin{array}{cc}1-g_1\left(\mathit{x}\right)& -g_2\left(\mathit{x}\right)\\ -g_2\left(\mathit{x}\right)& 1+g_1\left(\mathit{x}\right)\end{array}\right).$$

(8)

$\kappa \left(\mathit{x}\right)$ is related to the deflection potential by the Poisson equation $\Delta \psi \left(\mathit{x}\right)=2\kappa \left(\mathit{x}\right)$. $\gamma \left(\mathit{x}\right)$ can be physically interpreted as the deflection caused by parts of the mass density surrounding $\mathit{x}$. The shear thus represents the non-local part of the deflection potential that is not included in the Poisson equation. It enters $\psi \left(\mathit{x}\right)$ via the boundary conditions, as detailed in [27]. $\mathit{g}\left(\mathit{x}\right)=(g_1\left(\mathit{x}\right),g_2\left(\mathit{x}\right))$ is the reduced shear, which is the quantity constrained by observables, as detailed in Section 2.3. Given $A\left(\mathit{x}\right)$ either as parametrised by Equations (7) or (8), it is possible to reconstruct vectors in the source plane in the vicinity of $\mathit{x}$ up to an overall scale factor, as also detailed in Section 2.3.

The set of degenerate critical points ${\mathit{x}}_0$ for which $\nabla \varphi ({\mathit{x}}_0,\psi )=0$ and $\det \left(A\left({\mathit{x}}_0\right)\right)=0$ are called the critical curves. Using Equations (7) and (8), they are given by all ${\mathit{x}}_0$ fulfilling

$$\det \left(A\left({\mathit{x}}_0\right)\right)={\varphi }_{11}\left({\mathit{x}}_0\right){\varphi }_{22}\left({\mathit{x}}_0\right)-{\varphi }_{12}{\left({\mathit{x}}_0\right)}^2={\left(1-\kappa \left({\mathit{x}}_0\right)\right)}^2-\gamma {\left({\mathit{x}}_0\right)}^2=0.$$

(9)

As shown in [53], the critical curves described by Equation (9) consist of only two kinds of critical points, so-called fold points, ${\mathit{x}}_{\mathrm{f}}$, and cusp points, ${\mathit{x}}_{\mathrm{c}}$. Irrespective of the specific form of $\psi \left(\mathit{x}\right)$, ${\mathit{x}}_{\mathrm{f}}$ and ${\mathit{x}}_{\mathrm{c}}$ are determined by a Taylor expansion of the Fermat potential of order four around each critical point.

Thus, we have two possibilities to determine local lens properties without assuming a specific form of $\psi \left(\mathit{x}\right)$: first, we assume that $A\left(\mathit{x}\right)$ is constant in the vicinity of the centre of light of a multiple image. This approximation is valid if the extensions of a multiple image are small compared to the scale on which the deflecting mass density changes. A comprehensive synopsis of the validity of these approximations for extended sources and finite beam sizes in gravitational lensing can be found in [65]. In [66] the most general equations to constrain local lens properties from Equation (5) were derived. In [23], we developed this idea further by assuming a constant $A\left(\mathit{x}\right)$ in the area spanned by the quadrupole of an unresolved multiple image or in the area spanned by the convex hull of all identifiable reference points of a resolved multiple image. Section 2.3 details the ansatz and shows the local lens properties that are determined in the vicinity of the multiple images. They are calculated using the quadrupoles of at least two unresolved images in a fold configuration. Alternatively, the ansatz leads to the same local lens properties for at least three identifiable reference points in at least two images in a fold configuration.

Second, we perform the Taylor expansion to fourth order around a critical point and employ the observables of the multiple images to infer local properties of the critical curves, as detailed in Section 2.4. Consequently, local lens properties in the vicinity of the multiple images and approximations to the critical curves can be determined independently of the specific global form of $\psi \left(\mathit{x}\right)$, i.e., without any prior assumptions about the decomposition of the deflecting gravitational potential.

2.3. Local Lens Properties from a Taylor Expansion around the Centre of Light of the Multiple Images

The work in Ref. [62] already described the possibility to extract the local lens properties contained in $A\left(\mathit{x}\right)$ by transforming multiple images onto each other. In [23,47], we realised this idea for unresolved and resolved multiple images that contain at least two non-parallel vectors, as indicated by the black arrows in the figures of Table 1. We assume to have m reference points within each multiple image i which are located at positions ${\mathit{x}}_{i\alpha }$, $\alpha =1,\dots ,m$. Using the quadrupole of an unresolved image as observables, the three reference points are the centre of light and the end points of the semi-major and semi-minor axis of the quadrupole. Furthermore, we assume that the multiple images have a very small extent compared to the scale on which the mass density changes. Section 3.2 and Section 3.3 show the validity of this assumption in two practical applications on cluster and on galaxy scale. This implies that the entries in $A\left(\mathit{x}\right)$ are approximately constant over the area spanned by the ${\mathit{x}}_{i\alpha }$ in each multiple image. Since the vectors in all multiple images ${\mathit{x}}_{i\alpha }-{\mathit{x}}_{i\beta }$, $\alpha ,\beta =1,\dots ,m$, originate from the same, corresponding vectors ${\mathit{y}}_{\alpha }-{\mathit{y}}_{\beta }$ in the source plane, we note that for each of these vectors

$${\mathit{y}}_{\alpha }-{\mathit{y}}_{\beta }=A\left({\mathit{x}}_j\right)\left({\mathit{x}}_{j\alpha }-{\mathit{x}}_{j\beta }\right)=A\left({\mathit{x}}_i\right)\left({\mathit{x}}_{i\alpha }-{\mathit{x}}_{i\beta }\right),\forall i,j=A,B,C,\dots \forall \alpha ,\beta =1,2,\dots ,m$$

(10)

holds. As the distortion matrix is assumed to be constant over the area spanned by the reference points, ${\mathit{x}}_i$ and ${\mathit{x}}_j$ represent an arbitrary point within this area. Reformulating the part on the right-hand side of Equation (10) yields a linear transformation $T_{ij}$ between vectors in image j to vectors in image i as

$${\mathit{x}}_{j\alpha }-{\mathit{x}}_{j\beta }=A^{-1}\left({\mathit{x}}_j\right)A\left({\mathit{x}}_i\right)\left({\mathit{x}}_{i\alpha }-{\mathit{x}}_{i\beta }\right)\equiv T_{ij}\left({\mathit{x}}_{i\alpha }-{\mathit{x}}_{i\beta }\right),\forall i,j=A,B,C,\dots \forall \alpha ,\beta =1,2,\dots ,m.$$

(11)

Using Equation (11) to solve for the entries in $A\left(\mathit{x}\right)$ given the reference points ${\mathit{x}}_{i\alpha }$ is a total-least-squares parameter estimation problem. It differs from the standard least-squares parameter estimation because the vectors between the reference points on both sides of the equation are subject to measurement uncertainties. To account for this, we introduce one latent variable ${\mathit{x}}_i$ per multiple image, the so-called anchor points in [23,47].

Solving the total-least-squares parameter estimation yields ratios of convergences, $f_{ij}$, between pairs of multiple images i and j and the reduced shear $g\left(\mathit{x}\right)$ at the positions of the multiple images

$$f_{ij}^{\left(\kappa \right)}\equiv {\displaystyle \frac{1-\kappa \left({\mathit{x}}_i\right)}{1-\kappa \left({\mathit{x}}_j\right)}},g_i\equiv g\left({\mathit{x}}_i\right)={\displaystyle \frac{\gamma \left({\mathit{x}}_i\right)}{1-\kappa \left({\mathit{x}}_i\right)}}\forall i,j=A,B,C,\dots .$$

(12)

It also yields the positions of the anchor points ${\mathit{x}}_i$. Alternatively, expressed in derivatives of the Fermat potential, we can infer the analogous quantities

$$f_{ij}^{\left(\varphi \right)}={\displaystyle \frac{{\varphi }_{11}\left({\mathit{x}}_i\right)}{{\varphi }_{11}\left({\mathit{x}}_j\right)}},{\tilde{\varphi }}_{12}\left({\mathit{x}}_i\right)={\displaystyle \frac{{\varphi }_{12}\left({\mathit{x}}_i\right)}{{\varphi }_{11}\left({\mathit{x}}_i\right)}},{\tilde{\varphi }}_{22}\left({\mathit{x}}_i\right)={\displaystyle \frac{{\varphi }_{22}\left({\mathit{x}}_i\right)}{{\varphi }_{11}\left({\mathit{x}}_i\right)}}\forall i,j=A,B,C,\dots .$$

(13)

In [23], we performed a comparison between the different parametrisations of $A\left(\mathit{x}\right)$ by Equations (12) and (13) on simulated data. We found that both parametrisations have their advantages and drawbacks, so that it depends on the image configuration, which one yields more accurate, precise, and robust results. Furthermore, we showed that the larger the area of each multiple image from which the observables are extracted, the smaller the confidence bounds on the inferred local lens properties. On the other hand, the area should be so small that $A\left(\mathit{x}\right)$ is still approximately constant for all $\mathit{x}$ within the area. In any case, the approach needs at least three non-collinear reference points ${\mathit{x}}_{i\alpha }$, $\alpha \ge 3$, in each multiple image i and at least three images to return a unique solution for the local lens properties as given by Equations (12) or (13). The system of equations is under-determined, if $T_{ij}$ is diagonal, so that we cannot determine local lens properties from aligned multiple-image configurations created by an axisymmetric lens (see [23] for a detailed explanation). In addition, $\kappa \left({\mathit{x}}_j\right)=1$ should be avoided in the denominator of $f_{ij}$ in Equation (12), such that $f_{ij}$ does not become singular. Apart from these restrictions, it handles any configuration of at least three images including central maxima or counter images of cusp configurations. Figure 2 visualises the principle to obtain Equations (12) and (13). An implementation to solve for Equation (12) is publicly available at https://github.com/ntessore/imagemap and details about it can be found in [47].

As mentioned in Section 2.2, Equation (10) shows that vectors in the source plane, ${\mathit{y}}_{\alpha }-{\mathit{y}}_{\beta }$ can be determined from the back-projection of the respective vectors in the image plane, ${\mathit{x}}_{i\alpha }-{\mathit{x}}_{i\beta }$, using $A\left({\mathit{x}}_i\right)$ as constrained by the local lens properties in Equations (12) or (13). Comparing Equations (12) and (13) with Equations (7) and (8), we find that ${\mathit{y}}_{\alpha }-{\mathit{y}}_{\beta }$ are only determined up to an overall scale factor because the observables do not constrain $\kappa \left(\mathit{x}\right)$ but only ratios of convergences, $f_{ij}$. An example source reconstruction using this method is performed in Section 3.2 and clearly shows that the reconstructed morphology of the source can give useful information about the relative shapes and sizes of the star-forming regions of a multiply-imaged galaxy. A more detailed analysis and comparison to lens-model-based source reconstructions is currently in preparation.

From the local lens properties given by Equations (12) and (13), we calculate the relative image parity information between the multiple images i and j,

$${\mathcal{J}}_{ij}\equiv \det T_{ij}=\det {\left(A\left({\mathit{x}}_j\right)\right)}^{-1}\det \left(A\left({\mathit{x}}_i\right)\right)\forall i,j=A,B,C,\dots .$$

(14)

For resolved multiple images, the relative parity can be read off the telescope observation by visual inspection as a consistency check. In addition, ${\mathcal{J}}_{ij}$ is also the magnification ratio between the two images i and j. Hence, it can be compared to the measured flux density ratios, if available, to indicate extinction due to dust, micro-lensing, or higher order lensing effects beyond convergence and shear.

The Fold Case

In order to uniquely solve Equation (11) for a fold configuration, we have to reduce the amount of variables in Equations (12) and (13). Therefore, we approximate $f_{ij}\approx 1$ in Equations (12) and (13), assuming that the convergence or the ${\varphi }_{11}\left(\mathit{x}\right)$ are equal for the areas covered by the two fold images. Thus, only the reduced shear can be determined for the two images at a fold.

The validity of the approximation $f_{ij}\approx 1$ was investigated in [23] for two analytic mass density models (the ones discussed in Appendix A, for choices of parameters, see [23]). We found that the decrease in accuracy of the approximation $f_{ij}\approx 1$ depends on the underlying mass density. Even for our simulated extreme case of a highly elliptical lens, at least one of the parameterisations in Equations (12) and (13) yielded local lens properties that were within the 2$-\sigma$ confidence intervals around the true values.

2.4. Local Lens Properties from a Taylor Expansion around a Critical Point

Next, we Taylor-expand the Fermat potential around a critical point ${\mathit{x}}_0$, whose source is located at the so-called caustic point ${\mathit{y}}_0$. Without loss of generality, we choose the coordinate systems in the lens and the source plane such that the critical and the caustic point are at the origin, ${\mathit{x}}_0={\mathit{y}}_0=\mathbf{0}$. Inserting Equation (5) into the Fermat potential, the Taylor series to order four reads

$$\begin{array}{ll}{\varphi }_{\mathrm{T}}(\mathit{x},\psi )\equiv \hfill & \varphi ({\mathit{x}}_0,\psi )+\frac{1}{2}{\mathit{y}}^2-\mathit{x}\mathit{y}+\frac{1}{2}{\varphi }_{11}\left({\mathit{x}}_0\right)x_1^2+{\varphi }_{12}\left({\mathit{x}}_0\right)x_1{x}_2+\frac{1}{2}{\varphi }_{22}\left({\mathit{x}}_0\right)x_2^2\hfill \\ & +\frac{1}{6}{\varphi }_{111}\left({\mathit{x}}_0\right)x_1^3+\frac{1}{2}{\varphi }_{112}\left({\mathit{x}}_0\right)x_1^2{x}_2+\frac{1}{2}{\varphi }_{122}\left({\mathit{x}}_0\right)x_1{x}_2^2+\frac{1}{6}{\varphi }_{222}\left({\mathit{x}}_0\right)x_2^3\hfill \\ & +\frac{1}{24}{\varphi }_{1111}\left({\mathit{x}}_0\right)x_1^4+\frac{1}{6}{\varphi }_{1112}\left({\mathit{x}}_0\right)x_1^3{x}_2+\frac{1}{4}{\varphi }_{1122}\left({\mathit{x}}_0\right)x_1^2{x}_2^2+\frac{1}{6}{\varphi }_{1222}\left({\mathit{x}}_0\right)x_1{x}_2^3+\frac{1}{24}{\varphi }_{2222}\left({\mathit{x}}_0\right)x_2^4.\hfill \end{array}$$

(15)

As ${\mathit{x}}_0$ is a critical point, Equation (9) must be fulfilled for the second-order derivatives in Equation (15). Using Equation (15) as the local Fermat potential, the time-delay difference between light rays arriving from two multiple images i and j that are located close to ${\mathit{x}}_0$ is given by

$$\begin{array}{c}{\tau }_{ij}=t_i-t_j=\Gamma {\varphi }_{\mathrm{T}}({\mathit{x}}_i,\psi )-\Gamma {\varphi }_{\mathrm{T}}\left({\mathit{x}}_j\psi \right).\hfill \end{array}$$

(16)

$$\nabla {\varphi }_{\mathrm{T}}(\mathit{x},\psi )=0$$

(17)

yields the lensing equation to determine the positions of multiple images in the vicinity of ${\mathit{x}}_0$. From Equations (16) and (17), a system of equations is set up that contains $\mathit{y}$ and derivatives of the Fermat potential at ${\mathit{x}}_0$ as unknowns. The observables, introduced in Section 2.1, are employed to determine the angular positions $\mathit{x}$ of the multiple images with respect to ${\mathit{x}}_0$, as detailed in the following sections.

Depending on the type of critical point, ${\mathit{x}}_{\mathrm{f}}$ or ${\mathit{x}}_{\mathrm{c}}$, different further approximations to Equation (15) are necessary to reduce the number of unknowns, such that the system can be uniquely solved.

Table 3. Synopsis of approximations to Equation (15) and embedding of the observables for fold and cusp critical points, ${\mathit{x}}_{\mathrm{f}}$ and ${\mathit{x}}_{\mathrm{c}}$, in the coordinate system of Equation (18). The critical point is marked by a black dot in the figures. Degeneracies arise in the sign of the slope of the critical curve at ${\mathit{x}}_{\mathrm{f}}$. Depending on the parity of image A, there are two possible configurations for a cusp.

Fold	Cusp
	(Positive)	(Negative)
coordinate system specifications ${\varphi }_{222}\left({\mathit{x}}_0\right)>0$,	coordinate system specifications ${\varphi }_{122}\left({\mathit{x}}_0\right)<0$, ${\varphi }_{2222}\left({\mathit{x}}_0\right)>0$, ${\varphi }_{122}{\left({\mathit{x}}_0\right)}^2-\frac{1}{3}{\varphi }_{11}\left({\mathit{x}}_0\right){\varphi }_{2222}\left({\mathit{x}}_0\right)\ne 0$
leading order ${\varphi }_{\mathrm{T}}(\mathit{x},\psi )=$ $-\mathit{x}\mathit{y}+\frac{1}{2}{\varphi }_{11}\left({\mathit{x}}_0\right)x_1^2+\frac{1}{2}{\varphi }_{112}\left({\mathit{x}}_0\right)x_1^2{x}_2$ $+\frac{1}{2}{\varphi }_{122}\left({\mathit{x}}_0\right)x_1{x}_2^2+\frac{1}{6}{\varphi }_{222}\left({\mathit{x}}_0\right)x_2^3$	leading order ${\varphi }_{\mathrm{T}}(\mathit{x},\psi )=$ $-\mathit{x}\mathit{y}+\frac{1}{2}{\varphi }_{11}\left({\mathit{x}}_0\right)x_1^2$ $+\frac{1}{2}{\varphi }_{122}\left({\mathit{x}}_0\right)x_1{x}_2^2+\frac{1}{24}{\varphi }_{2222}\left({\mathit{x}}_0\right)x_2^4$

shows these approximations to Equation (17) to obtain the lensing equations for Section 2.4.1 and Section 2.4.2. The approximations are inspired by axisymmetric and elliptical lens models (see Appendix A for two examples) and can be motivated by the fact that Newtonian gravity (as also employed for the description of the gravitational lens), is a central force. The goodness of the approximating assumptions will increase with decreasing distance of the multiple images to the centre of the potential of the deflecting galaxy or galaxy cluster and with increasing depth of this potential with respect to neighbouring gravitationally interacting objects.

Diagonalising $A\left({\mathit{x}}_0\right)$ introduces a special coordinate system, which we choose, as in [50],

$${\varphi }_1\left({\mathit{x}}_0\right)={\varphi }_2\left({\mathit{x}}_0\right)={\varphi }_{12}\left({\mathit{x}}_0\right)={\varphi }_{22}\left({\mathit{x}}_0\right)=0,{\varphi }_{11}\left({\mathit{x}}_0\right)\ne 0$$

(18)

without loss of generality. The transformation into this coordinate frame simplifies the equations, but also complicates the combination of several multiple-image configurations in a mutual telescope observation because each multiple-image system has its own frame of reference.

Table 3 summarises the specifications that we employ to embed the observables of Section 2.1 into the coordinate system given by Equation (18).

Similarly to Section 2.3, the observables allow for determining the distortion matrices at the image positions, $A\left({\mathit{x}}_i\right)$, $i=A,B,C$, as described in Section 2.4.1 and Section 2.4.2. Here, ${\mathit{x}}_i$ is the centre of light of a multiple image. Again, these distortion matrices can be employed to reconstruct vectors in the source plane from corresponding vectors in the image plane close to ${\mathit{x}}_i$ up to an overall scale factor.

2.4.1. The Fold Case

The work in Ref. [50] derived that, to leading order, the fold configuration consists of two images that are symmetrically located around $x_{\mathrm{f}}$ as shown in Table 3 (left figure). Employing ${\varphi }_{\mathrm{T}}(\mathit{x},\psi )$ in the coordinate system as introduced in Table 3, we derive the time delay difference between light rays arriving from the two images by Equation (16) and the lensing equation by Equation (17). The distortion matrix as the second derivative of ${\varphi }_{\mathrm{T}}(\mathit{x},\psi )$ with respect to $\mathit{x}$ reads

$$A\left({\mathit{x}}_i\right)={\varphi }_{11}\left({\mathit{x}}_{\mathrm{f}}\right)\left(\begin{array}{cc}1& \frac{{\varphi }_{122}\left({\mathit{x}}_{\mathrm{f}}\right)}{{\varphi }_{11}\left({\mathit{x}}_{\mathrm{f}}\right)}{x}_{i2}\\ \frac{{\varphi }_{122}\left({\mathit{x}}_{\mathrm{f}}\right)}{{\varphi }_{11}\left({\mathit{x}}_{\mathrm{f}}\right)}{x}_{i2}& \frac{{\varphi }_{222}\left({\mathit{x}}_{\mathrm{f}}\right)}{{\varphi }_{11}\left({\mathit{x}}_{\mathrm{f}}\right)}{x}_{i2}\end{array}\right),i=A,B.$$

(19)

Hence, we make the approximation that ${\varphi }_{11}\left({\mathit{x}}_i\right)\approx {\varphi }_{11}\left({\mathit{x}}_{\mathrm{f}}\right)$, $i=A,B$ and interpolate $A\left(\mathit{x}\right)$ linearly between ${\mathit{x}}_{\mathrm{f}}$ and ${\mathit{x}}_i$. We factor $A\left({\mathit{x}}_i\right)$ into an overall scale, ${\varphi }_{11}\left({\mathit{x}}_{\mathrm{f}}\right)$ and a reduced distortion matrix, analogous to Equation (8). Comparing Equation (19) to Equation (13), we arrive again at the approximation made in Section 2.3.1. In addition, the mirror symmetry implies that the images are of different parity, which can easily be verified by inserting $x_{i2}$ into Equation (19).

Due to the symmetry, $x_{A2}={\delta }_{AB}/2$ and $x_{B2}=-{\delta }_{AB}/2$, in which the distance ${\delta }_{AB}$ between ${\mathit{x}}_A$ and ${\mathit{x}}_B$ is an observable. Further observables are the quadrupole of each image around the centre of light, ${\mathit{x}}_i$, parametrised by the axis ratio $r_i$ of the semi-minor to the semi-major axis and the orientation angle of the semi-major axis to ${\delta }_{AB}$, ${\phi }_i$, $i=A,B$. The validity of the symmetry assumption can thus be checked by comparing the observables of the two quadrupoles with each other. A measured flux ratio close to one between images A and B additionally corroborates the symmetry assumption, if dust extinction and micro-lensing are negligible effects (see also Section 2.3). If the symmetry assumption is applicable, the images are close enough to the critical curve to assume that the source properties are negligible. In the appendix of [22], we generally showed that the influence of the intrinsic quadrupole of the source decreases for decreasing distance to the critical point. Consequently, we assume that the observed quadrupole of the multiple images is caused by the gravitational lensing effect. We set it equal to the reduced distortion matrix in Equation (19).

Altogether, we can determine the following leading-order local lens properties from a fold configuration of resolved or unresolved multiple images with a measured time delay difference of ${\tau }_{AB}$

$$\begin{array}{}\mathrm{(20)}& \hfill {\varphi }_{222}\left({\mathit{x}}_{\mathrm{f}}\right)& ={\displaystyle \frac{12{\tau }_{AB}}{\Gamma {\delta }_{AB}^3}},\hfill \mathrm{(21)}& \hfill \left|{\tilde{\varphi }}_{122}\left({\mathit{x}}_{\mathrm{f}}\right)\right|& \equiv \left|{\displaystyle \frac{{\varphi }_{122}\left({\mathit{x}}_{\mathrm{f}}\right)}{{\varphi }_{11}\left({\mathit{x}}_{\mathrm{f}}\right)}}\right|={\displaystyle \frac{2\left|{\tilde{\varphi }}_{12}\left({\mathit{x}}_i\right)\right|}{{\delta }_{AB}}}={\displaystyle \frac{2}{{\delta }_{AB}}}·\left|{\displaystyle \frac{\tan \left({\phi }_i\right)(r_i-1)}{r_i{\tan }^2\left({\phi }_i\right)+1}}\right|,i=A,B,{\phi }_i\in \left[-\frac{\pi }{2},\frac{\pi }{2}\right],\hfill \mathrm{(22)}& \hfill {\tilde{\varphi }}_{222}\left({\mathit{x}}_{\mathrm{f}}\right)& \equiv {\displaystyle \frac{{\varphi }_{222}\left({\mathit{x}}_{\mathrm{f}}\right)}{{\varphi }_{11}\left({\mathit{x}}_{\mathrm{f}}\right)}}={\displaystyle \frac{2{\tilde{\varphi }}_{22}\left({\mathit{x}}_i\right)}{{\delta }_{AB}}}={\displaystyle \frac{2}{{\delta }_{AB}}}·{\displaystyle \frac{r_i+{\tan }^2\left({\phi }_i\right)}{r_i{\tan }^2\left({\phi }_i\right)+1}}\approx {\displaystyle \frac{2r_i}{{\delta }_{AB}}},\hfill \mathrm{(23)}& \hfill m_{\mathrm{f}}& \equiv {\displaystyle \frac{\left|{\varphi }_{122}\left({\mathit{x}}_{\mathrm{f}}\right)\right|}{{\varphi }_{222}\left({\mathit{x}}_{\mathrm{f}}\right)}},\hfill \end{array}$$

in which the slope of the critical curve at ${\mathit{x}}_{\mathrm{f}}$, $m_{\mathrm{f}}$, is determined by the ratio of Equations (21) and (22). The sign of the slope cannot be determined. This requires a third image, as discussed in Section 2.4.2. In [22], we found that the approximation assuming ${\phi }_i=0$ on the right-hand side of Equation (22) yields the most accurate results for this ratio of derivatives.

We note that Equations (21)–(23) do not contain any information about the angular diameter distances to the lens or the source plane. As derived in [25] and explained in [27], the former fact arises because the Fermat potential is an angular, projected potential, and the multiple-image shapes are the corresponding angular observables on the celestial sphere. The entire geometric information is contained in $\Gamma$ (see Equation (6)) and therefore only relevant in Equation (20). Hence, the gravitational lensing geometry, and properties of the underlying cosmology, can only be constrained with measured time delay differences by this model-independent approach.

Furthermore, Equations (21)–(23) only constrain ratios of derivatives. This is a consequence of the fact that the deflecting mass density is an inhomogeneity defined on top of a background mass density. Rescaling the background mass density and scaling the deflecting mass density accordingly is a choice of parametrisation and does not change the observables. As detailed in [27], Equation (20) breaks this degeneracy for fixed $\Gamma$, i.e., fixing the background cosmology and thereby the background mass density at ${\mathit{x}}_A$ and ${\mathit{x}}_B$.

Linearising Equation (5) and applying it to fields of single weakly distorted galaxies, [67] already found that weak-gravitational-lensing observations only constrain the reduced shear $\mathit{g}\left(\mathit{x}\right)$ and not $\gamma \left(\mathit{x}\right)$ because the weak-lensing equations are also subject to the mass-sheet degeneracy (MSD) that is intrinsic in the lensing formalism, [68]. This finding and the MSD can now be understood with the results as stated above. Further details can be found in Section 3.6 and in [27].

2.4.2. The Cusp Case

Analogously to the derivations performed for Section 2.4.1, the leading-order local lens properties for the cusp configuration as embedded into the coordinate system shown in Table 3 (centre and right figure) can be determined. The distortion matrix which is set equal to the quadrupoles of the multiple images is given by

$$A\left({\mathit{x}}_i\right)={\varphi }_{11}\left({\mathit{x}}_{\mathrm{c}}\right)\left(\begin{array}{cc}1& \frac{{\varphi }_{122}\left({\mathit{x}}_{\mathrm{c}}\right)}{{\varphi }_{11}\left({\mathit{x}}_{\mathrm{c}}\right)}{x}_{i2}\\ \frac{{\varphi }_{122}\left({\mathit{x}}_{\mathrm{c}}\right)}{{\varphi }_{11}\left({\mathit{x}}_{\mathrm{c}}\right)}{x}_{i2}& \frac{{\varphi }_{122}\left({\mathit{x}}_{\mathrm{c}}\right)}{{\varphi }_{11}\left({\mathit{x}}_{\mathrm{c}}\right)}{x}_{i1}+\frac{1}{2}\frac{{\varphi }_{2222}\left({\mathit{x}}_{\mathrm{c}}\right)}{{\varphi }_{11}\left({\mathit{x}}_{\mathrm{c}}\right)}{x}_{i2}^2\end{array}\right),i=A,B,C.$$

(24)

The $x_2$-axis of the special coordinate system of the fold is aligned with the relative distance between the two images, such that the observables directly correspond to the variables in the equations. Contrary to that, the coordinate system for the cusp configuration is constructed such that its $x_1$-axis is the symmetry axis of the parabolic approximation of the critical curve around ${\mathit{x}}_{\mathrm{c}}$. Therefore, the quadrupoles of the multiple images have to be rotated into this coordinate system first, i.e., calculating the ${\phi }_i$ with respect to the positive $x_1$-axis from the relative angles ${\tilde{\phi }}_i$ as shown in Table 1. Appendix B shows the equation to obtain the rotation angle which transforms a configuration as shown in Table 1 into the coordinate system shown in Table 3. In this appendix, we also give the coordinates of the image closest to the cusp, ${\mathit{x}}_A$. This determines the position of ${\mathit{x}}_{\mathrm{c}}$ with respect to ${\mathit{x}}_A$ from the observables, analogous to ${\mathit{x}}_{\mathrm{f}}$ being half the distance between image A and B in Section 2.4.1.

Without loss of generality, we denote the two images that are closest to ${\mathit{x}}_{\mathrm{c}}$ by A and B, as shown in Table 3. Any combination of two images from $A,B,C$ can be used to determine the following leading-order local lens properties. The most accurate results are obtained for images A and B closest to the expansion point of the Taylor approximation, as we showed in [22]. Summarising the calculations of [21,22,23], we obtain as leading-order lens properties from the relative distances between the multiple images, the quadrupoles, and a time delay difference of ${\tau }_{AB}$ between images A and B

$$\begin{array}{}\mathrm{(25)}& \hfill {\varphi }_{2222}\left({\mathit{x}}_{\mathrm{c}}\right)& ={\displaystyle \frac{8{\tau }_{AB}}{\Gamma {\delta }_{AB2}^4}},\hfill \mathrm{(26)}& \hfill {\tilde{\varphi }}_{122}\left({\mathit{x}}_{\mathrm{c}}\right)& \equiv {\displaystyle \frac{{\varphi }_{122}\left({\mathit{x}}_{\mathrm{c}}\right)}{{\varphi }_{11}\left({\mathit{x}}_{\mathrm{c}}\right)}}={\displaystyle \frac{{\tilde{\varphi }}_{12}\left({\mathit{x}}_A\right)-{\tilde{\varphi }}_{12}\left({\mathit{x}}_B\right)}{{\delta }_{AB2}}}={\displaystyle \frac{1}{{\delta }_{AB2}}}\left({\displaystyle \frac{\cot \left({\phi }_A\right)(r_A-1)}{r_A{\cot }^2\left({\phi }_A\right)+1}}-{\displaystyle \frac{\cot \left({\phi }_B\right)(r_B-1)}{r_B{\cot }^2\left({\phi }_B\right)+1}}\right),\hfill \mathrm{(27)}& \hfill {\tilde{\varphi }}_{2222}\left({\mathit{x}}_{\mathrm{c}}\right)& \equiv {\displaystyle \frac{{\varphi }_{2222}\left({\mathit{x}}_{\mathrm{c}}\right)}{{\varphi }_{11}\left({\mathit{x}}_{\mathrm{c}}\right)}}={\displaystyle \frac{-6}{{\delta }_{AB2}^2}}\left({\displaystyle \frac{{\delta }_{AB1}}{{\delta }_{AB2}}}{\tilde{\varphi }}_{12}\left({\mathit{x}}_B\right)+{\tilde{\varphi }}_{22}\left({\mathit{x}}_A\right)\right)\hfill & & ={\displaystyle \frac{-6}{{\delta }_{AB2}^2}}\left({\displaystyle \frac{{\delta }_{AB1}}{{\delta }_{AB2}}}{\displaystyle \frac{\cot \left({\phi }_B\right)(r_B-1)}{r_B{\cot }^2\left({\phi }_B\right)+1}}\stackrel{-}{+}{\displaystyle \frac{r_A+{\cot }^2\left({\phi }_A\right)}{r_A{\cot }^2\left({\phi }_A\right)+1}}\right),{\phi }_A,{\phi }_B\in \left[0,\pi \right],\hfill \mathrm{(28)}& \hfill m_{\mathrm{c}}& \equiv -{\displaystyle \frac{\frac{1}{2}{\tilde{\varphi }}_{2222}\left({\mathit{x}}_{\mathrm{c}}\right)-{\tilde{\varphi }}_{122}{\left({\mathit{x}}_{\mathrm{c}}\right)}^2}{{\tilde{\varphi }}_{122}\left({\mathit{x}}_{\mathrm{c}}\right)}}.\hfill \end{array}$$

As in Section 2.4.1 and for the same reasons, we only obtain ratios of derivatives from the quadrupole observables and, as before, the time delay difference breaks this degeneracy. To leading order, the cusp configuration is a linear superposition of two fold configurations, as investigated in detail in [22,50]. This also becomes apparent when we compare Equations (21) and (26)1. We see that Equation (26) is the sum of the $\left|{\tilde{\varphi }}_{12}\left(\mathit{x}\right)\right|$ of images A and B, while, for Equation (21), images A and B are described as mirror images of each other, so that $\left|{\tilde{\varphi }}_{12}\left(\mathit{x}\right)\right|$ of one image, A or B, is sufficient to determine ${\tilde{\varphi }}_{122}\left({\mathit{x}}_0\right)$. Treating images A and B as a fold configuration, the sign of the slope in Equation (21) can be determined from the position of image C, as the sign of the slope at ${\mathit{x}}_{\mathrm{f}}$ between images A and B should be consistent with the slope of the parabola given by $m_{\mathrm{c}}$. The leading-order Taylor expansion of the cusp configuration requires a rotation into a coordinate system that is not directly aligned with an observable and it needs a fourth order derivative, i.e., an expansion to one order more than the fold configuration. Consequently, the local lens properties determined at ${\mathit{x}}_{\mathrm{c}}$ are less accurate than the ones at ${\mathit{x}}_{\mathrm{f}}$, see [22,23] for details.

The remaining degeneracy is shown in Equation (27). Consistently with Section 2.3 and Section 2.4.1, only the relative parities between the three images are known. Hence, there are two possible pathways for the critical curve around the three images, as further detailed in [53]. The first with the minus-sign in Equation (27) encloses image A in the parabolic approximation to the critical curve (called a positive cusp, as shown in Table 3 (centre)). The second with the plus-sign in Equation (27) encloses images B and C, such that image A lies on the negative side of the $x_1$-axis (called a negative cusp, as shown in Table 3 (right)).

2.4.3. The Giant-Arc Case

“Golden lenses” are highly (axi-)symmetric galaxies that show multiple images in form of giant arcs (see Table 2 (right)) close to the (almost) circular critical curve, or even complete Einstein rings, see e.g., [69] or [70]. Giant arcs and Einstein rings as highly distorted and extended images provide non-local information in the sense that they trace the critical curves in a larger region compared to less distorted and therefore locally more constrained multiple images that are located at larger distances to the critical curves. As noted in Section 2.2, potential biases in gravitational lens modelling decrease for decreasing distance to the model-invariant critical curve, such that Einstein rings are model-independent probes of the total enclosed mass. For static sources, it is advantageous that the intrinsic source properties become negligible with decreasing distance to the critical curve as well. Multiple images of a time-varying source close to the critical curve embedded in a host galaxy that forms giant arcs or an Einstein ring are ideal configurations to constrain $H_0$, see e.g., [14], or to calculate the total enclosed mass with Equation (16). Example applications are given in Section 3.6.

Symmetric multiple-image configurations with a small perturbation also constrain the properties of the perturbing mass. Using lens models, upper bounds on the temperature and the sub-galactic distribution of dark matter can be inferred, see e.g., [71,72,73,74].

In [24], we investigated which properties of a perturber to an axisymmetric mass-density distribution could be determined from the observables shown in Table 2 (right) without a lens model. Due to the underlying axisymmetry, we used polar coordinates $(r,\phi )$ and perturbed the symmetric deflection potential, ${\psi }_{\mathrm{a}}\left(r\right)$, around the Einstein radius $r_{\mathrm{E}}$ with a small deflection potential ${\psi }_{\mathrm{p}}(r,\phi )$, $\left|{\psi }_{\mathrm{p}}(r,\phi )\right|\ll \left|{\psi }_{\mathrm{a}}\left(r\right)\right|$, such that the total deflection potential $\psi (r,\phi )$ reads

$$\psi (r,\phi )=\sum _{n=0}^{\infty }\left(a_n+p_n\left(\phi \right)\right){(r-r_{\mathrm{E}})}^n$$

(29)

with the coefficients $a_n$ and $p_n\left(\phi \right)$ determined as derivatives of the two potentials at $r=r_{\mathrm{E}}$

$$\begin{array}{c}{a}_n={\left.{\displaystyle \frac{1}{n!}}{\displaystyle \frac{{\mathrm{d}}^n{\psi }_{\mathrm{a}}\left(r\right)}{\mathrm{d}{r}^n}}\right|}_{r=r_{\mathrm{E}}},p_n\left(\phi \right)={\left.{\displaystyle \frac{1}{n!}}{\displaystyle \frac{{\partial }^n{\psi }_{\mathrm{p}}(r,\phi )}{\partial r^n}}\right|}_{r=r_{\mathrm{E}}}.\hfill \end{array}$$

(30)

To leading order, the lens equation yields the model-independent local lens properties

$$\begin{array}{cc}{r}_i-r_j\hfill & ={\displaystyle \frac{{\alpha }_{\mathrm{p},r}(r_{\mathrm{E}},{\phi }_i)-{\alpha }_{\mathrm{p},r}(r_{\mathrm{E}},{\phi }_j)}{2(1-{\kappa }_{\mathrm{a}}\left(r_{\mathrm{E}}\right))}},\hfill \end{array}$$

(31)

$$\begin{array}{cc}{\phi }_i-{\phi }_j\hfill & ={\displaystyle \frac{{\alpha }_{\mathrm{p},\phi }(r_{\mathrm{E}},{\phi }_i)}{r_i^2}}-{\displaystyle \frac{{\alpha }_{\mathrm{p},\phi }(r_{\mathrm{E}},{\phi }_j)}{r_j^2}}.\hfill \end{array}$$

(32)

${\kappa }_a\left(r_{\mathrm{E}}\right)$ denotes the convergence of the axisymmetric part of the potential at $r_{\mathrm{E}}$. ${\mathit{\alpha }}_{\mathrm{p}}=({\alpha }_{\mathrm{p},r},{\alpha }_{\mathrm{p},\phi })$ is the deflection angle caused by the perturber. $r_i-r_j$ can be the difference between the observed radial positions of arcs A and B. Alternatively, the difference between the radial positions of the centre-close and centre-far isophote delineating one arc can also be inserted as $r_i-r_j$. Dividing Equation (31) for the radial differences of the delineating isophotes of arc A by the same equation for arc B, we arrive again at Equation (12). Since this procedure maps the two arcs onto each other, it is not surprising that the ratio of Equation (31) of arc A to B is a special case of Equation (12) for lenses with almost axisymmetric mass-density distributions.

Equation (31) implies a degeneracy between the deflection angle of the perturber and the convergence of the axisymmetric lens. As further detailed in [24], this degeneracy can be interpreted as an exact source-position transformation, introduced in [75]. As explained in [27], this degeneracy originates from an a priori arbitrary split of the total deflection potential into a main lens and a perturber. Without additional assumptions or information on the mass-density profile of the axisymmetric lens at $r_{\mathrm{E}}$ or on the deflection angle of the perturber, Equations (31) and (32) are under-constrained. By means of simulations, this was noted in [76]. Applying Equations (31) and (32) to their data, as detailed in [24], the observed degeneracies can be explained.

At best, different models based on a single smooth deflecting mass density or including perturbing inhomogeneities can be ranked with respect to each other in a Bayesian framework, see e.g., [31] or [77]. In this way, the most likely explanation for the perturbed, symmetric multiple-image configuration is determined from a certain class of models. In Section 3.4, we apply Equations (31) and (32) to an example, which is further detailed in [24]. Assuming a specific model for the perturber and the main lens and inserting them into the model-independent equations, we can compare the total mass of the perturber with existing estimates obtained by the Bayesian approach of [31].

2.5. Lensing Distance Ratios in a General Friedmann Universe

SNe are standardisable candles which allow for a data-based distance measure that is independent of any specific parametrisation of the FLRW metric. Replacing the distances based on the FLRW metric (see Equation (1), using Equation (2) as $E\left(z\right)$) by data-based distance measures makes the angular diameter distances in Equation (3) independent of any parametrisation of the FLRW metric. This means that the data-based distances are agnostic about any partition of the total cosmic energy content into radiation, matter, or additional contributions (as done for Equation (2)). Thus, using these SNe-based distance measures, we free our lens-model-independent approach from a parametrisation of the FLRW metric as well. In particular, we avoid setting up particular assumptions about the cosmological constant or potential dark energy models.

In the standardisation of the SNe, an overall distance scale for an ensemble of observed SNe has to be fixed by complementary data. It can either be given by the absolute magnitude of a standard SNe (e.g., from measurements in our local cosmic neighbourhood, see [78]), or by $H_0$ (as determined in our local cosmic neighbourhood, [79], or from the cosmic microwave background, [64]). This implies that data sets of SNe determine $E\left(z\right)$ and leave $H_0$ as a free parameter. The Pantheon sample, as detailed in [41], is the most recent compilation of 1048 unscaled SNe with redshifts $z<2.3$. This redshift range covers a large fraction of the universe, in which gravitational lensing effects occur. Therefore, we determined angular diameter distances based on the Pantheon sample in a flat FLRW metric.

First, we expanded the scale-free, observable luminosity distances ${\tilde{D}}_{\mathrm{L}}(0,z)$ into a complete, orthonormal set of analytic basis functions $b_j$, $j=0,\dots ,n_{\mathrm{b}}$, the so-called Einstein–de-Sitter polynomial basis as introduced in [42]

$${\tilde{D}}_{\mathrm{L}}(0,z)\equiv {\displaystyle \frac{H_0}{c}}{D}_{\mathrm{L}}(0,z)=(1+z)\underset{0}{\overset{z}{\int }}{\displaystyle \frac{\mathrm{d}x}{E\left(x\right)}}\equiv \sum _{j=0}^{n_{\mathrm{b}}}{c}_j{b}_j\left(z\right).$$

(33)

From the Pantheon sample of SNe, four coefficients for the first four basis functions can be significantly determined given the measurement precision and the covariances in the sample. Inserting these scale-free luminosity distances represented by the Einstein–de-Sitter basis into the equation

$${\displaystyle \frac{\mathrm{d}{\tilde{D}}_{\mathrm{L}}(0,z)}{\mathrm{d}z}}={\displaystyle \frac{{\tilde{D}}_{\mathrm{L}}(0,z)}{1+z}}+{\displaystyle \frac{1+z}{E\left(z\right)}}$$

(34)

and requiring that $E\left(0\right)=1$, we obtained a data-based cosmic expansion function, $E\left(z\right)$, as first noted by [80]. Subsequently, we inserted this $E\left(z\right)$ into Equation (1) together with a value for $H_0$, e.g., as constrained by [64] or [79] to obtain a distance measure.

Replacing the model-based angular diameter distances in Equation (3) by the data-based ones, we obtain the data-based lensing distance ratio, as detailed in [26]. Leaving $H_0$ as a free parameter, we can use time delay differences to constrain its value, if the difference of the Fermat potentials between the image positions is known. Alternatively, with a fixed value for $H_0$, the time delay differences between two images constrain the difference between their Fermat potentials (see Equation (16)).

In Section 3.5, we compare the relative precision of the lensing distance ratio of angular diameter distances as obtained from the Pantheon sample with the lensing distance ratio based on the angular diameter distances of the Friedmann model as established by [64]. Alternatively, given the Fermat potential difference between two images, a measured time delay difference between these images infers a value for $H_0$. Constraints on its precision are estimated in Section 3.6.

3. Results

In this section, we discuss example applications for our approach, theoretically outlined in Section 2, to observational data. In Section 3.2, we show the application of the approach of Section 2.3 for five resolved images of a background galaxy behind the galaxy cluster CL0024 observed with the Hubble Space Telescope (HST). Further details about this example can be found in [47]. Section 3.3 subsequently treats the application of our approach of Section 2.3 to four unresolved, extended images of a quasar background source behind the spiral galaxy B0128, as further detailed in [48]. For this application, very-long-baseline interferometry (VLBI) observations in the radio bands by [63], were employed. In Section 3.4, we discuss the degeneracies detailed in Section 2.4 that arise for the two giant arcs of a background galaxy caused by the early-type galaxy B1938+666. We review the precision up to which the lensing distance ratio in Equation (6) can be determined from the Pantheon sample of SNe type Ia, [41], in Section 3.5. The details to set up the angular diameter distance measure from the sample are found in [26]. Section 3.6 sets up theoretical limits on the precision to which $H_0$ can be determined by Equation (16). This section is based on the findings in [27] and supported by data from quasar time delay monitoring, see e.g., the paper series of [81], and FRB observations, discussed in [49].

In [22,23], we used simulation data based on simple lens models to test the accuracy of the approximations detailed in Section 2.3 and Section 2.4. The relative image distances and the parameters employed for the lens models considered in these simulations were extreme cases compared to observations and lens model parameters as inferred by observations. Despite this, we obtained inaccuracies of about 10% and mostly less for the ratios of third order derivatives at folds and cusps. Most of the $f_{ij}$ and $g\left({\mathit{x}}_i\right)$, $i,j=A,B,C$, lay within the 1$-\sigma$ confidence bounds around the true value. As a consequence, we expect our approach to yield more accurate results in average-case scenarios considered in this section, so that the local lens properties and pathways of the critical curves in the vicinity of multiple images can be constrained more tightly.

3.1. Data Preprocessing

Before any observables of Table 1 and Table 2 can be extracted, the observational data have to be preprocessed. Effects of the data acquisition, such as shot noise or the impact of the point spread function (PSF), have to be taken into account. In addition, astrophysical effects like stray light from or partial occlusion by neighbouring objects must be handled, see e.g., [82] or [83] for a detailed explanation of all effects.

In order to keep our approach model-independent, the data preprocessing should not rely on any model. Data processing methods that employ a model-based separation of the foreground deflecting galaxy from the multiple images for galaxy-scale lensing can be found in [84,85,86]. Furthermore, we separate the extraction of the observables from the application of our approach to these data. In this way, the lensing analysis pipeline is kept as flexible as possible to use our approach for data of various ground-based and space-based instruments and the impact of different preprocessing methods on the local lens properties can be easily compared.

Since the PSF of the HST is of the order of one pixel and thus much smaller than the extensions of the multiple images, effects of the convolution with the PSF are negligible. For Section 3.2, the reference points for in the multiple images were identified manually directly from the FITS-file provided in the HST archive. The same applies to the observables in Section 3.4. The observables in the radio band, including the quadrupoles of the multiple images as required for the application in Section 3.3, were already measured by [63] using the National Radio Astronomy Observatory Astronomical Image Processing Software package AIPS.

3.2. Comparison to Lens Modelling Approaches for the Cluster-Scale Lens CL0024

CL0024+1654, abbreviated as CL0024, is a massive galaxy cluster at redshift $z_{\mathrm{l}}=0.39$. It consists of two components merging along the line of sight (as detailed e.g., in [87,88,89]) and it magnifies a blue galaxy at redshift $z_{\mathrm{s}}$ = 1.675 into five highly-distorted, well-resolved images, as first described by [90] (see images 1.1 to 1.5 in the yellow ellipses in Figure 3). In Figure 3, the yellow points in the ellipses mark the coordinates of the multiple images when they are treated as point-like constraints of global lens reconstructions. In the details at the sides, the six reference points that can be identified in all multiple images are marked by red circles. We extracted the reference points in all multiple images as observables from an HST ACS/WFC picture in the F475W band2. Using these reference points, we applied our approach of Section 2.3 to these multiple images to constrain the local lens properties in Equation (12). World-coordinate-system coordinates of all points can be found in [47]. The work in Ref. [91] predicted ten additional multiple-image sets (MISs) from sources with photometrically determined mean redshifts $z_{\mathrm{s}}\in \left[0.25,4.16\right]$ (see blue points 2.1 to 11.4 in Figure 3, coordinates of these points can be found in [91]). From these additional ten MISs, we employed MIS 3, 4, 5, 8, and 10 as point-like multiple-image constraints to set up global mass-density reconstructions. We selected these MISs because they are scattered broadly in the region covered by the cluster on the sky so that we expected well-constrained global mass-density reconstructions. In addition, these MISs are located close to the resolved multiple images of set 1, so that a potential influence of constraints from neighbouring MISs on the local lens properties at MIS 1 should be observable.

As global lens reconstruction methods, we chose the parametric reconstruction method Lenstool, described in detail in [36] and references therein, and the free-form reconstruction approach Grale, as detailed in [34]. Appendix C contains a brief characterisation of both algorithms in terms of their peculiarities that are relevant for our comparison. Both approaches used the same data and yielded maps of the convergence, $\kappa \left(\mathit{x}\right)$, and the shear, $\gamma \left(\mathit{x}\right)$, in the region covered by the multiple images. We calculated the local lens properties of Equation (12) from the $\kappa \left(\mathit{x}\right)$- and $\gamma \left(\mathit{x}\right)$-maps of the global reconstruction approaches at the positions of the multiple images of MIS 1 marked by the yellow points in Figure 3. Then, we compared these local lens properties obtained by the global methods with the local lens properties at the same position obtained by our model-independent approach detailed in Section 2.3.

Implementing different constraints in the Lenstool and Grale reconstructions and comparing the results among each other and with the results of our approach, we investigated

• whether neighbouring MISs influence the reconstruction of local lens properties at the image positions of MIS 1 in Lenstool and Grale and thus introduce correlations between MISs that are not accounted for in our approach,
• whether the local lens properties at MIS 1 as obtained by our local approach and by the two global reconstructions coincide, or if image distortions beyond leading order play a significant role (apart from potential correlations between MISs mentioned under 1.),
• whether the light-traces-mass (LTM) assumption used in Lenstool is corroborated and, if so, which mass-to-light ratio (MLR) is favoured for the cluster member galaxies, and
• whether the small-scale fine-tuning of the mass density implemented in Grale overfits the free-form lens reconstruction to the multiple-image constraints.

The local lens properties that we obtained showed broad confidence bounds for the multiple images close to the $\kappa \left(\mathit{x}\right)=1$ isocontour, which was to be expected from Equation (12). The tightest confidence bounds were obtained using all six reference points which delineate the largest area in each multiple image in accordance with the statements made in Section 2.3. Relative sizes of the confidence bounds with respect to the quantities in Equation (12) ranged from 8% to 150% with an average relative size of 9.4%. The implementation of our approach took less than half a second to determine the quantities in Equation (12). Assuring that all reconstructions are of similar quality, we found a high degree of coincidence in the local lens properties obtained by all three approaches. From this result, we concluded that

• non-local influences of neighbouring MISs are negligible at the current precision of the quantities in Equation (12). Thus, the local lens properties are mainly derived from the local constraints of the respective multiple images. Hence, we can reconstruct the morphology of the source up to an overall scale as detailed in Section 2.3. The back-projections of all images to the source plane by their model-independently reconstructed $A\left({\mathit{x}}_i\right)$ and their pixel-wise averaged source are shown in Figure 4 (for comparison with model-based source reconstructions, see [90,92]; a physical analysis of the Lenstool-reconstructed source can be found in [93,94]).
• the leading-order local lens properties of our model-independent approach mostly agree to local lens properties as obtained by Lenstool and Grale within their 1$-\sigma$ confidence bounds.
• the LTM assumption is corroborated by the high degree of agreement of the local lens properties as obtained by Lenstool and the other two approaches that do not make any assumptions about the relation between dark and luminous matter. A comparison of the local lens properties of all approaches favours a constant MLR for the brightest cluster member galaxies over a non-constant MLR. This result was also found in [89,95]. A non-constant MLR that reproduces the fundamental plane (see [96] and references therein) yields a higher degree of agreement between the redshifts as estimated by Lenstool and the photometric redshifts of the additional MISs. This indicates that spectroscopic redshift observations are necessary to further constrain the redshifts and the MLR.
• the small-scale fine-tuning of the mass density in Grale does not significantly change the local lens properties, nor does it tighten the confidence bounds. As it does not significantly alter the pathways of the critical curves, either, it can be omitted for run-time optimisation.

Consequently, assuming that our results for CL0024 are representative, we demonstrated that a comparison of the local lens properties between the three conceptually different reconstruction ansatzes allows for corroborating or refuting assumptions of global reconstruction methods. Furthermore, we showed that the source object can also be reconstructed by a back-projection of the multiple images with the local distortion matrix. Hence, a full reconstruction of the deflecting mass density is not necessary to study the morphology of individual background objects. The back-projections shown in Figure 4 are the most general leading-order source reconstructions directly based on the standard lensing formalism and the observables in the multiple images. In addition, as our source reconstruction does not assume a lens model, studying star formation and stellar populations in faint and distant, otherwise unobservable galaxies in the early universe (even at redshifts as high as 9–10) without a potential lens-model bias becomes possible. The source reconstruction is still subject to the local MSD of this MIS (see [25] for this local generalisation of the global MSD as defined by [68]). However, the latter can be easily broken by rescaling the distortion matrices, as soon as the actual scaling is known, e.g., from time delay difference observations.

3.3. Transfer of the Method to Unresolved Multiple Images for the Galaxy-Scale Lens B0128

In the second example which is discussed in detail in [48], we demonstrate that our approach can also be applied to galaxy-scale gravitational lenses with multiple-image configurations as observed in the radio bands. Figure 5 (left) shows the configuration of four images of a quasar at $z_{\mathrm{s}}=3.124$ created by the lenticular or late-type galaxy B0128+437, B0128 for short. Most probably, its redshift is $z_{\mathrm{l}}=1.145$, [97]. The observation was made in the Multi-Element Radio Linked Interferometer Network (MERLIN) 5 GHz band, as detailed in [98]. The multiple images denoted by A, C, and D in Figure 5 are resolved into three clearly identifiable sub-components in VLBI observations in the 5 GHz and 8.4 GHz bands (see details on the left and right of the MERLIN observation), see [63]. Image B is most likely scatter-broadened. The coordinates of all observational data are summarised in [48].

Analogously to the example in Section 3.2, we used the example to compare the local lens properties obtained by our approach according to Section 2.3 with the local lens properties obtained by the global free-form lens reconstruction approach PixeLens, [29]. We investigated whether the lens reconstruction by previous global parametric lens reconstruction methods was oversimplified because no global parametric lens-modelling approach could explain the sub-components in images A, C, and D to their sub-milli-arcsecond precision, while the configuration of the four images shown in Figure 5 (left, large figure) could be described by a singular isothermal elliptical (SIE) mass-density profile with external shear, as detailed in [63]. According to the simulations performed in [73], it was found that the anomalous flux ratio between images A and B is not likely to be caused by small-scale dark-matter clumps. Therefore, [73] concluded that the scatter-broadening in image B and an oversimplified lens reconstruction could cause the discrepancies between the fits of smooth symmetric models and the observations.

Using the sub-components of images A, C, and D as reference points, we calculated the local lens properties according to Equation (12). In addition, as the sub-components 1 and 3 in these three images were characterised by their quadrupole moment, we also calculated the local lens properties for both sub-components using the centre of light and the end points of the semi-minor and the semi-major axis as reference points (see Section 2.1). Figure 5 (right) sketches the transformations on image and sub-component scale to obtain the respective local lens properties. Due to the almost parallel alignment of the sub-components, the local lens properties determined from the sub-components as reference points were subject to broad confidence bounds. Their relative sizes compared to the f- and g-values range from 8% to over 1000%. Comparing these confidence bounds to the ones obtained from the quadrupole moments, which mostly have relative sizes on the order of 20%, we found that the latter are smaller because the vectors spanned by the reference points are orthogonal to each other. Half of the local lens properties at the positions of the two sub-components agree within their confidence bounds, so that it remains undecided if the distortion matrix can be assumed constant over the area of all sub-components or not. Only a weak upper limit to the scale of potential higher order perturbations, like gradients, in the surface mass density at the position of the images A, C, and D may be inferred from this result. However, the result confirms the hypothesis of [73] that the global parametric lens reconstructions were indeed oversimplified. A comparison with our PixeLens lens reconstruction and an additional symmetry analysis according to the method based on relative opening angles between the multiple images as detailed in [99] both further corroborate the hypothesis. All results concordantly hint at a mass-density distribution that deviates from the simple elliptical models usually assumed.

3.4. Degeneracies Arising in the Detection of Perturbers for the Galaxy-Scale Lens B1938+666

B1938+666, detected in the Jodrell Bank Very Large Array Astrometric Survey, [100], is an early-type gravitational lens at redshift $z_{\mathrm{l}}=0.88$, [101]. The background source, a quasar embedded in its host galaxy at redshift $z_{\mathrm{s}}=2.06$, [102], is assumed to consist of two components, [103]. Multifrequency Very Large Array (VLA), MERLIN and VLBI radio observations revealed that the first source component is mapped into four images and the second into two images, [104]. As observed in [105] by HST observations and later on confirmed by observations with the NIRC2 camera on the Keck II Telescope, using the Laser Guide Star Adaptive Optics System in [86], the source galaxy forms an Einstein ring visible in the near infrared band. Figure 6 shows the observational data in the optical and near infra-red.

Perturbations in the symmetric brightness distribution of the Einstein ring were used to infer properties of a potential perturber, like its location and its mass in [86,106]. As there is no evidence for luminous perturbers equal and above the mass estimate determined in [106] (see below), [86], the existence and the properties of a potential perturber were solely inferred from the observable brightness profiles of the giant arcs. In [86,106], the Bayesian lens-reconstruction algorithm of [31] was employed to infer the most likely properties of a perturber. Using the HST data, [106] arrived at

$$m_{\mathrm{p}}=\left(1.9\pm 0.1\right)·{10}^8{\mathrm{M}}_{\odot }$$

(35)

assuming that the perturber is located close to the critical curve in the lens plane and the projected distance between the symmetry centre of the main lens and the perturber in the lens plane is the total separation. Hence, Equation (35) is a lower limit on the perturbing mass. As also confirmed in [86], the mass estimates derived from the H- and K’-band Keck data are consistent with the value in Equation (35). Furthermore, [86] followed up on these results with Keck observations and compared the precision to which the parameters of the most likely lens reconstruction can be determined to the precision obtainable when using the HST observations. For this highly symmetric case, the Keck observations were found to yield more precise parameter estimates. Nevertheless, due to the high symmetry of the main lens and the featureless, i.e., unresolved, giant arcs, degeneracies between the lens model parameters still remained large. This implies that significantly different models can fit the data equally well. To arrive at this result and test the robustness of the obtained properties of the perturber, [86,106] had to set up different models and probe the multi-dimensional parameter space of each model within given prior ranges.

Employing the approach outlined in Section 2.4.3 to estimate the mass of a potential perturber, we first fitted a circle to the Einstein ring of the HST data and obtained $r_{\mathrm{E}}={0.45}^{″}$. This agrees well with $r_{\mathrm{E}}={0.465}^{″}$, as determined in [105]. In [106], a lower value of $3.39\mathrm{kpc}$ was used, which corresponds to ${0.444}^{″}$ at $z_{\mathrm{l}}=0.88$, assuming a flat $\Lambda$CDM with ${\Omega }_m=0.25$, ${\Omega }_{\Lambda }=0.75$, $H_0=73$ km/s/Mpc. Fitting two circles to the upper westward and lower eastward arc, denoted as image A and B, respectively (see Figure 6 (left)), $r_{\mathrm{A}}={(0.50\pm 0.1)}^{″}$ and $r_{\mathrm{B}}={(0.41\pm 0.1)}^{″}$ were obtained. Since the signal-to-noise ratio is to low to resolve brightness features in the arcs like the ones observable in the multiple images in CL0024, Equation (32) cannot be applied to this case. Assuming that the perturber is a satellite galaxy with a singular isothermal sphere (SIS) as mass-density profile close to a SIS main lens, we inserted the ${\mathit{\alpha }}_{\mathrm{p}}\left(\mathit{x}\right)$ and ${\kappa }_{\mathrm{a}}\left(r_{\mathrm{E}}\right)$ for a SIS (see e.g., [50]) into Equation (31). Solving for the Einstein radius of the perturber, we obtained

$$m_{\mathrm{p}}=\left({5.89}_{-2.30}^{+2.91}\right)·{10}^8{\mathrm{M}}_{\odot }$$

(36)

as a mass estimate for the perturber, if it is optimally aligned, i.e., located on the line connecting the centre of light of the giant arc and the symmetry centre of the lens in the lens plane. Comparing Equations (35) and (36), we find that both are of the same order of magnitude.

Summarising the results for this example, as further detailed in [24], we showed that our simple Equations (31) and (32) to relate lens properties with observables allow us to infer lens properties that are in agreement with sophisticated lens reconstruction algorithms like [31], which implement more realistic classes of main lens and perturbing mass-density profiles. In addition, Equations (31) and (32) clearly show the most general transformations between the perturber and the main, axisymmetric lens that leave the observables of the giant-arc brightness profiles invariant. Similar equations could be analogously derived for main-lens mass profiles with other symmetries (e.g., ellipses). However, as also noted in [71], without additional non-lensing information on the properties of the perturber, like its location (on the sky and also along the line of sight) and its brightness profile, the properties of the potential perturbers remain highly degenerate.

3.5. Lensing Distance Ratios from the Pantheon Sample of Supernovae

As outlined in Section 2.5 and detailed in [26], we used the Pantheon sample of SNe type Ia of [41] and $H_0$ from [107] to set up data-based angular diameter distances for the lensing distance ratio in Equation (6). We abbreviate the lensing distance ratio by

$$D(z_{\mathrm{l}},z_{\mathrm{s}})\equiv {\displaystyle \frac{D_{\mathrm{A}}(0,z_{\mathrm{s}})D_{\mathrm{A}}(0,z_{\mathrm{l}})}{D_{\mathrm{A}}(z_{\mathrm{l}},z_{\mathrm{s}})}}.$$

(37)

To investigate the precision of the data-based $D(z_{\mathrm{l}},z_{\mathrm{s}})$ as established with the Pantheon sample, we implemented a Monte Carlo simulation of 1000 Pantheon-like data sets, each containing 1048 simulated, scale-free luminosity distances at the same redshifts as the Pantheon sample. Each of the 1048 data points was generated by drawing a random number from a Gauss distribution around the measured scale-free luminosity distance from the Pantheon sample with a standard deviation given by its measurement uncertainty. With these simulated data sets, the confidence bounds on $D(z_{\mathrm{l}},z_{\mathrm{s}})$ can be determined for any lens and source redshift in the range of the Pantheon sample ($z\le 2.3$), given a value for $H_0$. As confidence bounds, we calculated the standard deviation around the mean $D(z_{\mathrm{l}},z_{\mathrm{s}})$ and the confidence bounds based on percentiles, given by the $68\%$, $95\%$, and $99\%$ confidence levels. In principle, the precision of $D(z_{\mathrm{l}},z_{\mathrm{s}})$ also depends on $H_0$ and its uncertainty. However, due to the still ongoing debate over the tension in the value for $H_0$, [64,79], we only considered the relative precision of $D(z_{\mathrm{l}},z_{\mathrm{s}})$ and did not include the imprecision of $H_0$ in its uncertainty.

To compare the precision of the data-based $D(z_{\mathrm{l}},z_{\mathrm{s}})$ with the one obtained from the $\Lambda$CDM standard cosmological model as determined by [64], we developed a second Monte Carlo simulation. We drew 1000 different cosmological models from a Gaussian distribution around the cosmological parameters as determined by [107], ${\Omega }_m=(0.3089\pm 0.0062)$, ${\Omega }_{\Lambda }=(0.6911\pm 0.0062)$, with their measurement uncertainty as standard deviation3. For each of the 1000 cosmological models, we determined model-based luminosity distances at the redshifts of the Pantheon sample, given $H_0=67.74$ km/s/Mpc. With this catalogue of 1000 data sets containing 1048 data points each, we analogously calculated the relative precision of $D(z_{\mathrm{l}},z_{\mathrm{s}})$, leaving out the imprecision of $H_0$.

As $D(z_{\mathrm{l}},z_{\mathrm{s}})$ is dependent on $z_{\mathrm{l}}$ and $z_{\mathrm{s}}$, the best visualisation is obtained when fixing $z_{\mathrm{l}}$ and plotting the relative imprecision of $D(z_{\mathrm{l}},z_{\mathrm{s}})$ for all $z_{\mathrm{s}}$, as shown in Figure 7. In addition to the comparison at $z_{\mathrm{l}}=0.5$ established in [26], we also show the relative imprecisions for $D(z_{\mathrm{l}},z_{\mathrm{s}})$ at $z_{\mathrm{l}}=0.25$ and $z_{\mathrm{l}}=1.0$ as given by the confidence bounds of the simulated data sets.

By comparing the left- and the right-hand side of Figure 7, we find that there are one to two orders of magnitude difference in the imprecision for the data-based (left) and Planck-model-based lensing distance ratio (right). The precision to which the data-based lensing distance ratio is determined is constrained by the measurement uncertainties of the SNe and the choice of the set of basis functions. Among the analytic sets of basis functions that we tested in [26], the Einstein–de-Sitter basis with four basis functions was found to yield the tightest confidence bounds and the most accurate reconstruction of a flat $\Lambda$CDM model. Its inaccuracies, determined in a simulation of a flat $\Lambda$CDM model, are contained within the confidence bounds of the imprecisions to avoid biases. In addition, the number of basis functions was determined in a ${\chi }^2$-parameter estimation with a reduced ${\chi }^2$-value of 0.9 to avoid overfitting. Details about the quality test of all basis functions are given in [26]. With further SNe data or a higher measurement precision, the confidence bounds can be further reduced and the number of basis functions increased.

The confidence bounds of the lensing distance ratio in a Planck-model-based cosmology are tighter and negligible compared to other sources of uncertainty in the time-delay equation (see Section 3.6). However, the tighter confidence bounds come at the cost that the model may be currently incomplete or it may become biased for an increasing amount of future data. While the extensions of the $\Lambda$CDM model that can be included to amend the bias have a direct physical interpretation, e.g., an additional curvature term or any forms of dark energy, it is not clear which extension approach should be followed. Contrary to that, an extension of the Einstein–de-Sitter basis with four basis functions is straightforward. Assigning the basis functions a physical interpretation is the difficult part in this approach.

3.6. Usage of the Time-Delay Equation to Infer $H_0$

As shown in [27], $H_0$ is uniquely determined by a time-delay-difference measurement in a FLRW background cosmology if the total integrated deflecting mass along the two light paths is known. However, observations to tightly constrain the total mass are rare, especially if we include the small-scale perturbing masses not associated with the most massive deflecting object along a light path. In addition, parametrisations of the FLRW background cosmology are still subject to debates about potentially required extensions to the current concordance $\Lambda$CDM standard model (see Section 3.5). Therefore, current approaches to determine $H_0$ first solve the time-delay equation

$${\tau }_{ij}={\displaystyle \frac{(1+z_{\mathrm{l}})}{c}}D(z_{\mathrm{l}},z_{\mathrm{s}})\left(\varphi ({\mathit{x}}_i,\psi )-\varphi ({\mathit{x}}_j,\psi )\right),$$

(38)

for the lensing distance ratio $D(z_{\mathrm{l}},z_{\mathrm{s}})$, Equation (37). For ${\tau }_{ij}$, the observed time-delay difference between light rays arriving from the multiple images i and j is inserted and ${\varphi }_{ij}\equiv \varphi ({\mathit{x}}_i,\psi )-\varphi ({\mathit{x}}_j,\psi )$ is derived from a lens model. Since $D(z_{\mathrm{l}},z_{\mathrm{s}})$ contains the full information about the cosmological geometry of the multiple-image configuration, subsequently inferring $H_0$ from it conveniently allows to investigate different parametrisations of the FLRW background cosmology without repeating the lens reconstruction.

In practice, reconstructions of ${\varphi }_{ij}$ are subject to degeneracies. In our approach as outlined in Section 2, a local MSD at ${\mathit{x}}_i$ and ${\mathit{x}}_j$ remains to be broken, i.e., ${\varphi }_{ij}$ between the two light rays can be scaled by a factor and $H_0$ can be rescaled accordingly to keep the observed time delay difference invariant. As, by definition, $H_0$ is linked to the overall energy density of our universe today, this freedom of rescaling can be interpreted as redefining the split between the cosmological background energy density and the deflecting mass density of gravitational lensing on top of it. Other approaches usually split the total mass along the line of sight into a main deflecting mass and small-scale perturbers along the line of sight. In this case, the factor $\lambda$ appearing in an MSD consists of two parts: one part that can be associated with the small-scale perturbers and neighbouring masses of the main deflector, which is usually called external convergence ${\kappa }_{\mathrm{ext}}$, see e.g., [77]. The second part of the factor is the remaining MSD far away from the lens, in the region where only the embedding cosmological background density remains. This explains the argument stated in [76] that the MSD cannot be solely caused by ${\kappa }_{\mathrm{ext}}$.

On the galaxy-cluster scale, [108] systematically investigated the relative precision to which $H_0$ can be determined from the well-resolved system of multiple images of the host galaxy of SN Refsdal and a plethora of multiple-image configurations from different background sources at different redshifts alone and in combination with different non-lensing observables and assumptions. They found that a combination of a flexible free-form with a consistent parametric lens reconstruction is capable of alleviating all occurring degeneracies to obtain $H_0$ to a relative imprecision of less than 10%. The precision of the inferred $H_0$-value greatly depends on the configuration of the multiple-image sets, the measured time-delay differences but even more on the additional assumptions being made.

Galaxies as lenses are usually described by comparably simple models, yet the high symmetry and the sparsity of multiple-image configurations also require additional non-lensing observables, like the velocity dispersions along the line of sight, to be included in the model, see e.g., [109]. On both scales, potential biases due to the choice of a specific model (class) and degeneracies between the model parameters remain, see [27] and references therein. To further decrease the imprecision in $H_0$, several multiple-image configurations with measured time delay differences can be jointly analysed, as done e.g., in [110,111], or [112].

Having obtained a value for $D(z_{\mathrm{l}},z_{\mathrm{s}})$, the definition of model-based angular diameter distances, Equation (1) usually using Equation (2) (or $E\left(z\right)$ for an extension of the $\Lambda$CDM model), is inserted to relate $D(z_{\mathrm{l}},z_{\mathrm{s}})$ to $H_0$. To determine $H_0$ from the time-delay equation in approaches based on Bayesian parameter inference, see e.g., [15,40], parameters of the expansion function are either fixed, as in Section 3.5, or considered as nuisance parameters and are marginalised over. Marginalising implies that the imprecision of $D(z_{\mathrm{l}},z_{\mathrm{s}})$ for Planck-model-based distances, is larger than the ones shown in Figure 7. As found in [40], the currently lowest imprecision for $H_0$ inferred from a single galaxy-scale lens is 3–4% at 68% confidence level. On galaxy-cluster level, the lowest imprecision for $H_0$ is 4–5% at 68% confidence level obtained from the multiple images of SN Refsdal in MACS 1149, [15].

To estimate the precision to which $H_0$ can be theoretically determined with a single time-delay-difference measurement when using the Pantheon-data-based distances, we rearrange Equation (38), such that

$$H_0={\displaystyle \frac{\tilde{\Gamma }}{{\tau }_{ij}}}{\varphi }_{ij}$$

(39)

with $\tilde{\Gamma }=H_0\Gamma$ (see Equation (6)). Then, we determine the relative imprecision of $H_0$ as

$${\displaystyle \frac{\Delta H_0}{H_0}}=\sqrt{{\left({\displaystyle \frac{\Delta D(z_{\mathrm{l}},z_{\mathrm{s}})}{D(z_{\mathrm{l}},z_{\mathrm{s}})}}\right)}^2+{\left({\displaystyle \frac{\Delta {\tau }_{ij}}{{\tau }_{ij}}}\right)}^2+{\left({\displaystyle \frac{\Delta {\varphi }_{ij}}{{\varphi }_{ij}}}\right)}^2}$$

(40)

assuming that the relative imprecision in $\tilde{\Gamma }$ is dominated by the imprecision in $D(z_{\mathrm{l}},z_{\mathrm{s}})$ and the measurement uncertainty in $z_{\mathrm{l}}$ is negligible. From Figure 7, we read off the relative imprecision of the data-based lensing distance ratio of about 2% at 68% confidence level. Constraints on the second term for different time-varying sources in galaxy-cluster-scale lensing can be found in [49]. For a quasar, we found a relative imprecision in the time delay difference of 1.5%. In [113] for quasars multiply-imaged by galaxy-scale lenses, it is stated that 3% imprecision are obtained by the best algorithms. For FRBs, the second term is negligible with absolute measurement imprecisions on the order of milli-seconds compared to time delay differences of the order of days or longer.

However, by far the largest source of uncertainties is the third term, see e.g., [108,114,115]. A systematic investigation for lens-modelling accuracy and precision on galaxy scale is currently being pursued by [12]. The work in Ref. [112] gives estimates for the accuracy and precision of models potentially achievable with HST in the 160 W band of the Wide Field Camera 3 for exposure times of about 10,000 s and an astrometric precision of the order of sub-milliarcseconds in the angular positions, assuming the time-varying source is an FRB. Including the line-of-sight external convergence to the lens modelling, their optimum estimate of relative imprecisions for the difference in the Fermat potential is about 2–3%. Estimating an imprecision of 5% for the difference of the Fermat potential4 for currently available data, we obtain

$$\begin{array}{cc}{\left(\frac{\Delta H_0}{H_0}\right)}_{\mathrm{cluster}}\hfill & =\sqrt{{\left(0.02\right)}^2+{\left(0.015\right)}^2+{\left(0.05\right)}^2}\approx 5.6\%,\hfill \end{array}$$

(41)

$$\begin{array}{cc}{\left(\frac{\Delta H_0}{H_0}\right)}_{\mathrm{galaxy}}\hfill & =\sqrt{{\left(0.02\right)}^2+{\left(0.030\right)}^2+{\left(0.05\right)}^2}\approx 6.1\%.\hfill \end{array}$$

(42)

At the current measurement and lens-reconstruction precision, the imprecision of the Pantheon-data-based angular diameter distances is of the same order of magnitude as the measurement uncertainty of the time-delay-difference and smaller than the imprecision caused by the lens reconstruction. Hence, the cost for avoiding a parametrisation of $E\left(z\right)$ in terms of specific ${\Omega }_i$ is an increase in the imprecision in $H_0$ up to a factor of two for galaxy-scale lenses and by a factor of 1.4 for galaxy-cluster-scale lenses. With the currently increasing number of observations of multiple-image configurations with measured time-delay-differences, this disadvantage can be compensated with a joint analysis of several configurations.

Vice versa, solving Equation (38) for ${\varphi }_{ij}$ for a given $H_0$, we found that the differences in the Fermat potential can be determined up to 2.7–3.7%, if $H_0$ can be measured with a relative imprecision of 1%. We showed in [49] by means of a simulation that even a single time-delay-difference observation can greatly increase the accuracy in the global lens reconstruction of a galaxy-cluster-scale lens, as it fixes the global MSD.

4. Discussion and Conclusions

Thus far, constraints from observables of multiple images on the total deflecting mass-density distribution only sparsely cover the area of the latter. Attempting to reconstruct the total deflecting mass-density distribution in the entire lensing region on the sky is thus an under-constrained optimisation problem. To solve it, additional assumptions are usually employed which capture the global shape and the overall smoothness of the deflecting mass-density distribution. One way to capture the shape uses physically motivated parametric models of the mass-density profile, like a Navarro–Frenk–White profile, or a superposition of parametric mass-density profiles. Another common way decomposes the deflecting mass density into a set of orthonormal basis functions and determines the weights of the leading basis functions by means of the observables. In both methodological approaches, a regularisation has to be introduced to avoid overfitting. The regularisation also sets the smoothness of the mass-density distribution. In Section 1, we gave examples of such approaches for galaxy-scale and galaxy-cluster-scale lenses.

Contrary to these global ansatzes that treat the reconstruction of the deflecting mass-density distribution as a global model fit and intertwine data-based evidence and model assumptions, we develop an approach to locally characterise strong gravitational lenses solely by observables of the multiple images that they create. Our approach exploits the symmetries and correlations between the observables in the multiple images and thereby avoids introducing specific lens models. During its development, we were able to mathematically derive all degeneracies, like the mass-sheet degeneracy, treating the single-lens-plane gravitational-lensing formalism as partial differential equations. Subsequently, we found a physical interpretation for the respective invariance transformations that explains why they appear.

As our method does not reconstruct the global underlying mass density, it can be applied to any multiple-image configuration on any lens scale in the same way. Thus far, we have employed relative distances between multiple images, features in the brightness profiles of resolved multiple images, multipole moments of the brightness profiles up to the quadrupole of unresolved multiple images, and time delay differences, if the source has a time-varying component. From these observables, we calculate local lens properties that are not subject to any degeneracies, if time delay differences in a fixed background metric are used. For static sources, the local equivalent of a mass-sheet degeneracy remains to be broken. This is a consequence of setting up the standard gravitational lensing formalism as a theory of light deflection caused by inhomogeneities on top of a cosmic background metric and only observing angular quantities on the celestial sphere.

To leading order, our approach determines ratios of deflecting mass densities, the reduced shear, and the linear and parabolic approximations of the critical curves in the vicinity of the multiple images. We showed the respective equations in Section 2. From simulations, we found inaccuracies of the order of and less than 10% for extreme cases of large deflecting masses and highly elliptical density profiles. For realistic deflecting mass densities, we expect them to be lower and in the range of a few percent, which will be investigated in N-body simulated gravitational lenses. Imprecisions for some observational cases were found to be much larger, but the imprecision greatly depends on the type of multiple image and the positions of the multiple images with respect to the scaled mass density isocontour $\kappa \left(\mathit{x}\right)=1$ or to the isocontour of ${\varphi }_{11}\left(\mathit{x}\right)=0$, depending on the parametrisation that is chosen.

The set of local lens properties is the maximum leading-order local information in which all lens models coincide and that is directly constrained by the observables in the multiple images in the standard single-lens-plane gravitational-lensing formalism. Formulated differently, our model-independent approach partitions the lensing region on the sky into areas where local lens properties are directly inferred from observables, and the remaining areas where global lens reconstructions are forced to make assumptions about the mass-density distribution because no constraints from multiple images are available. With the continuously increasing number of observed multiple-image configurations in galaxy-cluster-scale lenses, areas where global assumptions are necessary are shrinking. As a consequence, the confidence bounds on the local lens properties of our approach will tighten. Hence, our data-driven approach will become increasingly useful to determine the mass-density distribution of galaxy clusters. As multiple images of different sources rarely overlap, the local lens properties for newly discovered multiple-image sets can be easily added to the already known set of multiple images and their local lens characterisations. It is particularly advantageous for galaxy-cluster-scale lenses with their highly complex mass-density distribution to solicit properties of the deflecting mass distribution at certain positions that do not rely on global symmetry assumptions and that do not jointly utilise the information of several multiple-image configurations.

On galaxy-scale, the increasing amount of observed resolved multiple images of quasars in several radio bands reveals a more asymmetric mass-density distribution of their lenses than parametric, symmetric lens models can provide, so that the local lens properties determined by our approach are useful to sketch the asymmetries in the mass-density distribution. In addition, with increasing measurement precision to extract the shapes of the individual parts of the quasar, applying our approach on the scale of a quasar-image and on the scale of the individual parts of a quasar-image, small-scale lens properties may become detectable.

Beyond constraining local lens properties, comparing the model-independent approach to global lens reconstructions corroborates or refutes the global assumptions that are employed in the lens reconstructions—for instance, the hypothesis of whether luminous matter traces dark matter can be checked. Alternatively, our local lens properties can be integrated as constraints in the global lens reconstructions, so that global lens reconstructions can be envisioned that extrapolate the local lens properties, i.e., the information that all models agree upon, into areas where no multiple images are observed. Compared to existing approaches that directly intertwine observational evidence with model assumptions, the extrapolating global lens reconstructions could separate the evidence from the assumptions. In this way, updating existing lens models with new observations becomes very efficient. In contrast, if some of the multiple images are not spectroscopically confirmed, the global non-extrapolating lens reconstruction has to be redone.

In Section 3, we provided proof-of-principle example applications. They showed that our method determines local lens properties for galaxy-scale and galaxy-cluster-scale lenses and in optical and radio bands alike. For the galaxy-cluster-scale example, CL0024, we found that our approach obtained local lens properties in agreement with the local lens properties obtained by two different global lens reconstructions. For the galaxy-scale example, B0128, this comparison is still outstanding at the time of writing, but we already found that our approach and a third global lens reconstruction algorithm consistently hinted at an oversimplified parametric lens model that was assumed in previous works.

On the galaxy-cluster scale, we found that the model-independent approach is also of great use for the reconstruction of the source from the observables of the respective multiple images. Our approach first determines the local lens properties in the vicinity of the multiple images and then infers the source by a local back-projection of the multiple images. Hence, the reconstructed source does not depend on a lens model. This enables us to study the morphologies of faint and distant galaxies to high redshifts free of potential biases due to a lens model.

For galaxy-scale lenses that show giant arcs, we could establish a simple set of equations that clearly show the degeneracies between an axisymmetric main lens and a smaller perturbing mass density. These equations explain the degeneracies between the main lens and the perturber which were previously discovered in simulations. Furthermore, they show that most of the observed examples of perturbed, symmetric giant-arc image configurations do not allow us to simultaneously constrain the mass-density distribution in the main lens and properties of the perturber, like its position along the line of sight or its total mass. With Bayesian lens-modelling algorithms, the most likely composition of the deflecting mass-density distribution can be determined. Prior to an elaborate modelling and Monte Carlo simulations, our set of model-independent equations is an efficient means to exclude classes of models and parameter sets that are incompatible with the observations.

Strong gravitational lenses are also used to infer the Hubble-Lemaître constant, $H_0$. As the ratio of angular diameter distances in the time-delay equation is dependent on the cosmological background model, calculating $H_0$ from this equation implies to either fix the other parameters of the cosmological model by complementary observations, or marginalise over them. Both ways are difficult to realise in particular due to the unknown form of the parameter attributed to a cosmological constant or dark energy. To avoid this problem, we set up a ratio of angular diameter distances that is based on the standardisable supernovae type Ia from the recent Pantheon sample. The ratio of angular diameter distances is thus determined up to an overall scale factor, $H_0$, and all distances (up to this overall scale) are expressed in terms of a basis-function decomposition. Purely data-based distance measures cannot be established due to the sparsity of the supernovae. However, the flexible set of orthonormal, analytic basis functions captures all characteristics of the data set. In addition, it avoids overfitting as the number of basis functions is limited by the measurement uncertainties of the supernovae. Using these data-based angular diameter distances in the lensing time-delay equation, the imprecision on $H_0$ only grows about a factor of two for the current measurement uncertainty and imprecision in the lens reconstruction compared to using angular diameter distances based on the currently best-fit cosmological standard model. This imprecision can be compensated for by increasing the number of jointly evaluated lensing configurations. If multiply-imaged fast radio bursts are observed in the future, they can replace the commonly used quasar sources, which will also yield an increase in the precision of $H_0$ inferred from time-delay difference measurements.

On the whole, we have introduced a new approach to strong gravitational lensing with the aim to separate model-based assumptions from data-based evidence which is a complementary and useful tool to analyse large data sets in the most efficient and objective way. A detailed discussion of the results achievable with current observations can be found in Section 3. We plan to further extend the approach including higher-order features of the multiple images and additional non-lensing observables.

Figure 8 summarises our approach (grey boxes) with its prerequisites (yellow boxes) as outlined in Section 2. It shows the underlying working principle and the resulting local lens properties that can be obtained from the single-lens-plane standard gravitational-lensing formalism without using a specific lens model. The remaining degeneracies and the way to break them are also listed in addition to those cases to which the algorithm cannot be applied.