© ESO 2018

1 Introduction

High-redshift galaxies can be magnified by strong gravitational lensing such that their multiple images show resolved features, as has been observed by e.g. Colley et al. (1996), Donnarumma et al. (2011), Sharon et al. (2012). In this way, the properties and evolution of these faint and otherwise barely observable galaxies can be studied. When resolved features within several of the multiple images can be identified, Tessore (2017) established a general expression to link the transformations between any configuration of more than two multiple images with properties of the lens mapping. Information about the ratios of convergences and the reduced shear at the position of the images can then be obtained without assuming any lens model. In this work, we investigate which model-independent information about the gravitational lens can be retrieved from the transformation between multiple images with resolved features in terms of a Taylor-expanded lens potential, as first introduced in Wagner & Bartelmann (2016), in fold and cusp configurations close to the critical curve.

In Sect. 2 we revise the approach based on the Taylor-expanded lensing potential (Wagner & Bartelmann 2016; Wagner 2017), and the approach established in Tessore (2017) in a unified notation. We show that both approaches are equivalent for fold and cusp configurations close to the critical curve. Subsequently we compare and combine them to extend the model-independent knowledge retrievable from multiple image configurations. After the theoretical derivations, we discuss the accuracy and precision achievable, and in Sect. 3 briefly introduce the algorithmic implementation to obtain the lens parameters, i.e. the ratios of convergences, reduced shear components, or ratios of potential derivatives, from the linear transformation between resolved features in multiple images. This is detailed further in Wagner et al. (2018).

Section 4 then shows an example application: a simulated set of four multiple images in a singular isothermal elliptical (SIE) lens model. By means of this simulation we analyse how accurately the lens parameters can be determined for the different choices of variables of Tessore (2017) and Wagner (2017) in more detail and compare the results with the results that are gained when using the quadrupole moment of the individual images as observable instead of the transformation matrix. We also investigate the influence of the size of the multiple images and the influence of detection noise on the accuracy of the lens parameter reconstruction. Section 5 summarises the results and gives an outlook to the observational cases that can be analysed with this approach.

2 Transformation matrix for a Taylor-expanded lensing potential

2.1 Definitions and notations

Let ϕ(x, y) be the gravitational lensing potential that defines the lens mapping between in the image plane and in the source plane by ∇_xϕ(x, y) = 0. The critical curves are all points x₀ for which the lens mapping becomes singular. Mapping them into the source plane, we obtain the caustic points y ₀. As already introduced in Schneider et al. (1992), Wagner & Bartelmann (2016), and Wagner (2017), the most convenient coordinate system for characterising the lens mapping in the vicinity of the critical curve using a Taylor-expansion of ϕ around a point x₀ on the critical curve is given by the conditions (1)

A subscript i of ϕ denotes the partial derivative in the direction of x_i, i = 1, 2; the superscript (0) indicates that this variable is evaluated at the singular point. Analogously, superscripts of capital letters ((A), (B), etc.) denote that the variable is considered at the centre of light of the respective image. Furthermore, we define the relative distance between the centres of light of image I and J, x _I and x _J, as .

Given two images A and B and their distortion matrices, often also called magnification matrices M_A(x) and M_B(x), we define the transformation matrix, also called relative magnification matrix T, which maps image A onto image B by (2)

as also stated in Gorenstein et al. (1984), and Narayan & Bartelmann (1996). The transformation matrix between the two images can be interpreted as a linear approximation of the map that first projects points of image A, x _{A} into the source plane by y (x _{A}) ≈ M_A(x_{A}) and subsequently maps the source onto the points in image B by , as Fig. 1 shows.

In the coordinate system given by Eq. (1) and using the Taylor-expanded lensing potential, the distortion matrices are (3)

with the matrix entries

in the vicinity of a fold or a cusp critical point x₀. The fold and cusp singular points are defined by their back-projections into the source plane, i.e. the respective points on the caustics: back-projecting a fold singular point, the respective point on the caustic has a non-vanishing tangential vector, while the tangent to the cusp singular point is mapped to zero when back-traced into the source plane. For the cusp point to be a stable singularity, the normal vector must not vanish. Practically, we locate folds and cusps by the multiple images we observe in their vicinity: a fold lies on the connection line between the centres of light of two multiple images that are mirror images of each other (see Fig. 1 (left)). A cusp is the vertex point of the parabola that lies between three multiple images as shown in Fig. 1 (right).

To be consistent with the notation given above, we reformulate the transformation equations, Eqs. (17)–(20) in Tessore (2017), for n multiple images of a gravitationally lensed source as

where the subscripts of T i, j = 1, 2 denote the entries of the transformation matrices between the reference image, called image A without loss of generality, and the remaining n − 1 multiple images I. As usual, κ denotes the convergence and g = γ∕(1 − κ) is the reduced shear, for which a subscript i = 1, 2 denotes its components, as defined in Bartelmann & Schneider (2001).

