© ESO, 2017

1. Introduction

Strong gravitational lensing provides a highly valuable tool to obtain mass properties of galaxies and galaxy clusters (see, e.g., Bartelmann 2010; Kochanek 2006, and references therein). In particular, multiple image systems yield strong constraints on the mass distribution. The mass enclosed within the Einstein radius presents the most robust galaxy mass estimate currently available. Furthermore, the shape of the mass distribution (e.g., ellipticity, orientation) is well defined.

However, mass estimates for radii smaller or larger than the Einstein radius are less accurate. If only a finite set of individual lensed compact images is observed, too few observational constraints are available and certainly no unique radial mass profile can be found. The situation changes somewhat if extended source components are lensed where the constraints on the mass distribution are much more stringent. Nonetheless, even if we could find a mass model which reproduces all constraints perfectly, such a mass model would not be unique either. The reason for this degeneracy is known since 1985 (Falco et al. 1985) and is called the mass-sheet transformation (MST). If a given surface mass density κ(θ) reproduces all observational constraints, then the whole family of mass models, ${\mathit{\kappa }}_{\mathit{\lambda }}\mathrm{\left(}{\theta }\mathrm{\right)}\mathrm{=}\mathit{\lambda \kappa }\mathrm{\left(}{\theta }\mathrm{\right)}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\lambda }\mathrm{)}\mathit{,}$(1)will do the same. In particular, the MST leaves all observables invariant except the time delay¹. The transformation (1) modifies the slope of the density profile with a constant factor λ. This affects mass measurements outside the Einstein radius θ_E and determination of the Hubble constant H₀ directly.

Schneider & Sluse (2013, hereafter SS13) presented two mass profiles (namely, a Hernquist profile plus a modified Navarro, Frank and White profile, as well as a power-law mass profile) which showed almost the same imaging properties, although they are not exactly related through an MST. Following this unexpected result it became apparent that an even more general invariance transformation than the MST exists. The so-called source-position transformation (SPT) was finally introduced in Schneider & Sluse (2014, hereafter SS14).

For isolated individual images many ambiguities for the lens equation exist. Local transformations of the lensing mass distribution, which still reproduce the positional constraints from the lensed images, lead to an infinite number of mass models (see e.g., Saha & Williams 1997; Diego et al. 2005; Coe et al. 2008; Liesenborgs & De Rijcke 2012). The MST as given in Eq. (1) is a global transformation and equivalent to an isotropic uniform stretching of the source plane by a constant factor λ. The SPT is based on a more general (global) transformation of the source plane coordinates. Such transformations , where denotes the transformed source position, give rise to a new deflection law . The new deflection law will in general not be a gradient field and thus cannot be obtained from the deflection caused by a lens. However, if the curl component of is sufficiently small, then one may find a lensing mass distribution which yields a deflection law which is very close to , so close that it cannot be observationally distinguished from . In this paper we will explore this possibility, which of course depends on the SPT . In particular, if this deformation is “too strong”, then the resulting cannot be approximated with the deflection due to a lens – this will restrict the freedom in choosing transformations .

The outline of the paper is as follows. In Sect. 2 we will recapitulate the principle of the SPT. We characterize the deviation of the deflection law from a gradient field quantitatively in Sect. 3 by finding a gravitational potential $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\psi }\end{array}$ such that $\begin{array}{c}{}_{{˜}}\\ {\alpha }\end{array}\mathrm{=}\mathrm{\nabla }\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\psi }\end{array}$ is as close as possible to the SPT-transformed deflection law . To do so, we will start from a variational principle and show that the modified deflection potential $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\psi }\end{array}$ has to fulfill Neumann boundary conditions. Those can be solved using a Green’s function, and the solution will be given explicitly for a circular region. Furthermore, a numerical approach will be presented to find degenerate deflection laws and their corresponding mass profiles. By considering a specific deformation function and assuming a positional accuracy on lensed image positions typical of the Hubble Space Telescope (HST), we will present in Sect. 4 the implications of the “allowed” SPTs on current mass profile determinations, regarding the radial mass profile and the angular structure of the lens. Different diagnostics for the change of the mass profile by an SPT, and how it can be distinguished from an MST, will be explored in Sect. 5 in terms of the aperture mass. Finally, we will discuss our findings in Sect. 6.

2. The principle of the source position transformation

In the following we will describe the principle of the SPT and its properties. For a more detailed account the reader is referred to SS14. We use standard gravitational lensing notation throughout this paper (see, e.g., Schneider 2006).

In general, a surface mass density distribution κ(θ) gives rise to a deflection law α(θ), where θ is the angular position in the lens plane, i.e., the observer’s sky. The mass distribution or convergence κ is defined as the ratio of projected surface mass density to the critical surface mass density, where the latter depends only on the angular diameter distances of lens and source. If that mass distribution is sufficiently concentrated (i.e., typically κ(θ) ≳ 1 for some region in the lens plane) a source may have multiple images, depending on its position relative to the deflector on the sky. Then, the source located at the (unobservable) position β will have its images at locations described by the solutions θ_i = β + α(θ_i) of the lens equation. Since multiple images are from the same source, we can deduce the constraints on the deflection law α(θ) as ${{\theta }}_{\mathit{i}}\mathrm{-}{\alpha }\mathrm{\left(}{{\theta }}_{\mathit{i}}\mathrm{\right)}\mathrm{=}{{\theta }}_{\mathit{j}}\mathrm{-}{\alpha }\mathrm{\left(}{{\theta }}_{\mathit{j}}\mathrm{\right)}\mathit{,}$(2)or likewise for an alternative deflection law as ${{\theta }}_{\mathit{i}}\mathrm{-}\hat{{\alpha }}\mathrm{\left(}{{\theta }}_{\mathit{i}}\mathrm{\right)}\mathrm{=}{{\theta }}_{\mathit{j}}\mathrm{-}\hat{{\alpha }}\mathrm{\left(}{{\theta }}_{\mathit{j}}\mathrm{\right)}\mathit{,}$(3)for all i<j, such that α(θ) as well as yield exactly the same sets of multiple images. If such equivalent deflection laws exist, they will correspond to source positions β = θ−α(θ) or , respectively (see Fig. 1).

Fig. 1 Illustration of the source position transformation. A source at β causes multiple images θ in the lens plane under the deflection law α. The same multiple images are obtained from a source at , provided the deflection law is changed to , according to Eq. (4).

We can now consider a one-to-one mapping that connects the original source coordinates to the new ones. This allows us to define the transformed deflection law as $\hat{{\alpha }}\mathrm{\left(}{\theta }\mathrm{\right)}\mathrm{=}{\alpha }\mathrm{\left(}{\theta }\mathrm{\right)}\mathrm{+}{\beta }\mathrm{-}\hat{{\beta }}\mathrm{\left(}{\beta }\mathrm{\right)}\mathrm{=}{\theta }\mathrm{-}\hat{{\beta }}\mathrm{(}{\theta }\mathrm{-}{\alpha }\mathrm{(}{\theta }\mathrm{\left)}\mathrm{\right)}\mathit{,}$(4)where in the last step we inserted the original lens equation β = θ−α(θ).

Hence, any bijective (i.e., one-to-one) function leads to an SPT which leaves the condition (2) invariant. Moreover, as can be deduced from the Jacobian of the modified lens equation, the relative magnification matrices and the relative image shapes between image pairs of the same source remain unchanged. However, the Jacobian will not be symmetric in general, and therefore cannot be written as the gradient of a deflection potential (i.e., is not a curl-free field). This implies that no corresponding mass distribution exists that yields a deflection angle , in general. However, it was shown in SS14 that the asymmetric part of the Jacobian can be small in realistic cases; this will be explored more quantitatively in Sect. 3. In the special case that the lens is axisymmetric and the transformation corresponds to a radial stretching of the form ${\beta ̂}\mathrm{=}\mathit{f}\mathrm{\left(}\mathrm{\right|}{\beta }\mathrm{\left|}\mathrm{\right)}{\beta }\mathit{,}$(5)the SPT is an exact invariance transformation: in this case, the Jacobian is symmetric, and for every transformation (5) and its corresponding deflection law there exists a corresponding axi-symmetric mass distribution .

Provided the curl component of is small, then we expect that there exists a mass distribution $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\kappa }\end{array}$ whose corresponding deflection law $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\alpha }\end{array}$ will be very similar to , in the sense that their difference is smaller than the astrometric accuracy of current observations. In this case, the SPT will be, for all practical purposes, a global invariance transformation for lenses.

3. The transformed mass distribution

3.1. The general method

Since the deflection law $\hat{{\alpha }}$ (4) is not a gradient field, it does not correspond to a deflection field caused by a gravitational lens. However, if the curl component of $\hat{{\alpha }}$ is sufficiently small, one may be able to find a deflection potential $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\psi }\end{array}$ and a corresponding deflection law $\begin{array}{c}{}_{\mathrm{˜}}\\ {\alpha }\end{array}\mathrm{=}\mathrm{\nabla }\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\psi }\end{array}$ such that the difference between $\hat{{\alpha }}\mathrm{a}\mathit{nd}\begin{array}{c}{}_{\mathrm{˜}}\\ {\alpha }\end{array}$ is small, for example smaller than the astrometric accuracy of current observations. Since only the region of the lens plane where multiple images occur is constrained by lensing observations, the difference $\hat{{\alpha }}\mathrm{-}\begin{array}{c}{}_{\mathrm{˜}}\\ {\alpha }\end{array}$ needs to be small only in a finite region, which we denote as .

We thus consider the “action” $\mathit{S}\mathrm{=}{\mathrm{\int }}_{\mathrm{𝒰}}{\mathrm{d}}^{\mathrm{2}}\mathit{\theta }{\left|\mathrm{\nabla }\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\mathrm{-}\hat{{\alpha }}\right|}^{\mathrm{2}}\mathit{,}$(6)for which we want to find a minimum.

Using this particular variational principle is just one possibilty of finding $\begin{array}{c}{{˜}}_{}\\ {\alpha }\end{array}$. One could also apply the Helmholtz theorem and decompose into its irrotational (curl-free) and solenoidal (divergence-free) part. This would lead to similar but not identical results for $\begin{array}{c}{{˜}}_{}\\ {\alpha }\end{array}$, thus not changing the main conclusions of this paper². Another possible ansatz would be to find a gradient deflection angle such that its maximum deviation from would be minimized; however, the solution of this problem seems to be much more difficult to find than our variational principle.