Fig. 1 Transformation between image A and B at a fold singular point (left) assembled from projecting A back to the source using MA and projecting the source onto B using (left, centre). Transformations between three images at a cusp singular point, projecting the central image A onto B using and projecting Aonto C using (right). The transformation between B and C, e.g. , is not shown.

2.2 Equivalence of both approaches

The equivalent set of equations using the derivatives of the lensing potential as variables reads

with the abbreviation (16)

Given at least three multiple images of one source with their transformation matrices and introducing the variable (17)

as in Tessore (2017), both systems of equations can be solved, one for , the other for . Both solutions are transformable into one another and are unique.

If all observable images are aligned, both systems of equations are underdetermined. This is the case for the three images of a source that is mapped by a Navarro-Frenk-White (NFW) profile, for example. In general, every axisymmetric lensing configuration, regardless of the specific density profile, generates multiple images that are positioned on a line through the symmetry centre of the lens, which means that the method cannot be applied to these cases. To demonstrate this degeneracy, consider the transformations between the three multiple images generated by an NFW profile. Their quadrupole moments are aligned with each other and extend tangentially or orthogonally to the critical curve, implying that their transformation matrices are diagnoal. Hence, the systems of equations, Eqs. (8)–(11) and Eqs. (12)–(15) for the three images A, B, C, denoting the reference image as A, reduce to four equations (two equations of T₁₁ from the combination (A, B) and (A, C) and two equations of T₂₂ from the same combinations) to solve for five variables, or . Hence, the system is underdetermined. Analogously, the transformation matrix also cannot be employed for two images that straddle a fold and are oriented orthogonally to the critical curve.

In the approach of Tessore (2017), no special coordinate system is required, while the system established in Eqs. (12)–(15) requires a coordinate system in which the images are extended along the x₂ -axis. Furthermore, the direction of the semi-major axis of the image quadrupole moments remains unchanged by the coordinate transformation, i.e. the sign of the off-diagonal entries in the magnification matrix is kept. Appendix A shows the solution to the system of equations in Eqs. (12)–(15) in these coordinates for three images, in analogy to the system derived in Tessore (2017). Yet, as detailed in Sect. 3, we pursue another way to solve for the lens parameters here. Measuring the quadrupole moment of an image that need not have resolved features, the reduced shear and second-order potential derivatives at its centre of light are determined by its axis ratio and orientation angle, as shown in Wagner (2017).

2.3 Transformations close to the critical curve

For two images A and B close to a point x₀ on the critical curve, in the special coordinate system of Eq. (1), (18)

holds, which eliminates one of the five variables in Eqs. (8)–(11) or Eqs. (12)–(15), so that the system is exactly solved using the transformation matrix between the two images. If additional images are present, the validity of the limit in Eq. (18) can be tested by calculating f_ϕ from Eqs. (12) to (15). Small deviations from the limit imply that the two images are so close to the critical curve that the intrinsic source ellipticity and orientation to the caustic are negligible and observable image ellipticities lead to the results deduced in Wagner (2017). Analogously, small deviations from the limit imply that the following ratios of second-order derivatives at the image positions can be determined when using the transformation matrix between the two images, as observable instead of their quadrupole moments

where we denote the image with positive parity as A and assume T₂₁≠0. For better readability, we drop the superscripts of T when considering only one pair of images. Inserting these relations into the lensing equations, we obtain

for the two images close to a fold.

For three images at a cusp A, B, C, with A being the reference image closest to the cusp, we obtain

and analogous relations using images A and C. Equation (18) has to be valid for all three images. Determining the , I = B, C, shows the goodness of approximation of Eq. (18).

Details about the derivations of Eqs. (23)–(26) are given in Wagner (2017), as the calculations employing the quadrupole moment are equivalent until the ratios of second-order derivatives at the image positions are replaced by the observables. Hence, we obtain the same ratios of potential derivatives as results which determine the shape of the critical curve in the vicinity of the images, as derived in Wagner (2017).

In addition, we simplify the reconstruction of the relative image position of A established in Wagner (2017) using

which is directly derived from Eqs. (5) and (7). Having determined x _A, we know the position of the cusp critical point.

If T₂₁ = 0, , which implies that both are zero, considering Eq. (5). From this it follows that T₁₁ = 1, T₁₂ = 0, and . As a consequence, we infer that the images are oriented orthogonally to the critical curve and that the third equation is underdetermined, so that we cannot retrieve .

Hence, if a transformation with T₂₁ = 0 occurs that does not also have T₁₁ = 1 and T₁₂ = 0, the considered image configuration is inconsistent with our approach, which may hint at local asymmetries in the lensing potential, microlensing, dust extinction, or may suggest that the two images do not originate from the same source, as already detailed in Wagner (2017). Given more than two images, i.e. an overdetermined system of equations, we can also use the surplus constraints to detect such anomalies.