Equation (6) can be minimized by considering small variations of $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\psi }\end{array}\mathrm{\to }\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\psi }\end{array}\mathrm{+}\mathit{\delta }\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\psi }\end{array}$, and finding the conditions for which the action is stationary for all variations $\mathit{\delta }\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\psi }\end{array}$. Up to linear terms in $\mathit{\delta }\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\psi }\end{array}$, we find $\begin{array}{ccc}\mathit{S}\mathrm{+}\mathit{\delta S}& \mathrm{=}& {\mathrm{\int }}_{\mathrm{𝒰}}{\mathrm{d}}^{\mathrm{2}}\mathit{\theta }{\left|\mathrm{\nabla }\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\mathrm{+}\mathrm{\nabla }\mathrm{\left(}\mathit{\delta }\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\mathrm{\right)}\mathrm{-}\hat{{\alpha }}\right|}^{\mathrm{2}}\\ & \mathrm{=}& \mathit{S}\mathrm{+}\mathrm{2}{\mathrm{\int }}_{\mathrm{𝒰}}{\mathrm{d}}^{\mathrm{2}}\mathit{\theta }\mathrm{\nabla }\mathrm{\left(}\mathit{\delta }\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\mathrm{\right)}\mathrm{(}\mathrm{\nabla }\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\mathrm{-}\hat{{\alpha }}\mathrm{)}\\ & \mathrm{=}& \mathit{S}\mathrm{+}\mathrm{2}{\mathrm{\int }}_{\mathit{\partial }\mathrm{𝒰}}\mathrm{d}\mathit{s\delta }\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\mathrm{(}\mathrm{\nabla }\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\mathrm{-}\hat{{\alpha }}\mathrm{)}\mathrm{·}{n}\\ & & \mathrm{-}\mathrm{2}{\mathrm{\int }}_{\mathrm{𝒰}}{\mathrm{d}}^{\mathrm{2}}\mathit{\theta \delta }\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\mathrm{(}{\mathrm{\nabla }}^{\mathrm{2}}\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\mathrm{-}\mathrm{\nabla }\mathrm{·}\hat{{\alpha }}\mathrm{)}\mathit{,}\end{array}$(7)where we made use of Gauß divergence theorem. The boundary curve of is denoted as , ds is the line element of the boundary curve, and n(s) the outward directed normal vector. Requiring δS = 0 leads to the Neumann problem ${\mathrm{\nabla }}^{\mathrm{2}}\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\mathrm{=}\mathrm{\nabla }\mathrm{·}\hat{{\alpha }}\mathrm{=}\mathrm{:}\mathrm{2}\mathit{\kappa ̂}\mathrm{;}\mathrm{and} \mathrm{\nabla }\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\mathrm{·}{n}\mathrm{=}\hat{{\alpha }}\mathrm{·}{n}\mathit{,}$(8)where the first equation is required for all points , and the second one for all points on the boundary . The solution $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\psi }\end{array}$ of Eq. (8) is specified only up to an additive constant, since a constant in the deflection potential does not affect the deflection angle.

In order to solve the system (8), we can either use numerical standard methods for such boundary problems, or we can obtain the solution by means of a Green’s function. Both methods will be explored in this section.

3.2. Solving the Neumann problem numerically

We defined the convergence of the transformed deflection law to be . The curl component of is reasonably small if the closest curl-free approximation to (which is $\begin{array}{c}{{˜}}_{}\\ {\alpha }\end{array}$) is smaller than a chosen astrometric accuracy ε_acc$\mathrm{|}\hspace{0.17em}{\alpha ̂}\mathrm{(}{\theta }\mathrm{)}\mathrm{-}\begin{array}{c}{˜}\\ {\alpha }\end{array}\mathrm{(}{\theta }\mathrm{)}\hspace{0.17em}\mathrm{|}\mathrm{=}\mathrm{\left|}\mathrm{\Delta }{\alpha }\mathrm{\right(}{\theta }\mathrm{\left)}\mathrm{\right|}\mathit{<}{\mathit{\epsilon }}_{\mathrm{acc}}$(9)for all . To solve the system (8) numerically, we set up a successive overrelaxation method (SOR; Press et al. 1996, their Sect. 19.5) on a square grid to calculate $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\psi }\end{array}$. An SOR is a converging iterative process based on the extrapolation of the Gauß-Seidel method, and it is a standard method to solve boundary value problems (see, e.g., Seitz & Schneider 2001). Using a second-order accurate finite differencing scheme, the deflection law $\begin{array}{c}{{˜}}_{}\\ {\alpha }\end{array}$ is then derived from the deflection potential $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\psi }\end{array}$.

The lens is located at the center of the grid, chosen to be also the origin of the coordinate system. The grid has a length of 4 θ_E to cover the relevant area in which multiple images occur, i.e., it covers the region within 2θ_E from the lens center.

The SOR involves the calculation of a weighted average between the previous iterate ${\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\psi }\end{array}}_{\mathit{i,k}}^{\mathrm{(}\mathit{m}\mathrm{-}\mathrm{1}\mathrm{)}}$ and the computed Gauß-Seidel iterate $\stackrel{}{{{\mathrm{\Psi }}_{\mathit{i,k}}^{\mathrm{\left(}\mathit{m}\mathrm{\right)}}}_{\mathrm{˜}}}$ successively for each component ${\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}}_{\mathit{i,k}}^{\mathrm{\left(}\mathit{m}\mathrm{\right)}}\mathrm{=}\mathit{\omega }\hspace{0.17em}\stackrel{}{{{\mathrm{\Psi }}_{\mathit{i,k}}^{\mathrm{\left(}\mathit{m}\mathrm{\right)}}}_{\mathrm{˜}}}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\omega }\mathrm{)}\hspace{0.17em}{\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}}_{\mathit{i,k}}^{\mathrm{(}\mathit{m}\mathrm{-}\mathrm{1}\mathrm{)}}\mathit{,}$(10)where ${\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\psi }\end{array}}_{\mathit{i,k}}^{\mathrm{\left(}\mathit{m}\mathrm{\right)}}$ is the value of $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\psi }\end{array}$ for the grid point (i,k) in iteration m, and ω is the extrapolation parameter. The parameter ω is chosen such that it accelerates the rate of convergence of the iterative variable to the solution; in this work $\mathit{\omega }\mathrm{=}\frac{\mathrm{2}}{\mathrm{1}\mathrm{+}\mathit{\pi }\mathit{/}\mathrm{(}\mathrm{𝒩}\mathrm{-}\mathrm{1}\mathrm{)}}\mathit{,}$(11)is applied, where is the total number of grid points. Initially, all ${\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\psi }\end{array}}_{\mathit{i,k}}$ are set to zero. In each iteration m, the Gauß-Seidel iterate $\stackrel{}{{{\mathrm{\Psi }}_{\mathit{i,k}}^{\mathrm{\left(}\mathit{m}\mathrm{\right)}}}_{\mathrm{˜}}}$ is calculated as follows (a fourth-order accurate finite differencing is used) $\begin{array}{ccc}\stackrel{}{{{\mathrm{\Psi }}_{\mathit{i,k}}^{\mathrm{(}\mathit{m}\mathrm{+}\mathrm{1}\mathrm{)}}}_{\mathrm{˜}}}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathrm{60}}\hspace{0.17em}({\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}}_{\mathit{i}\mathrm{+}\mathrm{2}\mathit{,k}}^{\mathrm{\left(}\mathit{m}\mathrm{\right)}}\mathrm{+}{\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}}_{\mathit{i}\mathrm{-}\mathrm{2}\mathit{,k}}^{\mathrm{\left(}\mathit{m}\mathrm{\right)}}\mathrm{+}{\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}}_{\mathit{i,k}\mathrm{+}\mathrm{2}}^{\mathrm{\left(}\mathit{m}\mathrm{\right)}}\mathrm{+}{\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}}_{\mathit{i,k}\mathrm{-}\mathrm{2}}^{\mathrm{\left(}\mathit{m}\mathrm{\right)}})\\ & \mathrm{+}& \frac{\mathrm{16}}{\mathrm{60}}\hspace{0.17em}({\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}}_{\mathit{i}\mathrm{+}\mathrm{1}\mathit{,k}}^{\mathrm{\left(}\mathit{m}\mathrm{\right)}}\mathrm{+}{\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{,k}}^{\mathrm{\left(}\mathit{m}\mathrm{\right)}}\mathrm{+}{\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}}_{\mathit{i,k}\mathrm{+}\mathrm{1}}^{\mathrm{\left(}\mathit{m}\mathrm{\right)}}\mathrm{+}{\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}}_{\mathit{i,k}\mathrm{-}\mathrm{1}}^{\mathrm{\left(}\mathit{m}\mathrm{\right)}})\\ & \mathrm{-}& \end{array}$(12)where h is the spacing of grid points. The divergence of is calculated with fourth-order accurate finite differencing method for each grid point, and for points on the boundary of the grid and the neighboring row and column, a second-order accurate finite differencing scheme is employed. Convergence is reached when two requirements are met: (i) at least iterations have been made; and (ii) the maximum difference $\mathrm{(}{\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\psi }\end{array}}_{\mathit{i,k}}^{\mathrm{\left(}\mathit{m}\mathrm{\right)}}\mathrm{-}{\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\psi }\end{array}}_{\mathit{i,k}}^{\mathrm{(}\mathit{m}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathrm{)}}_{\max }$ between two iterations increases. Typically, slightly more than are needed to reach convergence. If the process converges, the values of $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\psi }\end{array}$ at the four corners are calculated by extrapolation.

Fig. 2 Illustration of the extrapolation method used in the SOR method (Sect. 3.2) to calculate $\begin{array}{c}{{˜}}_{}\\ {\alpha }\end{array}$ is shown: Based on the calculation of $\begin{array}{c}{{˜}}_{}\\ {\alpha }\end{array}$ on two grids with indices (i,k) and (I,K), we can retrieve ${\begin{array}{c}{{˜}}_{}\\ {\alpha }\end{array}}^{\mathrm{true}}$ with a minimum accuracy Δα using the scheme described in the figure.

We consider that the typical accuracy on the image position of observed lens systems is of the order 5 mas, implying that ε_acc in Eq. (9) should be of the same order (this choice will be discussed in Sect. 4). Thus, the numerical error of our method has to be well below 1 mas ≈ 10^-3θ_E for typical galaxy scale lenses which is quite stringent. Increasing the grid size yields a strong increase in computational time, which scales roughly as . Therefore, we added an extrapolation method to the standard SOR to increase accuracy with a more reasonable increase in computational time. The principle of our extrapolation scheme is displayed in Fig. 2 and is based on the observation that the error | Δα | of the computed value $\begin{array}{c}{}_{{˜}}\\ {\alpha }\end{array}\mathrm{\left(}\mathit{h}\mathrm{\right)}$ scales as . This can be seen in the top panel of Fig. 3 where we applied our numerical scheme to the case of a non-singular isothermal sphere, i.e., where the true solution is known analytically. In this case, the deflection law is a pure gradient field, and thus ${\alpha }\mathrm{=}{\alpha ̂}\mathrm{=}\begin{array}{c}{}_{{˜}}\\ {\alpha }\end{array}$. Using this scaling behavior we can extrapolate to the true deflection ${\begin{array}{c}{{˜}}_{}\\ {\alpha }\end{array}}^{\mathrm{true}}$, which would be obtained in the limit h → 0, for every grid point ${\begin{array}{c}{˜}\\ {\alpha }\end{array}}_{\mathit{i,k}}\mathrm{\left(}\mathit{h}\mathrm{\right)}\mathrm{=}{\begin{array}{c}{˜}\\ {\alpha }\end{array}}_{\mathit{i,k}}^{\mathrm{true}}\mathrm{+}{\mathit{E}}_{\mathit{i,k}}^{\mathrm{num}}\mathrm{(}\mathit{h}{\mathrm{)}}^{\mathrm{2}}\mathit{,}$(13)where E^num is the numerical error³. However, the asymptotic deflection ${\begin{array}{c}{{˜}}_{}\\ {\alpha }\end{array}}^{\mathrm{true}}$ and the value of the numerical error E^num are unknown in general. We can determine the two unknowns by calculating $\begin{array}{c}{{˜}}_{}\\ {\alpha }\end{array}$ for two different values of h, i.e., for different . Hence, we calculate $\begin{array}{c}{{˜}}_{}\\ {\alpha }\end{array}$ on two grids, of and points. The coordinates of the first and second grid are denoted respectively with indices (I,K) and (i,k) and we have to match every grid point (i,k) with its corresponding position (I,K). Then we can obtain the true value ${\begin{array}{c}{{˜}}_{}\\ {\alpha }\end{array}}^{\mathrm{true}}$${\begin{array}{c}{˜}\\ {\alpha }\end{array}}_{\mathit{i,k}}^{\mathrm{true}}\mathrm{=}\frac{\mathrm{4}{\begin{array}{c}{˜}\\ {\alpha }\end{array}}_{\mathit{I,K}}\mathrm{(}{\frac{\mathit{h}}{\mathrm{2}}}^{\mathrm{)}}\mathrm{-}{\begin{array}{c}{˜}\\ {\alpha }\end{array}}_{\mathit{i,k}}\mathrm{\left(}\mathit{h}\mathrm{\right)}}{\mathrm{3}}\mathit{,}$(14)as indicated in Fig. 2.