2.4 Relative parity information

As can be easily derived from Eq. (2), (29)

holds for the magnification ratio μ_B∕μ_A. Hence, the sign of the determinant of the transformation yields the relative parities between the images. Given the parity of one image, all parities can be fixed. For instance, assuming that the faint image of a source very close to the lens centre is a maximum in the lens mapping and thus of positive parity, it can be used to determine the parities of the remaining images. The five multiple images of a source at z_s = 1.675 in the galaxy cluster Cl0024+1654, Colley et al. (1996) are an example for which the image parities can be determined in this way (see Fig. 2).

Deviations between the absolute values of det(T) and the observed magnification ratios may hint at microlensing or dust extinction.

Fig. 2 Cl0024+1654 with details of a five-image system. The arrows indicate images of opposite parity, i.e. which are mirror images of the other. Given that the central image 5 is of positive parity, the parities of images 1–4 can be fixed. (Image credits: NASA, ESA, M. J. Jee (Johns Hopkins University)).

2.5 Accuracy and precision

The results of Tessore (2017) are derived for the idealised case of a transformation between images at infinite resolution in the absence of noise. In realistic cases, assumptions have to be made to deal with finite image resolutions and noise. One possibility is to assume that the images are small compared to the scale on which the lens properties change, so that the entries of the transformation matrix are constant (Wagner et al. 2018).

As two examples to test the range of validity of this approximation and to show differences between the choice of variables, we chose NFW profiles and SIE lenses. We first simulate NFW profiles (30)

in spherical coordinates with radius r; scale density ρ_s; and scale radius r_s with masses M = 10¹⁵ M_⊙ and M = 0.5 × 10¹⁵ M_⊙ at redshift 0.5; and a concentration c_NFW ≡ r₂₀₀∕r_s = 3 (see Merten et al. 2015 for examples of cluster reconstructions in this parameter range), where r₂₀₀ denotes the radius at which the density of the NFW profile is 200 times the critical density of the universe at redshift 0.5. The source planeis assumed to be at redshift 2.0.

Next, we simulate a SIE lens given by the convergence for each point x = (r, φ) in polar coordinates in the image plane (31)

with scale radius r_s for the axis ratios f = 0.9 and f = 0.8 in the notation of Kormann et al. (1994). To fix r_s, we assume that the lens plane is at redshift 0.5, the source plane at 2.0, and the velocity dispersion is 1000 km s⁻¹, which is typical for galaxy cluster lenses.

Sources are placed at increasing distance from the fold and cusp points in the source plane and for the cusp, the sources are placed on the symmetry axis connecting the lens centre with the cusp.

For both profiles, we investigate the goodness of approximation of Eq. (18) dependent on the distance to the critical curve when using convergence and reduced shear and when using potential derivatives. As the distance to the critical curve is not observable, we plot the ratios of convergences and potential derivatives with respect to the observable relative image distance. Figure 3 shows the results using κ, as in Tessore (2017), and Fig. 4 the results using ϕ₁₁, as in Wagner & Bartelmann (2016); Wagner (2017). For both profiles, image A is the one with positive parity, image B the one with negative parity. In the cusp case, image B is the one closest to the cusp.

From a comparison of Fig. 3 with 4, we deduce that both choices of parametrisation are equally well suited for fold configurations of SIE lenses, the ansatz using the convergence is superior to the one with potential derivatives for cusp configurations of SIE lenses and vice versa for fold configurations of NFW profiles. The very small decreasing slope for the NFW profile (right plot in Fig. 4) explains the high accuracy of the Taylor approximation that we found in Wagner (2017). In that case, the quadrupole moments have to be used as observables because the system of equations using the transformation matrix is underdetermined.

We notethat ratios of ϕ₁₁ smaller than one for a SIE cusp configuration are achieved when interchanging A and B with respect to the choice of Fig. 3, as , while at the folds of both profiles holds.

The precision of the ratios of potential derivatives is given by the precision of the transformation and the relative distances between thecentres of light of the images. As an example, Gorenstein et al. (1984) obtained the eigenvalues of T with uncertainties below 5% and the direction of the semi-major axis of T with an uncertainty of 7% in a χ² -parameter estimation mapping image A in Q0957+561 onto image B in VLBI observations. Compared to these uncertainties for the T_ij, the relative uncertainty of the distance between A and B of the order of is negligible.

Fig. 3 Accuracy of the approximation stated in Eq. (18) for the variables used in Tessore (2017) close to asymmetric cusp configuration of a SIE lens (left), close to a fold configuration of a SIE lens (centre), and close to the fold configuration of an NFW profile (right). Distances between the images |δ AB | are determined as observed on the sky, i.e.rescaled to rs∕Dd with the angular diameter distance from the observer to the lens Dd. Black solid lines denote the accuracy for a SIE lens with f = 0.9 or an NFW with M = 1.0 × 1015M⊙, and red dashed lines for a SIE lens with f = 0.8 or an NFW with M = 0.5 × 1015 M⊙. (Notation for the models as defined in Kormann et al. 1994 and Bartelmann 1996.)