Incorporating this extrapolation method in the code decreases the numerical error for the grid point numbers that are used () roughly by a factor of 10³, as can be seen in the lower panel of Fig. 3, which also shows that the numerical error with this extrapolation scheme decreases much faster with than without. The largest numerical deviation (Δα)_max can be found near the center of the grid. This is as expected, since the deviation from ${\begin{array}{c}{{˜}}_{}\\ {\alpha }\end{array}}^{\mathrm{true}}$ depends on higher-order derivatives, which for the chosen lens model are largest near the center. However, multiple images near the center of the lens are usually highly demagnified and rarely observable for galaxies as lenses (see, e.g., Hezaveh et al. 2015; Winn et al. 2004) and are therefore not relevant.

Fig. 3 Maximum difference of $\mathrm{|}{\alpha ̂}\mathrm{-}\begin{array}{c}{}_{{˜}}\\ {\alpha }\end{array}\mathrm{|}\mathrm{=}\mathrm{\left|}\mathrm{\Delta }{\alpha }\mathrm{\right|}$ for a non-singular isothermal sphere with core radius θc = 0.1 θE as a function of the number of grid points number used in the numerical solution. Blue dots are the numerical results, whereas the curves present power-law fits to these points with h being the spacing of grid points. In the top panel, the results are shown for the “standard” method, where the numerical error scales as . Incorporating the extrapolation scheme, the numerical error decreases much faster with the number of grid point, as can be seen in the lower panel (note the different scale for the y-axis in the upper and lower panel). For the typical values used in the papers (), a gain in accuracy by three orders of magnitude is obtained with extrapolation with only a modest increase of computational cost (~ 25%). Since extrapolation includes calculating $\begin{array}{c}{{˜}}_{}\\ {\alpha }\end{array}$ twice with grid points 2N and N, values are only shown up to which correspond to in the graph above.

We have also tested the accuracy of the numerical implementation. One method is to check whether $\mathrm{\nabla }\mathrm{·}{\alpha ̂}\mathrm{\left(}{\theta }\mathrm{\right)}\mathrm{=}\mathrm{\nabla }\mathrm{·}\begin{array}{c}{}_{{˜}}\\ {\alpha }\end{array}\mathrm{\left(}{\theta }\mathrm{\right)}$ is valid for the whole grid for any deformation function . In all our calculations, deviations $\mathrm{|}\mathit{\kappa ̂}\mathrm{-}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\kappa }\end{array}\mathrm{|}$ are always smaller than 10^-4 if the corner regions, i.e., θ ≥ 2 θ_E, are excluded from our analysis. Thus, we only consider the behavior of the mass profile in the circular region | θ | /θ_E< 2, where numerical errors in $\mathrm{|}{\alpha ̂}\mathrm{-}\begin{array}{c}{}_{{˜}}\\ {\alpha }\end{array}\mathrm{|}$ do not exceed 10^-6θ_E.

3.3. Solution by means of a Green’s function

A different approach is to find a solution of Eq. (8) by means of a Green’s function. For that, we make use of Green’s second theorem, considering a function h(θ), ${\mathrm{\int }}_{\mathrm{𝒰}}{\mathrm{d}}^{\mathrm{2}}\mathit{\theta }\mathrm{[}\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\hspace{0.17em}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{h}\mathrm{-}\mathit{h}\hspace{0.17em}{\mathrm{\nabla }}^{\mathrm{2}}{\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}}^{\mathrm{]}}\mathrm{=}{\mathrm{\int }}_{\mathit{\partial }\mathrm{𝒰}}\mathrm{d}\mathit{s}\mathrm{[}\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\hspace{0.17em}\mathrm{\nabla }\mathit{h}\mathrm{·}{n}\mathrm{-}\mathit{h}\hspace{0.17em}\mathrm{\nabla }\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\mathrm{·}{{n}}^{\mathrm{]}}\mathit{.}$(15)We choose the function h(θ) ≡ H(ϑ;θ) depending on the vector ϑ such that it obeys the following equations: $\begin{array}{ccc}{\mathrm{\nabla }}_{\mathit{\theta }}^{\mathrm{2}}\mathit{H}\mathrm{\left(}{\vartheta }\mathrm{;}{\theta }\mathrm{\right)}\mathrm{=}\mathit{\delta }\mathrm{(}{\theta }\mathrm{-}{\vartheta }\mathrm{)}\mathrm{-}\frac{\mathrm{1}}{\mathit{A}}& \mathrm{for}& {\theta }\mathrm{\in }\mathrm{𝒰}\mathit{,}\\ {\mathrm{\nabla }}_{\mathit{\theta }}\mathit{H}\mathrm{\left(}{\vartheta }\mathrm{;}{\theta }\mathrm{\right)}\mathrm{·}{n}\mathrm{=}\mathrm{0}& \mathrm{for}& \end{array}$(16)where A is the area of , and ϑ is a point within . The term 1 /A in Eq. (16) is needed to satisfy Green’s divergence theorem applied to the vector field ∇h, which requires ${\mathrm{\int }}_{\mathrm{𝒰}}{\mathrm{d}}^{\mathrm{2}}\mathit{\theta }{\mathrm{\nabla }}^{\mathrm{2}}\mathit{h}\mathrm{=}{\mathrm{\int }}_{\mathit{\partial }\mathrm{𝒰}}\mathrm{d}\mathit{s}\mathrm{\nabla }\mathit{h}\mathrm{·}{n}\mathrm{;}$(17)the conditions (16) set both side of this equation to zero. With (16), Eq. (15) simplifies to $\begin{array}{ccc}\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\mathrm{\left(}{\vartheta }\mathrm{\right)}& \mathrm{=}& \mathrm{⟨}{\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}}^{\mathrm{⟩}}\mathrm{+}{\mathrm{\int }}_{\mathrm{𝒰}}{\mathrm{d}}^{\mathrm{2}}\mathit{\theta H}\mathrm{\left(}{\vartheta }\mathrm{;}{\theta }\mathrm{\right)}\hspace{0.17em}{\mathrm{\nabla }}^{\mathrm{2}}\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\mathrm{-}{\mathrm{\int }}_{\mathit{\partial }\mathrm{𝒰}}\mathrm{d}\mathit{sH}\mathrm{\left(}{\vartheta }\mathrm{;}{\theta }\mathrm{\right)}\hspace{0.17em}\mathrm{\nabla }\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\mathrm{·}{n}\\ & \mathrm{=}& \end{array}$(18)where ${}^{\mathrm{⟨}}\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\psi }\end{array}^{\mathrm{⟩}}$ is the average of $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\psi }\end{array}$ on , and we used Eq. (8) in the last step. We note that the integral $\mathit{f}\mathrm{\left(}{\theta }\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{𝒰}}{\mathrm{d}}^{\mathrm{2}}\mathit{\vartheta H}\mathrm{\left(}{\vartheta }\mathrm{;}{\theta }\mathrm{\right)}$(19)is a constant, since ∇²f(θ) = 0 and n·∇f = 0 on the boundary of . Therefore, if we integrate Eq. (18) over , the two integrals on the r.h.s. compensate each other, due to the divergence theorem, so that this solution is consistent.

Whereas for a general region it will be difficult to obtain a solution of Eq. (16) for H(ϑ;θ), such a solution is analytically known if is a circle of radius R. In this case, $\begin{array}{ccc}\mathit{H}\mathrm{\left(}{\vartheta }\mathrm{;}{\theta }\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\left[\ln \frac{{\left|{\vartheta }\mathrm{-}{\theta }\right|}^{\mathrm{2}}}{{\mathit{R}}^{\mathrm{2}}}\mathrm{+}\ln \left(\mathrm{1}\mathrm{-}\frac{\mathrm{2}{\vartheta }\mathrm{·}{\theta }}{{\mathit{R}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{\left|}{\vartheta }{\mathrm{|}}^{\mathrm{2}}\mathrm{\right|}{\theta }{\mathrm{|}}^{\mathrm{2}}}{{\mathit{R}}^{\mathrm{4}}}\right)\right]\\ & & \end{array}$(20)which has the properties that

$\begin{array}{ccc}& & {\mathrm{\nabla }}_{\mathit{\vartheta }}^{\mathrm{2}}\mathit{H}\mathrm{\left(}{\vartheta }\mathrm{;}{\theta }\mathrm{\right)}\mathrm{=}\mathit{\delta }\mathrm{(}{\vartheta }\mathrm{-}{\theta }\mathrm{)}\mathrm{-}\frac{\mathrm{1}}{\mathit{\pi }{\mathit{R}}^{\mathrm{2}}}\mathrm{=}{\mathrm{\nabla }}_{\mathit{\theta }}^{\mathrm{2}}\mathit{H}\mathrm{\left(}{\vartheta }\mathrm{;}{\theta }\mathrm{\right)}\mathrm{for}{\vartheta }\mathit{,}{\theta }\mathrm{\in }\mathrm{𝒰}\mathit{,}\\ & & {\mathrm{\nabla }}_{\mathit{\theta }}\mathit{H}\mathrm{\left(}{\vartheta }\mathrm{;}{\theta }\mathrm{\right)}\mathrm{·}{n}\mathrm{\left(}{\theta }\mathrm{\right)}\mathrm{=}\mathrm{0}\mathrm{for}{\theta }\mathrm{\in }\mathit{\partial }\mathrm{𝒰}\mathit{.}\end{array}$Hence, Eq. (20) indeed satisfies the conditions (16).

With this explicit solution, we can now calculate the deflection angle corresponding to the potential $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\psi }\end{array}$, $\begin{array}{c}{}_{\mathrm{˜}}\\ {\alpha }\end{array}\mathrm{=}\mathrm{\nabla }\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\psi }\end{array}$, by obtaining the gradient of H, ${\mathrm{\nabla }}_{\mathit{\vartheta }}\mathit{H}\mathrm{\left(}{\vartheta }\mathrm{;}{\theta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\left(\frac{{\vartheta }\mathrm{-}{\theta }}{\mathrm{|}{\vartheta }\mathrm{-}{\theta }{\mathrm{|}}^{\mathrm{2}}}\mathrm{-}\frac{{\vartheta }}{{\mathit{R}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{|}{\theta }{\mathrm{|}}^{\mathrm{2}}{\vartheta }\mathrm{-}{\mathit{R}}^{\mathrm{2}}{\theta }}{\left({\mathit{R}}^{\mathrm{4}}\mathrm{-}\mathrm{2}{\mathit{R}}^{\mathrm{2}}{\vartheta }\mathrm{·}{\theta }\mathrm{+}\mathrm{\left|}{\vartheta }{\mathrm{|}}^{\mathrm{2}}\mathrm{\right|}{\theta }{\mathrm{|}}^{\mathrm{2}}\right)}\right)\mathrm{·}$(21)Then, $\begin{array}{c}{˜}\\ {\alpha }\end{array}\mathrm{\left(}{\vartheta }\mathrm{\right)}\mathrm{=}\mathrm{2}{\mathrm{\int }}_{\mathrm{𝒰}}{\mathrm{d}}^{\mathrm{2}}\mathit{\theta }{\mathrm{\nabla }}_{\mathit{\vartheta }}\mathit{H}\mathrm{\left(}{\vartheta }\mathrm{;}{\theta }\mathrm{\right)}\hspace{0.17em}\mathit{\kappa ̂}\mathrm{\left(}{\theta }\mathrm{\right)}\mathrm{-}{\mathrm{\int }}_{\mathit{\partial }\mathrm{𝒰}}\mathrm{d}\mathit{s}{\mathrm{\nabla }}_{\mathit{\vartheta }}\mathit{H}\mathrm{\left(}{\vartheta }\mathrm{;}{\theta }\mathrm{\right)}\hspace{0.17em}\hat{{\alpha }}\mathrm{·}{n}\mathit{,}$(22)where we have to deal with a pole in the first term of Eq. (21). Using a conformal mapping as described in Appendix A, we can handle this pole numerically. In the third term the pole lies outside the circle and since there is no pole. However, if θ is on the circle (as occurs in the line integral in Eq. (22)), the third term can become quite large; hence, for points ϑ near the boundary, special care is needed to obtain an accurate solution.

This Green’s function approach has several advantages over using a SOR for a grid. First, the region on which the Neumann problem is defined can be chosen as a circle, instead of a rectangle for the SOR method. Second, the solution by means of the Green’s function yields higher accuracy. The reason for this is that the limited accuracy in finite differencing does not occur here. Third, if one is interested in the deflection only at a few points, this can be calculated much faster than with the SOR which necessarily calculated the solution over the whole region.

3.4. Interpretation

The expression (18) for the deflection potential, or the expression (22) for the deflection angle, contains quite a number of terms. In order to obtain a better understanding of the various terms, we consider again the case where the deflection angle $\hat{{\alpha }}$ is a pure gradient field, in which case $\begin{array}{c}{}_{\mathrm{˜}}\\ {\alpha }\end{array}\mathrm{=}\hat{{\alpha }}\mathrm{\equiv }{\alpha }$. Then the deflection angle at a point can be decomposed into a deflection α_in which is caused by matter inside , and one due to matter outside , denoted by α_out. Thus we expect that ${\alpha }\mathrm{\left(}{\vartheta }\mathrm{\right)}\mathrm{=}{{\alpha }}_{\mathrm{in}}\mathrm{\left(}{\vartheta }\mathrm{\right)}\mathrm{+}{{\alpha }}_{\mathrm{out}}\mathrm{\left(}{\vartheta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{𝒰}}{\mathrm{d}}^{\mathrm{2}}\mathit{\theta \kappa }\mathrm{\left(}{\theta }\mathrm{\right)}\hspace{0.17em}\frac{{\vartheta }\mathrm{-}{\theta }}{\mathrm{|}{\vartheta }\mathrm{-}{\theta }{\mathrm{|}}^{\mathrm{2}}}\mathrm{+}{{\alpha }}_{\mathrm{out}}\mathrm{\left(}{\vartheta }\mathrm{\right)}\mathit{.}$(23)Comparing the last Eqs. (23) to (22), we find that ${\alpha }\mathrm{\left(}{\vartheta }\mathrm{\right)}\mathrm{=}{{\alpha }}_{\mathrm{in}}\mathrm{\left(}{\vartheta }\mathrm{\right)}\mathrm{+}{A}\mathrm{\left(}{\vartheta }\mathrm{\right)}\mathrm{-}{{B}}_{\mathrm{in}}\mathrm{\left(}{\vartheta }\mathrm{\right)}\mathrm{-}{{B}}_{\mathrm{out}}\mathrm{\left(}{\vartheta }\mathrm{\right)}\mathit{,}$(24)where

$\begin{array}{ccc}& & {A}\mathrm{\left(}{\vartheta }\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{𝒰}}{\mathrm{d}}^{\mathrm{2}}\mathit{\theta \kappa }\mathrm{\left(}{\theta }\mathrm{\right)}\hspace{0.17em}\left(\mathrm{2}{\mathrm{\nabla }}_{\mathit{\vartheta }}\mathit{H}\mathrm{\left(}{\vartheta }\mathrm{;}{\theta }\mathrm{\right)}\mathrm{-}\frac{\mathrm{1}}{\mathit{\pi }}\frac{{\vartheta }\mathrm{-}{\theta }}{\mathrm{|}{\vartheta }\mathrm{-}{\theta }{\mathrm{|}}^{\mathrm{2}}}\right)\mathit{,}\\ & & \end{array}$(25)where we split the deflection angle on the boundary into terms due to matter inside and outside . Both of the terms A and B_in are due to matter inside , whose deflection is covered entirely by the first term α_in, so that we expect that ${A}\mathrm{\left(}{\vartheta }\mathrm{\right)}\mathrm{=}{{B}}_{\mathrm{in}}\mathrm{\left(}{\vartheta }\mathrm{\right)}\mathit{.}$(26)In Appendix B we show explicitly that this relation holds for the case of a circular region for which H is given by Eq. (21). Hence, Eq. (24) then provides a clean separation of the deflection angle coming from the inner mass distribution (α_in) and that coming from matter outside , given by B_out. This relation may be of practical relevance for the numerical calculation of the lensing properties from a complicated mass distribution, for which the lensing quantities are only needed inside a limited region. Instead of calculating, for every point inside , a two-dimensional integral of the surface mass density κ over the whole lens plane, one can proceed as follows: first, one can reduce the integration range over the region to get the contribution α_in. Second, one can calculate the contribution α_out for points on the boundary by integrating over the outer region of the lens in terms of a two-dimensional integral. Third, the contribution α_out for points inside the region can then be obtained by a one-dimensional integration over the boundary curve.

In general, if α is given on the boundary, it contains contributions from both the inner and the outer part. In other words, the split of B into B_in and B_out is not provided in that case. The term A then compensates for the contribution B_in of B.

4. Illustrative example – a quadrupole lens and an isotropic SPT

Our goal is to find criteria allowing us to assess whether an SPT-transformed deflection law is valid (i.e. deviates from a gradient field by less than ε_acc), using the methods explained in the previous section. Thus, we set an upper limit on how much the transformed deflection law is allowed to differ from its closest curl-free approximation $\begin{array}{c}{{˜}}_{}\\ {\alpha }\end{array}$ before leading to a non-negligible shift of the lensed images. Since observed lens systems are usually fit by simple mass models with only a small number of free parameters, we do not expect the fit to be perfect. We always have to deal with observational uncertainties as well as the presence of substructure (Xu et al. 2010; Bradač et al. 2004; Kochanek & Dalal 2004; Mao & Schneider 1998) and line-of-sight inhomogeneities (Xu et al. 2012; Metcalf 2005). Therefore, we cannot reproduce observed positions better than a few milliarcseconds with a smooth mass model. Hence, as long as | Δα(θ) | is less than the smallest angular scale on which modeling with a smooth mass model is still meaningful, differences are of no practical relevance (SS14).

We need to choose a lens model to explore how seriously the SPT may affect lens modeling. First, we consider a situation similar to SS14, namely a quadrupole lens with external shear γ_p${\alpha }\mathrm{=}\mathit{\kappa ̅}\mathrm{\left(}\mathrm{\right|}{\theta }\mathrm{\left|}\mathrm{\right)}\hspace{0.17em}{\theta }\mathrm{-}\left(\begin{array}{c}\\ {\mathit{\gamma }}_{\mathrm{p}}& \mathrm{0}\\ \mathrm{0}& \mathrm{-}{\mathit{\gamma }}_{\mathrm{p}}\end{array}\right){\theta }$(27)which is deformed by an SPT corresponding to a radial stretching, as in Eq. (5). Specifically, we choose ${\beta ̂}\mathrm{\left(}{\beta }\mathrm{\right)}\mathrm{=}\left(\mathrm{1}\mathrm{+}\frac{{\mathit{f}}_{\mathrm{2}}}{\mathrm{2}{{\mathit{\theta }}_{\mathrm{E}}}^{\mathrm{2}}}{\mathit{\beta }}^{\mathrm{2}}\right){\beta }\mathit{.}$(28)This deformation function is the lowest-order expansion of more general stretching functions, and its leading-order term is chosen such as to not yield an MST, to cleanly separate the effect of the MST from that of the more general SPT in this study. Furthermore, we choose as specific lens model a non-singular isothermal sphere (NIS), described by the mean convergence profile $\mathit{\kappa ̅}\mathrm{=}{\mathit{\theta }}_{\mathrm{E}}\hspace{0.17em}\frac{\mathrm{1}}{\sqrt{{{\mathit{\theta }}_{\mathrm{c}}}^{\mathrm{2}}\mathrm{+}{\mathit{\theta }}^{\mathrm{2}}}}\mathit{,}$(29)where θ_c is the core radius. For the rest of this paper, we fix the core to be θ_c = 0.1θ_E.

To get a quantitative estimate on how large deviations of $\begin{array}{c}{\mathrm{˜}}_{}\\ {\alpha }\end{array}$ from $\hat{{\alpha }}$ are tolerable before the lensing properties of the SPT deviate markedly from the original lens model, we take the Hubble Space Telescope (HST) as example. We estimate that the highest astrometric accuracy that can be achieved corresponds to about a tenth of a pixel in the ACS camera, or Δθ ≈ 5 mas ≈ 5 × 10^-3θ_E, where the last expression accounts for the fact that the typical Einstein radii of galaxy-scale lenses are of order one arcsecond. Thus, if the solution $\begin{array}{c}{\mathrm{˜}}_{}\\ {\alpha }\end{array}$ satisfies the condition (9) with ϵ_acc = 5 × 10^-3θ_E over the region | θ | ≤ 2 θ_E, we call the corresponding SPT allowed or valid.

4.1. Impact on the deflection law

The model we consider has two free parameters, the distortion parameter f₂ in the SPT (28), and the strength γ_p of the external shear. We start with exploring this parameter space to find the combination that yield allowed transformations, using the methods described in the previous section. In Fig. 4, we display the maximum deviation |Δα| _max as a function of these two parameters. It shows a wide range of allowed parameter combinations, where the allowed range of f₂ decreases with increasing external shear. The white regions in Fig. 4 denotes parameter combinations where |Δα| _max > 0.005 θ_E, and which are therefore not allowed according to our accuracy criterion.

Fig. 4 Values of |Δα| max are plotted against the parameters f2 from (28) and external shear strength γp. The colored region indicates allowed pairs of parameters that fulfill the |Δα| < 5 × 10-3θE-criterion. For obtaining this figure, we used the SOR method.

In SS14, we speculated that the curl of may yield a good indication for the deviation of the SPT-transformed deflection field from a gradient field. In this case, the curl , which describes the asymmetric part of the Jacobian, could be used as a proxy for | Δα |. For a quadrupole lens of the form (27) and the deformation law (28), the curl is given in Eq. (42) of SS14, $\begin{array}{ccc}\mathit{\kappa ̂}\mathrm{I}& \mathrm{\approx }& \\ & & \mathrm{\times }\mathrm{[}{{\mathit{\gamma }}_{\mathrm{p}}}^{\mathrm{2}}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\kappa ̅}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{\gamma }}_{\mathrm{m}}\mathrm{+}\mathrm{1}\mathrm{-}\mathit{\kappa ̅}\mathrm{)}\mathrm{+}\mathrm{2}{\mathit{\gamma }}_{\mathrm{m}}{\mathit{\gamma }}_{\mathrm{p}}\cos \mathrm{(}\mathrm{2}\mathit{\varphi }{\mathrm{)}}^{\mathrm{]}}\sin \mathrm{2}\mathit{\varphi ,}\end{array}$(30)where θ,ϕ describe polar coordinates in the lens plane and is the shear caused by the NIS lens.

Fig. 5 Values of are plotted against the parameters f2 from (28) and external shear strength γp. The colored region indicates allowed pairs of parameter that were chosen such that they roughly correspond to | Δα | < 5 × 10-3θE.

Figure 5 shows the maximum of $\mathit{\kappa ̂}{\mathrm{I}}_{}$ as a function of external shear γ_p and deformation “strength” f₂, which indeed is very similar to Fig. 4. The actual difference between those two approaches is seen in Fig. 6. An approximately linear correlation with an expected but modest scatter can be seen. In fact, from that figure we obtain for our specific model that $\mathrm{0.16}\hspace{0.17em}\mathrm{|}\mathit{\kappa ̂}\mathrm{I}{\mathrm{|}}_{\max }\lesssim \frac{\mathrm{|}\mathrm{\Delta }{\alpha }{\mathrm{|}}_{\max }}{{\mathit{\theta }}_{\mathrm{E}}}\lesssim \mathrm{0.3}\hspace{0.17em}\mathrm{|}\mathit{\kappa ̂}\mathrm{I}{\mathrm{|}}_{\max }\mathit{.}$(31)For other models, the relation between | Δα | _max and will be different; nevertheless, we see that the curl of indeed provides a useful indication for the validity of an SPT, since calculating $\mathit{\kappa ̂}{\mathrm{I}}_{}$ is much easier then obtaining the numerical solution for $\begin{array}{c}{{˜}}_{}\\ {\alpha }\end{array}$.

Fig. 6 For every allowed combination f2 and γp the values of | Δα | max (Fig. 4) are plotted against (Fig. 5). A clear correlation between these two quantities can be seen.

Figure 7 illustrates how a specific deflection law in a region | θ | ≤ 2 θ_E is affected by an SPT. It shows | Δα(θ) | for a quadrupole lens with external shear γ_p = 0.1 and deformation strength f₂ = 0.55, which is the highest allowed for this value of the external shear strength and thus is expected to show the largest deviations compared to the original mass profile. The figure shows that the largest deviations occur at an angle of 45° with respect to the external shear. This pattern, which is shown for one specific pair of f₂ and γ_p, is qualitatively the same for all f₂-γ_p-combinations.

Fig. 7 Map of | Δα(θ) | is shown for f2 = 0.55 and γp = 0.1. The strong changes in the corners, i.e. θ > 2 θE, are biased by large numerical uncertainty and should be neglected.

4.2. Implications for the convergence

We show in Fig. 8 the comparison between original (κ) and SPT-transformed mass distribution () for three different allowed pairs of parameters, f₂ = 0.55 and γ_p = 0.1 (the same combination of parameters as in Fig. 7), f₂ = −0.55 and γ_p = 0.1, and f₂ = 1.2 and γ_p = 0.05. The lower panel of Fig. 8 shows the change of the radial profile as .

Fig. 8 Upper panel shows the mass profile of the original NIS lens (solid curve), and that of three SPT-transformed lenses, with parameters f2 and γp indicated by the labels. For all of these three models, Δαmax ≈ εacc = 5 × 10-3θE. Since the transformed mass distributions have a finite ellipticity, the density is plotted as a function of the geometric mean of the major and minor semi-axis of the best-fitting ellipse to an isodensity contour, except for the case with negative f2, for which the outer isodensity contours are not closing around the lens center; in this special case, the x-axis corresponds to the θ1-axis. The convergence changes up to 28% for radii smaller than 1 θE, radii larger than that show a significantly smaller convergence for a positive f2. Negative f2 show an essentially mirrored behavior compared to positive f2. This leads to convergence that may not decrease monotonically. The lower panel shows the ratio between transformed and original mass profile.

The divergence of (i.e. ), was calculated analytically in SS14 (see their Eq. (41)) and can be used to compare our numerical results to the analytic solution. Specialized to our case, it reads $\begin{array}{ccc}\mathit{\kappa ̂}& \mathrm{=}& {\mathit{\kappa }}_{\mathrm{NIS}}\mathrm{+}\frac{{\mathit{f}}_{\mathrm{2}}}{\mathrm{2}}{\left(\frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{E}}}\right)}^{\mathrm{2}}\\ & & \mathrm{\times }({\mathit{\gamma }}_{\mathrm{m}}[\mathrm{2}{\mathit{\gamma }}_{\mathrm{p}}^{\mathrm{2}}\mathrm{+}\mathrm{3}(\mathrm{1}\mathrm{-}\mathit{\kappa ̅}{)}^{\mathrm{2}}]\mathrm{-}\mathrm{2}(\mathrm{1}\mathrm{-}\mathit{\kappa ̅})[(\mathrm{1}\mathrm{-}\mathit{\kappa ̅}{)}^{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathit{\gamma }}_{\mathrm{p}}^{\mathrm{2}}]{}^{}\\ & & \mathrm{+}[\mathrm{5}{\mathit{\gamma }}_{\mathrm{p}}(\mathrm{1}\mathrm{-}\mathit{\kappa ̅}{)}^{\mathrm{2}}\mathrm{-}\mathrm{6}{\mathit{\gamma }}_{\mathrm{p}}{\mathit{\gamma }}_{\mathrm{m}}(\mathrm{1}\mathrm{-}\mathit{\kappa ̅})\mathrm{+}{\mathit{\gamma }}_{\mathrm{p}}^{\mathrm{3}}]\hspace{0.17em}\cos \mathrm{2}\mathit{\varphi }\mathrm{+}{\mathit{\gamma }}_{\mathrm{p}}^{\mathrm{2}}{\mathit{\gamma }}_{\mathrm{m}}\hspace{0.17em}\cos \mathrm{4}\mathit{\varphi }{}^{})\mathit{.}\end{array}$(32)where again θ,ϕ describe polar coordinates in the lens plane. The change is proportional to the stretching parameter f₂, so that Δκ(−f₂) = −Δκ(f₂). This behavior can be seen in Fig. 8. Indeed, we checked that all numerically obtained deflection angles $\begin{array}{c}{{˜}}_{}\\ {\alpha }\end{array}$ are such that their corresponding surface mass densities agree with the analytical prediction (32). For example, the numerical result for the parameter combination γ_p = 0.1 and f₂ = 0.55 deviates by less than 3 × 10^-3 from the analytical solution.