Fig. 4 Accuracy of the approximation stated in Eq. (18) for the variables used in Wagner (2017) close to a symmetric cusp configuration of a SIE lens (left; here A and B are interchanged), close to a fold configuration of a SIE lens (centre), and close to thefold configuration of an NFW profile (right). Distances between the images |δ AB | are determined as observed on the sky, i.e.rescaled to rs∕Dd with the angular diameter distance from the observer to the lens Dd. Black solid lines denote the accuracy for a SIE lens with f = 0.9 or an NFW with M = 1.0 × 1015M⊙, and red dashed lines for a SIE lens with f = 0.8 or an NFW with M = 0.5 × 1015 M⊙. (Notation for the models as defined in Kormann et al. 1994 and Bartelmann 1996.)

3 Implementation

In this work, we identify resolved features in the intensity distributions of the images manually, and leave automated and more advanced techniques to further work.

The transformation matrix has four entries and thus four degrees of freedom that amount to a rotation angle, two stretching factors, and a parity inversion as physical degrees of freedom. Hence, we need at least three positions, i.e. two linearly independentvectors, in each image that can be mapped to their counterparts in another image to determine the transformation matrix.

In the vicinity of a fold in the leading order approximation we consider two points, i.e. one vector in each image, which is sufficient because there is always a parity inversion that is fixed and the only degrees of freedom are the rotation angle and one stretching factor to change the area and thus determine the absolute value of the magnification ratio caused by gravitational lensing.

If more than the necessary number of points are given, the system of equations is overconstrained. If the necessary number of images is given, the most efficient way to solve for the lens parameters is to parametrise the entries of the transformation matrices by using theright-hand sides of Eqs. (8)–(11) or Eqs. (12)–(15) and inserting them into the system of transformation equations (32)

where denote the reference points of the individual images to which the reference points in the reference image are mapped, x _I are their centres of light, and x_A the centre of light of the reference image. This system of equations is solved for the lens parameters in T^(A,I) as a non-linear least-squares parameter estimation problem (NLLSP), minimising the deviations between the right-hand side and the left-hand side. To account for different possibly correlated uncertainties of the reference point positions in each image, the deviations are weighted by the inverse variance of these uncertainties.

To determine the higher-order ratios of potential derivatives at the critical points, we further replace the second-order potential derivatives by Eqs. (4)–(7). Subsequently, the resulting lens parameters and their covariances are taken as input for an importance sampling Monte Carlo simulation. The median lens parameter values of all samples and their 1-σ confidence intervals (simply denoted as confidence intervals) are the final result. The number of samples is chosen to be at least 10 000 to assure that the weights of the importance sampling are well balanced.

If more than the necessary number of images is given, we do not parametrise T by means of the lens parameters, but employ a different more efficient parametrisation for T in the system of Eq. (32) that reduces the non-linearity of the NLLSP. This implementation detail is further described in Wagner et al. (2018) and is used to determine the lens parameters in Sect. 4.3.

4 Example application

4.1 Simulation of sources

In order to set up the transformation matrices between the images by the reference matching described in Sect. 3 and to estimate their uncertainties, we simulate a source galaxy as a superposition of four round Sérsic profiles with Sérsic index n = 4 to obtain four reference points in each image. The first source consists of a Sérsic profile with r_e = 15 px, denoted as S₀, and three adjacent Sérsic profiles with r_e = 5 px, denoted as S₁, S₂, S₃, at offsets (2 px, 0 px), (4 px, 4 px), and (0 px, 4 px) from the centre of S₀. For an easier identification of the points, we scale the intensity distributions of the three smaller Sérsics to 0.1, 0.3, and 0.7 of the large one. Figure 5 (left) shows the logarithmic intensity profile for this source.

The second source has the same Sérsic profiles as the first, but the three smaller Sérsic profiles are at offsets (3 px, 0 px), (5 px, 5 px), and (0 px, 5 px). Hence, the reference points in the second source are farther away from each other, so that the accuracy of the local lens properties dependent on the distance between the reference points within the images can be investigated.

Fig. 5 Logarithmic intensity distribution of a source consisting of four Sérsic profiles, described in Sect. 4.1, simulating a galaxy with resolved brightness features (left). Multiple images of this source generated by the SIE lens introduced in Wagner (2017) and described in Sect. 4.2 (right). The lens centre is indicated by a red cross; the small red circle marks the point through which the horizontal and vertical intensity profiles run (shown beneath andto the right of the picture).

4.2 Simulation of the lens