As seen from Eq. (32), the resulting mass distribution is no longer axi-symmetric, but that isodensity contours are nearly elliptical (i.e., a factor proportional to cos(2φ)) with a small boxiness (i.e., a factor proportional to cos(4φ)). Hence, we define the distance from the center generally as the geometric mean $\sqrt{{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}}$ using the 1- and 2-axis of the elliptical isodensity contours. However, for sufficiently negative f₂, the outer isodensity contours are no longer concentric, i.e., they are not closed curves around the center of the lens. In addition, for large negative values of f₂, the radial profile can become non-monotonic. We consider such a behavior as non-physical, i.e., such resulting models will be irrelevant in practice.

The ellipticity of the transformed mass profiles is non-negligible as shown in Fig. 9 where ϵ, defined as the axis ratio 1- over 2-axis, is plotted as a function of radius. Integrating the analytic representation (32) of up to 1 θ_E it can be shown that the mass enclosed within the Einstein radius of the original lens is conserved, independent of the chosen mass profile κ(θ) (see Appendix C).

Fig. 9 Radial dependence of the axis ratio ϵ. In the unperturbed case the isodensity contours are circular, i.e., $\mathit{ϵ}\mathrm{\left(}\sqrt{{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}}\mathrm{\right)}\mathrm{=}\mathrm{1}$. The SPTed mass distribution shows for radii $\sqrt{{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}}\mathit{<}\mathrm{1}\hspace{0.17em}{\mathit{\theta }}_{\mathrm{E}}$ deviations of up to 5% from circularity, whereas for larger radii the deviations can be up to 20%. The convergence map for f2 = −0.55 and γp = 0.1 does not show concentric isodensity contours for radii larger than 1.3 θE and therefore no ellipticity as a function radius can be determined.

5. Characterization of the modified mass distribution

The SPT leads to a modified deflection angle of the lens which yields exactly the same astrometric and photometric observational properties as the original mass distribution. For those modified profiles $\hat{{\alpha }}$ for which a deflection potential $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\psi }\end{array}$ can be found such that the differences between the corresponding Δα is sufficiently small, the modified surface mass density provides a viable alternative to the original mass model κ of the lens. In this section we want to consider a diagnostic for the change of the mass profile, both regarding the radial slope and the angular structure of the lens. Since the strong lensing properties of the lens can only be probed in the inner part of the mass distribution, we will apply these diagnostics only to those regions where multiple images can occur, i.e., | θ | ≲ 2θ_E.

5.1. Radial mass profile

The SPT changes the radial mass profile of the lens. We consider situations in which the original lens is described by a “simple” mass distribution, i.e., an NIS. Combined with a “mild” SPT the resulting remains simple, e.g., still shows closed, concentric isodensity contours.

The mass-sheet transformation is a special case of the SPT, and it is well known that the MST changes the radial profile of the lens. In order to highlight the new feature of the SPT not contained in the MST, we aim at a measure for the radial profile which is invariant under the MST. The MST transforms all derivatives of κ by a constant factor λ, hence it leaves the ratio of derivatives unchanged. Consequently, one possible diagnostic for the effect of the SPT is the radial profile of such ratios, e.g., ⟨κ⟩′′/⟨κ⟩′.

In particular, if the original mass profile is a power law, ⟨κ⟩(θ) ∝ θ^{− s}, then we have θ ⟨κ⟩′′/⟨κ⟩′ = −(s + 1); hence, any deviation from this constant value indicates the effect of the SPT on the modified mass profile . However, if there is no analytical expression of κ and , the ratio of derivatives is very sensitive to numerical noise, and therefore of little practical interest. We therefore consider hereafter alternative tests.

Noting that the MST yields a multiplication of 1−κ(θ) by λ, the ratio ${\mathit{R}}_{\mathit{\kappa }}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}\mathrm{-}⟨\mathit{\kappa }⟩\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{⟨\mathit{\kappa }⟩}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}\mathit{,}$(33)is well defined for monotonically decreasing mass profiles and invariant under the MST. Figure 10 shows R_κ for the NIS and various SPT transformed models. The variations of R_κ are particularly significant above one Einstein radius, in regions where the SPT-transformed profiles deviate also more strongly from the original profile. Despite the fact that deviates from κ_orginal by more than 20% within θ_E when f₂ = 1.2 (Fig. 8), the most extreme changes of R_κ(θ) reach no more than ~10% within one Einstein radius. As expected, R_κ deviates more strongly from the original profile when | f₂ | is large. Negative values of f₂ (not shown) are qualitatively similar (but mirrored w.r.t. ${\mathit{R}}_{\mathit{\kappa }}^{\mathrm{NIS}}$) to the situation encountered for positive f₂. However, at radii θ ~ 1.2 θ_E, stops being monotonically decreasing, and R_κ diverges.

Fig. 10 Top: quantity Rκ(θ) (Eq. (33)) calculated for an NIS with external shear γp (cf. Sect. 4; solid black) and for various SPT-transformed models with SPT of the form 1 + f2/ 2 (β/θE)2 (Eq. (28)). The range of positive values of f2 allowed by | Δ αmax | < 5 × 10-3θE (Fig. 4) is explored for two different choices of the shear: γp = 0.05 (blue) and γp = 0.1 (red). While Rκ is conserved under an MST, it is not under an SPT, with deviation that can reach tens of percents. Bottom: for each curve of the top panel, we show the difference between Rκ of the original NIS model and of the SPT transformed model.

Another possibility to characterize the radial profile change is through the aperture mass (see Schneider 1996). Consider a function $\mathit{U}\mathrm{\left(}\mathit{\theta }\mathrm{;}{\mathit{\theta }}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}{\mathit{\theta }}_{\mathrm{0}}^{-2}\hspace{0.17em}\mathit{u}\mathrm{(}\mathit{\theta }\mathit{/}{\mathit{\theta }}_{\mathrm{0}}\mathrm{)}$ such that u(x) is non-zero only for x ≤ 1; hence, θ₀ characterizes the range of support of U(θ;θ₀). Furthermore, we require that the filter function U has a vanishing two-dimensional integral over its support, which means that $\begin{array}{ccc}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\mathrm{d}\mathit{xxu}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathrm{0.}& & \end{array}$Then we define the aperture mass as ${\mathit{M}}_{\mathrm{ap}}\mathrm{\left(}{\mathit{\theta }}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}\mathit{\theta \kappa }\mathrm{\left(}{\theta }\mathrm{\right)}\hspace{0.17em}\mathit{U}\mathrm{\left(}\mathrm{\right|}{\theta }\mathrm{\left|}\mathrm{;}{\mathit{\theta }}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{\int }\mathrm{d}\mathit{\theta \theta }\hspace{0.17em}⟨\mathit{\kappa }⟩\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\hspace{0.17em}\mathit{U}\mathrm{\left(}\mathit{\theta }\mathrm{;}{\mathit{\theta }}_{\mathrm{0}}\mathrm{\right)}\mathit{.}$(34)The mass-sheet transformation (1) leads to multiplication of M_ap by a factor λ, whereas the additive term in Eq. (1) drops out, due to the compensated nature of the filter function U. Thus, if we consider the ratio of the aperture mass for two different scale lengths θ₀, the factor λ drops out, and this ratio M_ap(θ₁) /M_ap(θ₂) is invariant under the MST.

Consider again a power-law density profile, ⟨κ⟩(θ) = (1−s/ 2)(θ/θ_E)^{− s}, with 0 <s< 2, where θ_E is the Einstein radius in case of axi-symmetry. Then, $\begin{array}{ccc}{\mathit{M}}_{\mathrm{ap}}\mathrm{=}\mathrm{(}\mathrm{2}\mathrm{-}\mathit{s}\mathrm{)}\mathit{\pi }{\left(\frac{{\mathit{\theta }}_{\mathrm{0}}}{{\mathit{\theta }}_{\mathrm{E}}}\right)}^{\mathrm{-}\mathit{s}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\mathrm{d}\mathit{x}{\mathit{x}}^{\mathrm{1}\mathrm{-}\mathit{s}}\hspace{0.17em}\mathit{u}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{,}& & \end{array}$and M_ap(θ₁) /M_ap(θ₂) = (θ₁/θ₂)^{− s}. We thus define the effective slope ${\mathit{s}}_{\mathrm{eff}}\mathrm{:}\mathrm{=}\frac{\ln \mathrm{\left[}{\mathit{M}}_{\mathrm{ap}}\mathrm{\right(}{\mathit{\theta }}_{\mathrm{1}}\mathrm{)}\mathit{/}{\mathit{M}}_{\mathrm{ap}}\mathrm{(}{\mathit{\theta }}_{\mathrm{2}}{\mathrm{\left)}}^{\mathrm{\right]}}}{\ln \mathrm{(}{\mathit{\theta }}_{\mathrm{2}}\mathit{/}{\mathit{\theta }}_{\mathrm{1}}\mathrm{)}}\mathit{,}$(35)so that for a mass profile of the form ⟨κ⟩(θ) = a + bθ^{− s}, s_eff = s.

One can think of a number of appropriate weight functions u(x) and aperture scales θ_i to characterize the modified mass profile. The simplest form would be the sum of two delta functions, u(x) = δ(x−x₀)−x₀δ(x−1), with x₀< 1, for which M_ap(θ₀) = 2πx₀[⟨κ⟩(x₀θ₀)−⟨κ⟩(θ₀)]. Furthermore, choosing θ₂ = θ₁/x₀, the ratio of aperture masses becomes $\frac{{\mathit{M}}_{\mathrm{ap}}\mathrm{\left(}{\mathit{\theta }}_{\mathrm{1}}\mathrm{\right)}}{{\mathit{M}}_{\mathrm{ap}}\mathrm{(}{\mathit{\theta }}_{\mathrm{1}}\mathit{/}{\mathit{x}}_{\mathrm{0}}\mathrm{)}}\mathrm{=}\frac{⟨\mathit{\kappa }⟩\mathrm{\left(}{\mathit{x}}_{\mathrm{0}}{\mathit{\theta }}_{\mathrm{1}}\mathrm{\right)}\mathrm{-}⟨\mathit{\kappa }⟩\mathrm{\left(}{\mathit{\theta }}_{\mathrm{1}}\mathrm{\right)}}{⟨\mathit{\kappa }⟩\mathrm{\left(}{\mathit{\theta }}_{\mathrm{1}}\mathrm{\right)}\mathrm{-}⟨\mathit{\kappa }⟩\mathrm{(}{\mathit{\theta }}_{\mathrm{1}}\mathit{/}{\mathit{x}}_{\mathrm{0}}\mathrm{)}}\mathrm{·}$(36)In the case of (1−x₀) ≪ 1, the expression (35) becomes ${\mathit{s}}_{\mathrm{eff}}\mathrm{=}\mathrm{-}\mathrm{1}\mathrm{-}\frac{{\mathit{\theta }}_{\mathrm{1}}{⟨\mathit{\kappa }⟩}^{\mathrm{\prime \prime }}\mathrm{\left(}{\mathit{\theta }}_{\mathrm{1}}\mathrm{\right)}}{{⟨\mathit{\kappa }⟩}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{\theta }}_{\mathrm{1}}\mathrm{\right)}}\mathrm{+}\mathrm{𝒪}\mathrm{\left(}\mathrm{\right[}\mathrm{1}\mathrm{-}{\mathit{x}}_{\mathrm{0}}{\mathrm{]}}^{\mathrm{2}}\mathrm{)}\mathit{.}$(37)Hence, we see that in this case s_eff depends just on the ratio of second to first derivative, and reduces to s for a power-law mass profile with slope s.

More practical choices of u would be such that the profile is probed over an annulus around the Einstein radius θ_E. For example, one could use the compensated filter function $\mathit{u}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{for}{\mathit{x}}_{\mathrm{0}}\mathrm{\le }\mathit{x}\mathrm{\le }\mathrm{1}\\ & \mathrm{else}\mathit{.}\end{array}$(38)Figure 11 shows M_ap as a function of θ₀, fixing x₀ = 1/2 in Eq. (38). As expected, for two profiles transformed into each other via an MST, the ratio of aperture masses is independent of θ₀ and equals λ. Conversely, M_ap(θ₀) of the SPT-transformed profiles intersects the aperture mass “function” of the original profile (i.e. NIS), and the ratio between the two curves changes with θ₀. The radius at which the curves intersects is almost independent of the value of f₂. This can easily be deduced from the apparent self-similarity of the SPT-transformed mass density profiles (Fig. 8) for various values of f₂.

Figure 11 motivates a choice of radii corresponding to extrema of M_ap(θ₀) to calculate aperture mass ratios, such as θ₁ = 2 θ_E and θ₂ = θ_E. Then, M_ap(θ₁) will probe the annulus θ_E<θ< 2 θ_E, while for θ₂ ~ θ_E the annulus 0.5 θ_E ≤ θ ≤ θ_E would be probed. Figure 12 shows normalized aperture mass ratios M_ap(θ₁) /M_ap(θ₂) for the various SPT-transformed profiles studied in the previous section. We see that larger aperture mass ratios are found for larger values of f₂. In addition, the ratio depends only weakly of the shear amplitude γ_p, which means that the radial deformation of the mass profile produced by the SPT is mostly governed by the amplitude of f₂.

Fig. 11 Top: Map(θ0) (Eq. (34)) as a function of the “aperture” θ0. The filter function u(x) defined by Eq. (38), using x0 = 0.5, has been used such that for Map(θ0), the annulus [θ0/ (2θE),θ0/θE] is probed. The black curve shows Map for the NIS profile and the blue curves for the SPT-transformed profiles with γp = 0.05 and various values of f2. The green curve shows Map(θ0) for an MST transformed version of the NIS profile. Bottom: ratio between Map derived for the various transformed profiles and for the original NIS profile.

Fig. 12 Ratios of aperture mass between θ1 = 2 θE and θ2 = θE for SPT-transformed profiles with various values of f2. The ratios of aperture masses are normalized by the corresponding aperture ratios estimated for the original NIS profile (horizontal bar). The blue diamonds are for a shear γp = 0.05, and the red squares when γp = 0.1.

5.2. A specific lens model

Here, we apply the previous tests to three mass density profiles studied in SS13 and SS14, and used in those papers to illustrate degeneracies produced by the MST and the SPT. The reference model is a composite model constituted of the sum of a (spherically symmetric) Hernquist component to describe the baryonic component, and a (spherical) generalized Navarro-Frank-White (gNFW) density profile to describe the dark matter component of the galaxy. In addition, an external shear of amplitude γ_p = 0.1 is considered. Complex sets of lensed images from an ensemble of sources were generated with that model, and found to be all equally well reproduced by two (single) power-law profiles: a global power law, with an almost isothermal density slope γ′ = 1.98 (hereafter model M1), and a local power-law profile with a core radius θ_c = 0.1′′ and a slope γ′= 2.2 (hereafter model M2). We have applied the tests introduced in the previous subsection to these profiles to identify the nature of the degeneracy between the models. Figure 13 shows the difference R_κ between the original profile and transformed ones, i.e., ΔR_κ = R_κ(original)−R_κ(transformed). For comparison, we also show ΔR_κ obtained for an SPT with f₂ = 0.11. This figure suggests that indeed, this degeneracy is similar to an SPT.

The other diagnostic we present consists in calculating the aperture mass M_ap of the profiles. Figure 14 shows M_ap as a function of θ₀. In addition to the aperture mass calculated for the individual profiles, we also show the aperture mass corresponding to two different MST-transformed mass density profiles. As explained in SS13, model M1 is close to an MST transformed version of the composite model⁴ with λ = 0.84. On the other hand, Fig. 4 of SS14 shows that the MST contribution to M2 corresponds to λ = 0.932. Figure 14 is qualitatively similar to Fig. 11 but there is an offset of M_ap for M1 and M2 compared to the composite model. The reason is probably that the M1 and M2 profiles are transformed versions of the composite model via both an MST and an SPT. The MST contribution with λ = 0.84 is larger for M1 than for M2, for which λ ~ 0.93.