We use the SIE lens introduced in Wagner (2017) as it is the simplest lens for which all the discussed image configurations can be analysed. Figure 5 (right) shows the images generated by the SIE lens for the first source described in Sect. 4.1. Placing the second source at the same position in the source plane, the semi-major and semi-minor axes of the quadrupole moments of its multiple images are about 30% longer than the axes of the images of the first source. The image configuration (A, B, C) is an asymmetric cusp configuration in which the central image A has positive parity. In Wagner (2017) we analysed symmetric cusp configurations with a central image of negative parity in detail.

The extensions of these simulated multiple images match the scale of observed ones and, as already stated in Wagner (2017), the semi-major axis of the critical curve of the SIE lens has an extension of 23′′ and is thus similar to observed effective Einstein radii, as detailed in e.g. Zitrin et al. (2015), and the effective Einstein radius of the simulated Hera cluster, as detailed in Meneghetti et al. (2017).

4.3 Second-order lens parameters of all images

We first determine the ratios of convergences and reduced shear components from the observables of all multiple images fromthe first source by employing the algorithm outlined in Sect. 3. All weights in the NLLSP optimisation are set to 1; 10 000 samples are used for the Monte Carlo simulation; and image A, the positive parity image closest to the cusp, is chosen as reference image.

Figure 6 shows the resulting probability density distributions of the f_κ - and g_i -values, i = 1, 2, for all four multiple images at their centres of light. The confidence intervals for and are smaller than those for the remaining lens parameters because all remaining f_κ and g_i depend on and .

To compare the reconstructed f_κ and g_i with their true values at the centre of light, we show the true value in each probability density distribution as a red line. We observe that all true values are within the confidence intervals, except for .

4.4 Second-order lens parameters at a fold

Considering the image pair (A, B) in Fig. 5 as a fold configuration, we determine the reduced shear components under the assumption that f_κ ≈ 1 and for the same implementation specifications as described in Sect. 4.3. As can be read off the first two columns in Fig. 7, the true g₁ at the centre of light lie within the confidence intervals of both images, while the g₂ lie outside.

Subsequently, we calculate probability density distributions of the ratios of second-order potential derivatives and plot them in the last two columns in Fig. 7. This parametrisation leads to worse results compared to that using the reduced shear components because the true ratios of second-order potential derivatives for image A are farther away from the most likely ratios of the second-order potential derivatives obtained by our approach. The value of is not even contained in the 3-σ confidence interval that limits the plot range. Calculating and , we find that both values have a large deviation from Eq. (18), which explains why the reconstructed lens parameters show high inaccuracies for both parametrisations. In addition, we conclude that for increasing ellipticity of the SIE lens, the parametrisation using the convergence ratios and reduced shears yields more accurate results than that using the ratios of the potential derivatives.

Fig. 6 Comparison of the probability density distributions of the convergence ratios fκ and reduced shear components gi, i = 1, 2, of the SIE lens at the centres of light of the four multiple images shown in Fig. 5 with their true values (red lines). The dark grey shaded area, the grey shaded area, and the light grey shaded area delimit the 1-, 2-, and 3-σ confidence intervals, respectively.

Fig. 7 Comparison of the true lens parameters (red lines) at a fold with the probability density distributions of the lens parameters determined by our approach in the parametrisation of Tessore (2017) (first two columns) and in the parametrisation of Wagner (2017) (last two columns).

4.5 Third-order lens parameters at a fold point

Keeping the implementation specifications as described in the previous sections, we determine the probability density distributions for the third-order ratios of potential derivatives as shown in Fig. 8. As expected from the results of Sect. 4.4, the true values of the ratios of potential derivatives at the fold point lie outside the confidence intervals – the true value of lies beyond the plot range – and the reconstruction of the critical curve in the vicinity of the fold will consequently be very inaccurate.

Fig. 8 Comparison of the true third-order ratios of potential derivatives (red lines) at a fold with the probability density distributions of the ratios of potential derivatives determined by our approach.

4.6 Higher-order lens parameters at a cusp point

The accuracy of the reconstruction of the critical curve in the vicinity of the cusp is summarised in Fig. 9, employing the same implementation specifications as before. The first two columns are obtained using the image pair (A, B), the last two columns are obtained using the image pair (A, C) to determine the ratios of potential derivatives and the centre of light position of image A. From the first row, we conclude that the accuracy of the slope of the parabola that approximates the critical curve in the vicinity of a cusp is very low because the true values of the third- and fourth-order ratios of potential derivatives lie far off the 3-σ confidence intervals. Contrary to that, the true value of x_A1, the centre of light position of image A, lies within the confidence interval. Using image pair (A, B), we find that x_A2 lies within the confidence interval as well, so that the true cusp position is determined in more than 68% of the cases.

Fig. 9 Comparison of the true ratios of potential derivatives and image positions (red lines) at a cusp with the probability density distributions of the ratios of potential derivatives and image positions determined by our approach. The first two columns use the image pair (A, B) and the last two columns (A, C) for the reconstruction.

4.7 Comparison to the quadrupole moment as observable

To compare the results obtained by means of the transformation matrix with the results retrievable from the quadrupole moment as observable, we have to specify uncertainties for the quadrupole moment observables, the axis ratio r_I, and the orientation angle of the semi-major axis φ_I, I = A, B, C in order to generate the probability density distributions for the lens parameters. As a first estimate, we assume that all uncertainties in the observables follow a Gaussian distribution and are uncorrelated. The distance between the images is usually determined to high precision (see Sect. 2.5), so that a standard deviation of 10⁻⁴ δ_IJi, i = 1, 2, for the coordinate-wise distances between the images seems appropriate. For the axis ratio, we estimate a standard deviation of 0.1 r_I and for the orientation angle ± 5° from visual inspection of some HST observations. From these normal distributions of observables that are obtained for the first source in our SIE simulation, we generate 10 000 samples and determine the probability density distributions for the ratios of potential derivatives.

The results for the second-order ratios of potential derivatives at a fold are plotted in the first two colums of Fig. 10¹. They show that employing the quadrupole yields slightly higher probabilities of retrieving the true values of the second-order potential derivatives, compared to using the transformation matrix (see the last two columns in Fig. 7). However, the ratios of third-order potential derivatives, shown in the last two columns of Fig. 10, have larger confidence intervals and a very low probability of retrieving the true values. Therefore, using the quadrupole as observable yields worse results than employing the transformation matrix for the reconstruction of the critical curve at the fold.

Analogously, we determine the ratios of potential derivatives at the cusp point and the position of image A, once from the image pair (A, B) and once from the image pair (A, C), using the quadrupole moments of the images as observables. The results are plotted in the first two and the last two columns of Fig. 11, respectively. As for the transformation matrix, the image pair (A, B) has the highest probability of retrieving the true values for the ratios of potential derivatives and x _A. Comparing Fig. 11 with Fig. 9, using the quadrupole moment clearly yields higher probabilities for the true values than employing the transformation matrix, yet this comes at the cost of broader confidence intervals. Another disadvantage of the quadrupole moment as observable is that it may still contain a bias due to the intrinsic ellipticity of the source.

Fig. 10 Comparison of the true ratios of potential derivatives (red lines) at a fold with the probability density distributions of the ratios of potential derivatives determined by our approach using the quadrupole moment as observable instead of the transformation matrix. The first row of the last two columns uses image A and the bottom row uses image B for the reconstruction of the third-order potential derivatives.

Fig. 11 Comparison of the true ratios of potential derivatives and image position A (red lines) at a cusp with the probability density distributions of the ratios of potential derivatives and image position A determined by our approach using the quadrupole moment as observable instead of the transformation matrix. The first two columns use the image pair (A, B) and the last two colums use the image pair (A, C) for the reconstruction.

4.8 Influence of the image size

To investigate the influence of the size of the multiple images on the accuracy of the reconstruction, we repeat the analysis of Sect. 4.3 for the second source. Comparing the plots in Fig. 6 with those in Fig. 12, we observe that the confidence intervals in the former are larger, so that the true parameter values have a higher probability of being correctly determined for the smaller source. The reason for the decreasing accuracy with increasing image extensions is the decreasing validity of the assumption that the entries of the magnification matrix are constant over the extension of each multiple image. The semi-major axis of the quadrupole moment of image A is 9% of the distance between image A and B for the first source and 12% for the second, which gives an estimate of the scales on which higher-order moments should be taken into account.

Fig. 12 Comparison of the probability density distributions of the convergence ratios fκ and reduced shear components gi, i = 1, 2 of the SIE lens at the centres of light of the four multiple images of the second source with their true values (red lines). The dark grey shaded area, the grey shaded area, and the light grey shaded area delimit the 1-, 2-, and 3-σ confidence intervals, respectively.

4.9 Influence of noise

Taking into account the noise in the observations, the reference points may not be accurately determined at the pixel-level, even when setting them manually. Keeping all other implementation specifications as before, we increase the uncertainty in all reference point positions from 1 px to 2 px, while assuming that they are still uncorrelated. In addition, we have to increase the number of samples to 200 000, otherwise the number of effective samples in the importance sampling step falls below 100 samples. This leads to significant multi-modal artefacts in the probability density distributions.

Comparing the resulting lens parameter probability density distributions shown in Fig. 13 with those in Fig. 6, we observe that the confidence intervals are increased such that the true lens parameter values are more often correctly retrieved. Giving different uncertainties to the reference point positions in different images, we observe that the confidence intervals of the respective lens parameters are broadened accordingly, yet the confidence intervals can be as large as the parameter value itself, especially for the f_κ, so that the approach becomes unreliable for increasing noise.