Fig. 13 Difference between Rκ calculated for three different pairs of profiles: In blue, an NIS and an SPT-transformed model with f2 = 0.11 and γp = 0.1; in magenta, a composite Hernquist+gNFW model and a power-law model (M1); in red, Hernquist+gNFW model and a cored power-law model (M2). The shape of ΔRκ for the models presented in SS13 and SS14 are qualitatively similar to that observed for the fiducial SPT model presented in Sect. 4.

Fig. 14 Top: Map(θ0) (Eq. (34)) as a function of θ0, for the composite Hernquist+gNFW model (black), for the power law model M1 (blue), and the cored power-law M2 (red). Dashed red (blue) profile shows Map(θ0) for an MST transformed version of the composite model with λ = 0.93 (resp. 0.84). Bottom: ratio of Map(θ0) between the “transformed” models and the composite. The dashed curves correspond to MST-transformed versions of the composite model, and represent the contribution of the MST to M1 and M2. The solid red and blue curves suggest that the remaining of the degeneracy can be associated with an SPT.

6. Discussion – Implications of the SPT for strong lensing

In this paper we have studied several aspects of the SPT, an invariance transformation of the deflection angle that leaves all multiple images properties of gravitational lenses invariant. The central question, of whether there exists a gravitational lensing potential which gives rise to a deflection angle sufficiently close to the SPT-transformed one (which in general is not curl free) has been explored for a particular class of lens models, namely an NIS with external shear and an SPT given as a radial stretching of the source plane. The radial stretching deformation was chosen such that the classical MST did not contribute in altering the original deflection since we are only interested in higher-order effects that go beyond the well known MST. This example has shown that, for a large range of parameters pairs (external shear and distortion parameter of the radial stretching) there indeed exist lensing potentials whose associated deflection is sufficiently close to the one obtained from the SPT that these two cannot be distinguish observationally. We conducted this study by formulating an action as the integral over the squared difference of these two deflection angles, yield a Neumann problem. We gave a detailed description of how this problem can be solved; these methods are expected to be useful for future theoretical studies and applications of the SPT.

We have considered only one criterion for the validity of an SPT, namely that the corresponding gradient deflection field does not deviate from the SPT-transformed deflection by more than 5 × 10^-3θ_E. Changing this observationally motivated limit to a different value will modify the space of allowed parameter combinations. For the example considered here, we expect that the allowed range of the stretching parameter f₂ for a given external shear will be proportional to the allowed maximum deviation of these two deflection angles.

We point out that our method of obtaining a gradient deflection law in form of the variational principle (6) does not necessarily yield an “optimal” modified deflection. As briefly discussed in Sect. 3.1, one could imagine alternative constraints for finding a gradiant deflection “close to ”. In particular, finding a gradient field whose maximum deviation from is minimized over the region of relevance would be a promising ansatz which, however, is analytically challenging, if at all doable. Nevertheless, the solution obtained in this paper yields valuable insight in the freedom of lens model choices offered by the SPT.

The properties of the mass distribution resulting from an SPT were also studied in detail. In contrast to the MST, which is a special case of an SPT, the more general SPT gives rise to non-monotonic changes in the radial mass profile, and to the generation of a finite ellipticity even if the original mass distribution was axi-symmetric. Hence, the SPT offers a much larger range of mass profile modifications which leave all strong lensing observables invariant, than does the MST. This more complex class of invariance transformation is of particular interest because it may be of great relevance when trying to fit real lens system (which are expected to have a rather complex mass profile; see, e.g., Xu et al. (2016) with simple lens models. The fact that simple mass models yield satisfactory fits even in cases with a rich observed image structure may be a consequence of some SPT which transforms the deflection of the true mass distribution into that of a simple mass model. Ignoring the potential complexity of the real mass distribution, and thus the possibility that the SPT may be acting, may lead to biases in estimates of physical parameters of the lens system.

We have defined several diagnostic quantities which can distinguish a general SPT from a pure MST. Applying these diagnostics to the special case of nearly degenerate lens models studies in two earlier papers, we conclude that this degeneracy can to a large degree be accounted by an MST, but that a non-negligible contribution comes from a more general SPT. Hence, an SPT has been found “empirically”, even before the concept of the SPT was developed. In that sense, the SPT is not just a “theoretical possibility” for obtaining different but observationally equivalent mass models, but describes degeneracies which actually occur in real lens modeling.

Acknowledgments

We would like to thank Bastian Orthen for valuable discussions, Olivier Wertz for valuable comments on this paper, and the referee, Matthias Bartelmann, for his constructive comments and advice. Part of this work was supported by the German Deutsche Forschungsgemeinschaft, DFG project numbers SL 172/1-1 and SCHN 342/13-1. Sandra Unruh is a member of the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne. Dominique Sluse is supported by a Back to Belgium grant from the Belgian Federal Science Policy (BELSPO).

References

Appendix A: Practical integration of Eq. (22) in a circular region

Calculating the deflection angle $\mathrm{\nabla }\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\psi }\end{array}$ in Eq. (18) by integrating the product over the circle poses a challenge, due to the pole of the first term in Eq. (21). To integrate over this pole, polar coordinates centered on the pole position ϑ need to be chosen. This can be done by a translation of the integration variable to x = θ−ϑ, and integrating in the polar coordinates of x. However, the integration range of the polar angle will depend on | x |, according to the geometrical overlap of circles centered on the origin and those centered on θ.