Fig. 13 Comparison of the probability density distributions of the convergence ratios fκ and reduced shear components gi, i = 1, 2 of the SIE lens at the centres of light of the four multiple images of the first source with their true values (red lines) with an increased uncertainty in the reference positions of 2 px. The dark grey shaded area, the grey shaded area, and the light grey shaded area delimit the 1-, 2-, and 3-σ confidence intervals, respectively.

5 Conclusion

We investigated what model-independent information about a gravitational lens can be gained by linearly mapping multiple images with clearly resolved features onto each other close to fold or cusp critical points, and how accurately these properties can be determined. Our results can be stated as follows:

• The approach developed in Tessore (2017) can be re-parametrised in terms of ratios of potential derivatives to show that it yields the same information about the magnification matrices in the vicinity of the critical curve as the approach developed in Wagner (2017).

• The approach of Tessore (2017) determines the magnification matrices up to a scale factor for all multiple images (i.e. convergence ratios and reduced shear components) at the image positions, but it cannot be applied to images generated by axisymmetric lenses.

• The approach of Wagner (2017) determines the magnification matrices up to a scale factor at the image positions and an approximation to the critical curve in the vicinity of the images close to a fold or a cusp critical point, but it cannot determine the lens properties at the position of the counter image lying on the opposite side of the lens centre or the position of the central image close to the lens centre.

• Both approaches in their current implementation assume that the convergence and shear variations are negligible over the extensions of the individual multiple images, so that images can be mapped onto each other using a linear transformation.

• Combining the two approaches yields model-independent information about the magnification matrices at the positions of all multiple images, their relative parities, and allows the critical curve to be reconstructed in the vicinity of fold and cusp critical points.

• We simulate a galaxy-cluster scale SIE lens and a source galaxy consisting of four elliptical Sérsic profiles to mimic a galaxy that shows clearly resolved features in its intensity distribution. In each image, we need to identify at least three distinctive features from which we then assume that they are mapped onto each other by the transformation matrix. As we analyse the five-image configuration of a source at redshift z_s = 1.675 in the galaxy cluster Cl0024+1654 as a first observational test (Wagner et al. 2018), the scales and extensions of our simulated lens and source are chosen to be similar to that observation. Hence, the following results from analysing the simulated multiple-image configurations can be used as a calibration for the observational case.

• A comparison of the parametrisations from Tessore (2017) and Wagner (2017) shows that the symmetric approach using ratios of convergences and reduced shear components is more robust than that using ratios of potential derivatives because the former does not require a transformation into a special coordinate system. Furthermore, it yields higher probabilities of retrieving the true magnification matrices up to a scale factor.

• A comparison between using the transformation matrix between pairs of images to retrieve the lens parameters and using the quadrupole moments (characterised by their axis ratios and orientation angles, as detailed in Wagner 2017) of the individual images reveals that the approach using the quadrupole moments yields higher probabilities of retrieving the true lens parameters for the cusp and yields lower probabilities of accurately retrieving the lens parameters for the fold.

• The approximation that convergence and shear are constant over the image areas becomes inaccurate when the semi-major axis of the quadrupole moment of the images is of the order of 10% of the relative distance between the centres of light of the images.

• The confidence intervals are in the range of 10–40% of the estimated lens parameter values and increase to 50–100% when the uncertainties in the positions of the resolved features exceed a precision of one pixel.

On the whole, we simulated multiple images with realistic extensions and the distance scales that have been observed. Taking into account that multiple images with clearly resolved features seem to be more extended than unresolved multiple images, we conclude that the linear transformation matrix between pairs of these images may not be sufficient to yield accurate and precise information and the approach has to be extended including higher-order moments in order to become a reliable tool for the model-independent characterisation of gravitational lenses.

Acknowledgements

J.W. would like to thank MauricioCarrasco, Sven Meyer, Robert Reischke, Björn Malte Schäfer, Sebastian Stapelberg, Rüdiger Vaas, and the anonymous referee for helpful discussions, and gratefully acknowledges the support by the Deutsche Forschungsgemeinschaft (DFG) WA3547/1-1. N.T. acknowledges support from the European Research Council in the form of a Consolidator Grant, number 681431.

Appendix A Solution to the system of Eqs. (eqnlinking 12) to (eqnlinking 15)

Given an observation of three images A, B, C, with A being the reference image, in the coordinate system described in Sect. 2.2, the solution of the system of Eqs. (12)–(15) is calculated as

To arrive at these equations, we use (A.6)

as done analogously for the ratios of convergences and reduced shear components in Tessore (2017).

References

All Figures

	Fig. 1 Transformation between image A and B at a fold singular point (left) assembled from projecting A back to the source using MA and projecting the source onto B using (left, centre). Transformations between three images at a cusp singular point, projecting the central image A onto B using and projecting Aonto C using (right). The transformation between B and C, e.g. , is not shown.
In the text

	Fig. 2 Cl0024+1654 with details of a five-image system. The arrows indicate images of opposite parity, i.e. which are mirror images of the other. Given that the central image 5 is of positive parity, the parities of images 1–4 can be fixed. (Image credits: NASA, ESA, M. J. Jee (Johns Hopkins University)).
In the text

Fig. 3 Accuracy of the approximation stated in Eq. (18) for the variables used in Tessore (2017) close to asymmetric cusp configuration of a SIE lens (left), close to a fold configuration of a SIE lens (centre), and close to the fold configuration of an NFW profile (right). Distances between the images |δ AB | are determined as observed on the sky, i.e.rescaled to rs∕Dd with the angular diameter distance from the observer to the lens Dd. Black solid lines denote the accuracy for a SIE lens with f = 0.9 or an NFW with M = 1.0 × 1015M⊙, and red dashed lines for a SIE lens with f = 0.8 or an NFW with M = 0.5 × 1015 M⊙. (Notation for the models as defined in Kormann et al. 1994 and Bartelmann 1996.)

In the text

Fig. 4 Accuracy of the approximation stated in Eq. (18) for the variables used in Wagner (2017) close to a symmetric cusp configuration of a SIE lens (left; here A and B are interchanged), close to a fold configuration of a SIE lens (centre), and close to thefold configuration of an NFW profile (right). Distances between the images |δ AB | are determined as observed on the sky, i.e.rescaled to rs∕Dd with the angular diameter distance from the observer to the lens Dd. Black solid lines denote the accuracy for a SIE lens with f = 0.9 or an NFW with M = 1.0 × 1015M⊙, and red dashed lines for a SIE lens with f = 0.8 or an NFW with M = 0.5 × 1015 M⊙. (Notation for the models as defined in Kormann et al. 1994 and Bartelmann 1996.)

In the text

Fig. 5 Logarithmic intensity distribution of a source consisting of four Sérsic profiles, described in Sect. 4.1, simulating a galaxy with resolved brightness features (left). Multiple images of this source generated by the SIE lens introduced in Wagner (2017) and described in Sect. 4.2 (right). The lens centre is indicated by a red cross; the small red circle marks the point through which the horizontal and vertical intensity profiles run (shown beneath andto the right of the picture).

In the text

	Fig. 6 Comparison of the probability density distributions of the convergence ratios fκ and reduced shear components gi, i = 1, 2, of the SIE lens at the centres of light of the four multiple images shown in Fig. 5 with their true values (red lines). The dark grey shaded area, the grey shaded area, and the light grey shaded area delimit the 1-, 2-, and 3-σ confidence intervals, respectively.
In the text

	Fig. 7 Comparison of the true lens parameters (red lines) at a fold with the probability density distributions of the lens parameters determined by our approach in the parametrisation of Tessore (2017) (first two columns) and in the parametrisation of Wagner (2017) (last two columns).
In the text

	Fig. 8 Comparison of the true third-order ratios of potential derivatives (red lines) at a fold with the probability density distributions of the ratios of potential derivatives determined by our approach.
In the text

	Fig. 9 Comparison of the true ratios of potential derivatives and image positions (red lines) at a cusp with the probability density distributions of the ratios of potential derivatives and image positions determined by our approach. The first two columns use the image pair (A, B) and the last two columns (A, C) for the reconstruction.
In the text

Fig. 10 Comparison of the true ratios of potential derivatives (red lines) at a fold with the probability density distributions of the ratios of potential derivatives determined by our approach using the quadrupole moment as observable instead of the transformation matrix. The first row of the last two columns uses image A and the bottom row uses image B for the reconstruction of the third-order potential derivatives.

In the text

Fig. 11 Comparison of the true ratios of potential derivatives and image position A (red lines) at a cusp with the probability density distributions of the ratios of potential derivatives and image position A determined by our approach using the quadrupole moment as observable instead of the transformation matrix. The first two columns use the image pair (A, B) and the last two colums use the image pair (A, C) for the reconstruction.

In the text

	Fig. 12 Comparison of the probability density distributions of the convergence ratios fκ and reduced shear components gi, i = 1, 2 of the SIE lens at the centres of light of the four multiple images of the second source with their true values (red lines). The dark grey shaded area, the grey shaded area, and the light grey shaded area delimit the 1-, 2-, and 3-σ confidence intervals, respectively.
In the text

Fig. 13 Comparison of the probability density distributions of the convergence ratios fκ and reduced shear components gi, i = 1, 2 of the SIE lens at the centres of light of the four multiple images of the first source with their true values (red lines) with an increased uncertainty in the reference positions of 2 px. The dark grey shaded area, the grey shaded area, and the light grey shaded area delimit the 1-, 2-, and 3-σ confidence intervals, respectively.

In the text