A better method is obtained by a conformal mapping of the form $\mathit{x}\mathrm{=}\frac{\mathit{\theta }\mathrm{-}\mathit{\vartheta }}{\mathit{R}\mathrm{-}{\mathit{\vartheta }}^{\mathrm{\ast }}\mathit{\theta }\mathit{/}\mathit{R}}\mathit{,}$(A.1)where we now use complex notation, i.e., x, ϑ and θ are complex numbers with components ϑ = ϑ₁ + iϑ₂ etc. and an asterisk denotes complex conjugation. This conformal mapping maps the circle | θ | <R onto the unit circle | x | < 1, and the singularity point θ = ϑ is mapped onto the origin x = 0. For example, setting θ = Re^iϕ, we get $\begin{array}{ccc}\mathit{x}\mathrm{=}\frac{\mathit{R}{\mathrm{e}}^{\mathrm{i}\mathit{\varphi }}\mathrm{-}\mathit{\vartheta }}{\mathit{R}\mathrm{-}{\mathit{\vartheta }}^{\mathrm{\ast }}{\mathrm{e}}^{\mathrm{i}\mathit{\varphi }}}\mathrm{=}{\mathrm{e}}^{\mathrm{i}\mathit{\varphi }}\frac{\mathit{R}\mathrm{-}\mathit{\vartheta }{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\varphi }}}{{\left(\mathit{R}\mathrm{-}\mathit{\vartheta }{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\varphi }}\right)}^{\mathrm{\ast }}}\mathit{,}& & \end{array}$from which it is immediately seen that | x | = 1. The inverse of the transformation (A.1) is readily obtained, $\mathit{\theta }\mathrm{=}\frac{\mathit{Rx}\mathrm{+}\mathit{\vartheta }}{\mathrm{1}\mathrm{+}{\mathit{\vartheta }}^{\mathrm{\ast }}\mathit{x}\mathit{/}\mathit{R}}\mathit{,}$(A.2)from which one can easily check that the unit circle | x | = 1 is mapped onto the circle | θ | = R. In components, Eq. (A.2) reads $\begin{array}{ccc}{\mathit{\theta }}_{\mathrm{1}}& \mathrm{=}& \frac{\mathit{R}{\mathit{x}}_{\mathrm{1}}\mathrm{+}{\mathit{\vartheta }}_{\mathrm{1}}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{|}\mathit{x}{\mathrm{|}}^{\mathrm{2}}\mathrm{)}\mathrm{+}\mathrm{[}{\mathit{x}}_{\mathrm{1}}\mathrm{(}{\mathit{\vartheta }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}{\mathit{\vartheta }}_{\mathrm{2}}^{\mathrm{2}}\mathrm{)}\mathrm{+}\mathrm{2}{\mathit{\vartheta }}_{\mathrm{1}}{\mathit{\vartheta }}_{\mathrm{2}}{\mathit{x}}_{\mathrm{2}}\mathrm{]}\mathit{/}\mathit{R}}{\mathrm{1}\mathrm{+}\mathrm{2}{\vartheta }\mathrm{·}{x}\mathit{/}\mathit{R}\mathrm{+}\mathrm{\left|}{\vartheta }{\mathrm{|}}^{\mathrm{2}}\mathrm{\right|}{x}{\mathrm{|}}^{\mathrm{2}}\mathit{/}{\mathit{R}}^{\mathrm{2}}}\mathit{,}\\ {\mathit{\theta }}_{\mathrm{2}}& \mathrm{=}& \frac{\mathit{R}{\mathit{x}}_{\mathrm{2}}\mathrm{+}{\mathit{\vartheta }}_{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{|}\mathit{x}{\mathrm{|}}^{\mathrm{2}}\mathrm{)}\mathrm{+}\mathrm{[}\mathrm{2}{\mathit{\vartheta }}_{\mathrm{1}}{\mathit{\vartheta }}_{\mathrm{2}}{\mathit{x}}_{\mathrm{1}}\mathrm{-}{\mathit{x}}_{\mathrm{2}}\mathrm{(}{\mathit{\vartheta }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}{\mathit{\vartheta }}_{\mathrm{2}}^{\mathrm{2}}\mathrm{)}\mathrm{]}\mathit{/}\mathit{R}}{\mathrm{1}\mathrm{+}\mathrm{2}{\vartheta }\mathrm{·}{x}\mathit{/}\mathit{R}\mathrm{+}\mathrm{\left|}{\vartheta }{\mathrm{|}}^{\mathrm{2}}\mathrm{\right|}{x}{\mathrm{|}}^{\mathrm{2}}\mathit{/}{\mathit{R}}^{\mathrm{2}}}\mathrm{·}\end{array}$(A.3)The Jacobi determinant of the transformation x → θ, needed for the integration, is $\left|\frac{\mathit{\partial }{\theta }}{\mathit{\partial }{x}}\right|\mathrm{=}\frac{{\mathit{R}}^{\mathrm{2}}\mathrm{(}{\mathit{R}}^{\mathrm{2}}\mathrm{-}\mathrm{|}{\vartheta }{\mathrm{|}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{(}{\mathit{R}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{R}{\vartheta }\mathrm{·}{x}\mathrm{+}\mathrm{|}{\vartheta }{\mathrm{|}}^{\mathrm{2}}\mathrm{|}{x}{\mathrm{|}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}$(A.4)which is non-zero for all x inside the unit circle and | ϑ | <R. As a sanity check, we calculate the area of the circle in the transformed coordinates, $\begin{array}{ccc}{\mathit{R}}^{\mathrm{2}}\mathit{\pi }& \mathrm{=}& \int \mathrm{𝒰}{\mathrm{d}}^{\mathrm{2}}\mathit{\theta }\mathrm{=}{\mathrm{\int }}_{\mathrm{𝒞}}{\mathrm{d}}^{\mathrm{2}}\mathit{x}\left|\frac{\mathit{\partial }{\theta }}{\mathit{\partial }{x}}\right|\\ & \mathrm{=}& {\mathit{R}}^{\mathrm{2}}\mathrm{(}{\mathit{R}}^{\mathrm{2}}\mathrm{-}\mathrm{|}{\vartheta }{\mathrm{|}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\mathrm{d}\mathit{xx}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\frac{\mathrm{d}\mathit{\varphi }}{\mathrm{(}{\mathit{R}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{R}{\vartheta }\mathrm{·}{x}\mathrm{+}\mathrm{|}{\vartheta }{\mathrm{|}}^{\mathrm{2}}\mathrm{|}{x}{\mathrm{|}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}\mathrm{·}\end{array}$(A.5)The inner integral yields 2π(R² + | x | ² | ϑ | ²)/(R²− | x | ² | ϑ | ²)³, the outer integral then gives π/ (R²− | ϑ | ²)², and we re-obtain the area πR².

In complex notation, the singular term in Eq. (21) reads $\begin{array}{ccc}\frac{\mathrm{1}}{\mathrm{(}\mathit{\vartheta }\mathrm{-}\mathit{\theta }{\mathrm{)}}^{\mathrm{\ast }}}\mathrm{=}\frac{\mathit{x}\mathrm{+}\mathrm{|}\mathit{x}{\mathrm{|}}^{\mathrm{2}}\mathit{\vartheta }\mathit{/}\mathit{R}}{\mathit{R}\mathrm{|}\mathit{x}{\mathrm{|}}^{\mathrm{2}}\left(\mathrm{|}\mathit{\vartheta }{\mathrm{|}}^{\mathrm{2}}\mathit{/}{\mathit{R}}^{\mathrm{2}}\mathrm{-}\mathrm{1}\right)}\mathit{,}& & \end{array}$yielding $\left|\frac{\mathit{\partial }{\theta }}{\mathit{\partial }{x}}\right|\frac{{\vartheta }\mathrm{-}{\theta }}{\mathrm{|}{\vartheta }\mathrm{-}{\theta }{\mathrm{|}}^{\mathrm{2}}}\mathrm{=}\frac{{\mathit{R}}^{\mathrm{2}}\mathrm{\left(}\mathrm{\right|}{\vartheta }{\mathrm{|}}^{\mathrm{2}}\mathrm{-}{\mathit{R}}^{\mathrm{2}}\mathrm{)}}{{\left({\mathit{R}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{R}{\vartheta }\mathrm{·}{x}\mathrm{+}\mathrm{\left|}{\vartheta }{\mathrm{|}}^{\mathrm{2}}\mathrm{\right|}{x}{\mathrm{|}}^{\mathrm{2}}\right)}^{\mathrm{2}}}\left(\frac{\mathit{R}}{{\left|{x}\right|}^{\mathrm{2}}}{x}\mathrm{+}{\vartheta }\right)\mathit{.}$(A.6)We can check the consistency of this expression by calculating the deflection angle of a uniform disk with surface mass density κ₀, which reads $\begin{array}{ccc}{\alpha }\mathrm{\left(}{\vartheta }\mathrm{\right)}& \mathrm{=}& \frac{{\mathit{\kappa }}_{\mathrm{0}}}{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{𝒰}}{\mathrm{d}}^{\mathrm{2}}\mathit{\theta }\frac{{\vartheta }\mathrm{-}{\theta }}{\mathrm{|}{\vartheta }\mathrm{-}{\theta }{\mathrm{|}}^{\mathrm{2}}}\mathrm{=}\frac{{\mathit{\kappa }}_{\mathrm{0}}}{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{𝒞}}{\mathrm{d}}^{\mathrm{2}}\mathit{x}\left|\frac{\mathit{\partial }{\theta }}{\mathit{\partial }{x}}\right|\frac{{\vartheta }\mathrm{-}{\theta }}{\mathrm{|}{\vartheta }\mathrm{-}{\theta }{\mathrm{|}}^{\mathrm{2}}}\\ & \mathrm{=}& \frac{{\mathit{R}}^{\mathrm{2}}\mathrm{(}\mathrm{|}{\vartheta }{\mathrm{|}}^{\mathrm{2}}\mathrm{-}{\mathit{R}}^{\mathrm{2}}\mathrm{)}{\mathit{\kappa }}_{\mathrm{0}}}{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\mathrm{d}\mathit{xx}\\ & & \mathrm{\times }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\frac{\mathrm{d}\mathit{\varphi }}{{\left({\mathit{R}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{R}{\vartheta }\mathrm{·}{x}\mathrm{+}\mathrm{|}{\vartheta }{\mathrm{|}}^{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}\right)}^{\mathrm{2}}}\left(\frac{\mathit{R}}{{\mathit{x}}^{\mathrm{2}}}{x}\mathrm{+}{\vartheta }\right)\\ & \mathrm{=}& \frac{{\mathit{R}}^{\mathrm{2}}\mathrm{(}\mathrm{|}{\vartheta }{\mathrm{|}}^{\mathrm{2}}\mathrm{-}{\mathit{R}}^{\mathrm{2}}\mathrm{)}{\mathit{\kappa }}_{\mathrm{0}}}{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\mathrm{d}\mathit{xx}\hspace{0.17em}\left(\frac{\mathrm{-}\mathrm{2}\mathit{\pi }}{{\left({\mathit{R}}^{\mathrm{2}}\mathrm{-}\mathrm{|}{\vartheta }{\mathrm{|}}^{\mathrm{2}}{\mathit{x}}^{\mathrm{2}}\right)}^{\mathrm{2}}}\right){\vartheta }\\ & \mathrm{=}& \frac{{\mathit{R}}^{\mathrm{2}}\mathrm{(}\mathrm{|}{\vartheta }{\mathrm{|}}^{\mathrm{2}}\mathrm{-}{\mathit{R}}^{\mathrm{2}}\mathrm{)}{\mathit{\kappa }}_{\mathrm{0}}}{\mathit{\pi }}\left(\frac{\mathit{\pi }}{{\mathit{R}}^{\mathrm{2}}\left(\mathrm{|}{\vartheta }{\mathrm{|}}^{\mathrm{2}}\mathrm{-}{\mathit{R}}^{\mathrm{2}}\right)}\right){\vartheta }\mathrm{=}{\mathit{\kappa }}_{\mathrm{0}}\hspace{0.17em}{\vartheta }\mathit{,}\end{array}$(A.7)as expected.

Appendix B: Proof of the relation (26) for a circular region

In this section we use again the complex notation for two-dimensional vectors, in terms of which the vector field (21) reads ${\mathrm{\nabla }}_{\mathit{\vartheta }}\mathit{H}\mathrm{\left(}{\vartheta }\mathrm{;}{\theta }\mathrm{\right)}\mathrm{\to }\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\left(\frac{\mathrm{1}}{{\mathit{\vartheta }}^{\mathrm{\ast }}\mathrm{-}{\mathit{\theta }}^{\mathrm{\ast }}}\mathrm{-}\frac{\mathit{\theta }}{\mathrm{(}{\mathit{R}}^{\mathrm{2}}\mathrm{-}{\mathit{\vartheta }}^{\mathrm{\ast }}\mathit{\theta }\mathrm{)}}\mathrm{-}\frac{\mathit{\vartheta }}{{\mathit{R}}^{\mathrm{2}}}\right)\mathrm{·}$(B.1)Therefore, we obtain for the fields A and B_in the complex expressions $\begin{array}{ccc}\mathit{A}\mathrm{\left(}\mathit{\vartheta }\mathrm{\right)}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{𝒰}}{\mathrm{d}}^{\mathrm{2}}\mathit{\theta \kappa }\mathrm{\left(}{\theta }\mathrm{\right)}\hspace{0.17em}\left(\frac{\mathit{\theta }}{\mathrm{(}{\mathit{R}}^{\mathrm{2}}\mathrm{-}{\mathit{\vartheta }}^{\mathrm{\ast }}\mathit{\theta }\mathrm{)}}\mathrm{+}\frac{\mathit{\vartheta }}{{\mathit{R}}^{\mathrm{2}}}\right)\mathit{,}\\ {\mathit{B}}_{\mathrm{in}}\mathrm{\left(}\mathit{\vartheta }\mathrm{\right)}& \mathrm{=}& \end{array}$(B.2)where a_in = α_in·n, defined on the boundary of . We first calculate this product, noting that for a boundary point R e^iϕ of the circle, n = e^iϕ. Specializing the complex expression for α_in as given in Eq. (23) to a point on the boundary, ϑ = Rn, we obtain $\begin{array}{ccc}{\mathit{a}}_{\mathrm{in}}\mathrm{\left(}\mathit{Rn}\mathrm{\right)}& \mathrm{=}& {{\alpha }}_{\mathrm{in}}\mathrm{·}\mathit{n}\mathrm{=}\frac{{\mathit{\alpha }}_{\mathrm{in}}\mathit{/}\mathit{n}\mathrm{+}{\mathit{\alpha }}_{\mathrm{in}}^{\mathrm{\ast }}\mathit{n}}{\mathrm{2}}\\ & \mathrm{=}& \end{array}$(B.3)We now insert this expression into Eq. (B.2). The integral over the boundary is written as ds = R dϕ = −i R dn/n, and the integral over n extends over the unit circle. With θ = Rn, we then find $\begin{array}{ccc}{\mathit{B}}_{\mathrm{in}}\mathrm{\left(}\mathit{\vartheta }\mathrm{\right)}& \mathrm{=}& \frac{\mathit{R}}{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}_{\mathrm{𝒰}}{\mathrm{d}}^{\mathrm{2}}\mathit{\theta \kappa }\mathrm{\left(}{\theta }\mathrm{\right)}\hspace{0.17em}\frac{\mathrm{-}\mathrm{i}}{\mathrm{2}\mathit{\pi }}\mathrm{􏽉}\frac{\mathrm{d}\mathit{n}}{\mathit{n}}\left(\frac{\mathrm{1}}{\mathit{R}\mathrm{-}\mathit{n}{\mathit{\theta }}^{\mathrm{\ast }}}\mathrm{+}\frac{\mathit{n}}{\mathit{Rn}\mathrm{-}\mathit{\theta }}\right)\\ & & \end{array}$(B.4)The inner integrand is an analytic function of n inside the unit circle, except at the poles at n = 0 and at n = θ/R (ϑ and θ are both inside the circle). Applying the theorem of residue, the integral can thus be evaluated. The first pole yields a contribution −ϑ/R³, whereas the second pole results in the expression $\begin{array}{ccc}\frac{\mathit{\theta }}{\mathrm{(}\mathit{\theta }{\mathit{\vartheta }}^{\mathrm{\ast }}\mathrm{-}{\mathit{R}}^{\mathrm{2}}\mathrm{)}\mathit{R}}\mathrm{-}\frac{\mathit{\theta }}{\mathrm{(}{\mathit{R}}^{\mathrm{2}}\mathrm{-}\mathit{\theta }{\mathit{\vartheta }}^{\mathrm{\ast }}\mathrm{)}\mathit{R}}\mathrm{-}\frac{\mathit{\vartheta }}{{\mathit{R}}^{\mathrm{3}}}\mathrm{·}& & \end{array}$Adding up these two contributions then yields ${\mathit{B}}_{\mathrm{in}}\mathrm{\left(}\mathit{\vartheta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{𝒰}}{\mathrm{d}}^{\mathrm{2}}\mathit{\theta \kappa }\mathrm{\left(}{\theta }\mathrm{\right)}\left(\mathrm{-}\frac{\mathit{\theta }}{\mathrm{(}{\mathit{R}}^{\mathrm{2}}\mathrm{-}{\mathit{\vartheta }}^{\mathrm{\ast }}\mathit{\theta }\mathrm{)}}\mathrm{-}\frac{\mathit{\vartheta }}{{\mathit{R}}^{\mathrm{2}}}\right)\mathit{,}$(B.5)which we see agrees with the expression for A in Eq. (B.2). Thus we have shown explicitly that for the kernel function (20), the relation (26) holds.