© ESO 2018

1 Introduction

The use of the gravitational lensing phenomenon as a cosmological tool offers an independent way to probe the nature of the universe (for the early work, see Blandford & Narayan 1992). To date, numerous weak and strong lensing observations have been employed to infer the fundamental cosmological parameters with an increasingly competitive precision. In the strong lensing regime, Refsdal (1964) established that multiple-image systems can theoretically be used to infer the Hubble parameter H₀. The method relies upon the idea that the propagation time of light rays emitted from a background source (typically an active galactic nucleus, AGN) towards the observer differs from one lensed image to another. The corresponding difference in arrival times, known as the time delay, is inversely proportional to H₀. This idea lays the basis of the modern time delay cosmography, which has been extensively addressed in literature; see, for example, the recent review Treu & Marshall (2016) and references therein.

For sake of clarity, we recall few key results of the well-known theory of time delays. Relative to an unperturbed ray emitted by a source located at β, the extra light travel time T(θ) at an image position θ is formally defined by $$T(\theta )=\frac{D_{\text{Δ}t}}{c}\left\{\frac{1}{2}{\left[\theta -\beta (\theta )\right]}^2-\psi (\theta )\right\}\text{ =:}\frac{D_{\text{Δ}t}}{c}\tau (\theta ),$$(1)

where ψ(θ) is the deflection potential produced by a dimensionless surface mass density κ(θ) = ∇²ψ(θ)∕2, τ(θ) is known as the Fermat potential, and D_Δt is referred to as the time delay distance $$D_{\text{Δ}t}=(1+z_{\text{d}})\frac{D_{\text{d}}{D}_{\text{s}}}{D_{\text{ds}}}\propto H^{-1}{}_0,$$(2)

where z_d is the redshift of the deflector and D is the angular diameter distances between the observer and deflector (d), observer and source (s), and deflector and source (ds). In Eq. (1), the first term in brackets describes the geometrical deviation of the light ray due to the lens, whereas the second describes the time delay that a ray experiences as it crosses the deflection potential. The relative time delay Δ t_ij between a pair of lensed images θ_i and θ_j is obtained by differencing the corresponding extra light travel time: $$\text{Δ}{t}_{ij}=T({\theta }_i)-T({\theta }_j)=\frac{D_{\text{Δ}t}}{c}\left[\tau ({\theta }_i)-\tau ({\theta }_j)\right]\text{ =:}\frac{D_{\text{Δ}t}}{c}\text{Δ}{\tau }_{ij}.$$(3)

From Eq. (3), H₀ inference can be conceptually performed by constraining the time delay distance D_Δt, provided that both accurate time delay measurements and a mass model which predicts Δ τ_ij can be obtained. At present, a few percent precision time delays have been measured for several multiple-image systems based on different light curve analysis methods (see e.g., Vuissoz et al. 2008; Paraficz & Hjorth 2010; Courbin et al. 2011; Fohlmeister et al. 2013; Eulaers et al. 2013; Tewes et al. 2013; Rathna Kumar et al. 2013; Bonvin et al. 2017; Akhunov et al. 2017). In the foreseeable future, we can expect thousands of lensed quasars to be discovered by the next generation of instruments (Jean et al. 2001; Coe & Moustakas 2009; Oguri & Marshall 2010; Finet et al. 2012, 2015; Liao et al. 2015; Finet & Surdej 2016). Among them, numerous suitable candidates for robust time delay measurements should lead the time delay cosmography to the next level. However, constraining the lens mass distribution turns out to be as decisive as measuring time delays with high precision. Given a measured time delay between two lensed images, more concentrated mass distributions lead to shorter time delay distance estimations, hence to larger values of H₀ (Kochanek 2002). The Fermat potential difference Δτ_ij is primarily sensitive to the strong lensing effects produced by the main lens. However, a realistic time delay cosmography should also consider the lensing effects of any external mass structures located in the vicinity of the main lens, as well as along the line of sight (LOS; e.g., Seljak 1994; Bar-Kana 1996). If the LOS mass effects are sufficiently small, they can be approximated by an external shear and an external convergence, usually denoted as κ_ext, which need to be characterized¹ (see e.g., Keeton 2003; Fassnacht et al. 2006; Suyu et al. 2010, 2013; Wong et al. 2011, 2017). Otherwise, these external mass structures need to be explicitly included in the mass model, for instance by considering the full multi-plane lensing formalism (Schneider 2014b; McCully et al. 2014, 2016).

As firstshown in Falco et al. (1985), the dimensionless surface mass density κ(θ) and the class of mass models κ_λ(θ) defined as $${\kappa }_{\lambda }(\theta )=\lambda \kappa (\theta )+(1-\lambda ),$$(4)

along with the corresponding unobservable source rescaling β → λ β, lead to identical lensing observables, except for the time delays between pairs of lensed images which are transformed such that Δ t → λ Δt. If not broken, this degeneracy, referred to as the mass-sheet transformation (MST), may significantly affect cosmographic inferences, including H₀ (see e.g., Gorenstein et al. 1988; Saha 2000; Wucknitz 2002; Koopmans et al. 2003; Liesenborgs & De Rijcke 2012; Schneider & Sluse 2013; Schneider 2014a,b; Xu et al. 2016). It is worth mentioning that the external convergence κ_ext is based on physical effects whereas the MST (Eq. (4)) stems from a pure mathematical degeneracy (Schneider & Sluse 2013; hereafter SS13). Different solutions have been proposed to reduce the degeneracy induced by the MST in time delay cosmography (see e.g., Birrer et al. 2016; Treu & Marshall 2016, and references therein). A commonly used method consists in assuming a specific lens model, typically a power-law, and independently estimating the lens mass with the measurement of its velocity dispersion. However, SS13 have shown experimentally that two different classes of galaxy models with compatible velocity dispersions were able to reproduce equally well a set of image positions, but predicted significantly different time delays. Furthermore, because the time delay ratios were not constant, they suggested that the transformation between these two models was not exactly an MST but a more general one. This has naturally raised some concerns about the reliability of the H₀ determination from time delay cosmography.

Schneider & Sluse (2014; hereafter SS14) laid the theoretical basis for an approximate invariance transformation, the so-called source-position transformation (SPT), of which the MST is a special case. Unruh et al. (2017; hereafter USS17) explored further its properties, such as defining a criterion to determine whether an SPT is valid or not and exploring the density profile of SPT-modified mass distributions. They also pointed out that the degeneracy found experimentally in SS13 between the two models mimics an SPT, which thereby confirmed that it occurs in real lens modeling. To date, it is not clear whether the conclusions drawn in SS13 and SS14 about time delays and H₀ could be generalized to other mass distributions modified under exact SPTs or only reflect the behavior of a very special case. In this paper, we address this question by studying how time delays are sensitive to the effects of the SPT.

This paper is organized as follows. For readers who are not familiar with the SPT, we outline its basic principles in Sect. 2. In particular, we recall the importance of identifying a validity criterion. Owing to the valuable insight it offers for more general cases, we consider the SPT-modified time delays for axisymmetric lenses in Sect. 3. For nonaxisymmetric lenses, we discuss in Sect. 4 the relevance of the validity criterion defined in USS17 and analyze the SPT-modified time delays in detail. We summarize our findings and conclude in Sect. 5.

2 The principle of the source position transformation

This section focuses on the principle of the SPT and the most recent theoretical results. For a detailed discussion, the reader is referred to SS14 and USS17. All the analytical results presented in this paper have been implemented into a user friendly python package called pySPT. All the numerical results and figures have also been obtained from pySPT without using any extra software. For an overall description of the package, we refer the reader to the companion paper Wertz & Orthen (2018).

The basic idea underlying the SPT can be simply summarized as follows. A given general mass distribution κ(θ) defines a deflection law α(θ) which describes how the light paths are affected in the vicinity of the deflector. The n lensed image angular positions θ_i of a point-like source at unobservable position β are those which satisfy the lens equation β = θ_i −α(θ_i). Then, from astrometric observations, we can infer the constraints $${\theta }_i-\alpha ({\theta }_i)={\theta }_j-\alpha ({\theta }_j),$$(5)

for all 1 ≤ i < j ≤ n, leading to the mapping θ_i(θ₁) defined by the relative image positions of the same source. The SPT addresses the following question: can we define an alternative deflection law, denoted as $\widehat{\alpha }(\theta )$, which preserves the mapping θ_i(θ₁) for a uniquesource? If such a deflection law exists, it will necessarily be associated with the alternative source position $\widehat{\beta }={\theta }_i-\widehat{\alpha }({\theta }_i)$, defining a new lens mapping, in such a way that $$\theta =\beta +\alpha (\theta )=\widehat{\beta }+\widehat{\alpha }(\theta ).$$(6)

An SPT consists in a global transformation of the source plane formally defined by a mapping $\widehat{\beta }(\beta )$ which gives rise to the transformed deflection law: $$\widehat{\alpha }(\theta )=\alpha (\theta )+\beta -\widehat{\beta }(\beta )=\alpha (\theta )+\beta -\widehat{\beta }(\theta -\alpha (\theta )),$$(7)

where in the first step we used Eq. (6) and in the last step we inserted the original lens equation. The mapping $\widehat{\beta }(\beta )$ is chosen so that it satisfies $\text{det}(\partial \widehat{\beta }/\partial \beta )\ne 0$ for all β in the region of interest, hence $\widehat{\beta }(\beta )$ is one-to-one. This property of the source mapping guarantees the pairing of images to be conserved. With $\widehat{\alpha }$ defined thisway, Eq. (6) guarantees that all images of a given source β under the original deflection law α(θ) are also images of the source $\widehat{\beta }$ under the modified deflection law $\widehat{\alpha }(\theta )$. Therefore, the mapping θ_i(θ₁) is preserved for all source positions.

From the Jacobi matrix $\widehat{\mathcal{A}}(\theta )=\partial \widehat{\beta }/\partial \theta =(\partial \widehat{\beta }/\partial \beta )(\partial \beta /\partial \theta )$ of the modified lens mapping $\widehat{\beta }=\theta -\widehat{\alpha }(\theta )$, SS14 have shown that both the magnification ratios of image pairs and their relative shapes remain unchanged under an SPT. In general, the Jacobi matrix $\widehat{\mathcal{A}}(\theta )$ is not symmetric, which indicates that the modified deflection law $\widehat{\alpha }$ is not a curl-free field: $$|\nabla \times \widehat{\alpha }(\theta )|=\left|{\widehat{\mathcal{A}}}_{12}(\theta )-{\widehat{\mathcal{A}}}_{21}(\theta )\right|\ne 0,$$(8)

where the subscript indices refer to the matrix entries. Therefore, $\widehat{\alpha }$ cannot be in general expressed as the gradient of a deflection potential $\widehat{\psi }$ and does not correspond to the deflection produced by a gravitational lens. Thus, there exists no physical mass distribution $\widehat{\kappa }$ leading to the modified deflection law $\widehat{\alpha }$. The only cases for which $\widehat{\mathcal{A}}(\theta )$ is globally symmetric occur either when the SPT simply reduces to an MST $\widehat{\beta }(\beta )=\lambda \beta$, or when axisymmetric lenses are transformed under SPTs corresponding to a general radial stretching of the form $$\widehat{\beta }(\beta )=\left[1+f(|\beta |)\right]\beta ,$$(9)

where f is called the deformation function. For such cases, we can always define $\widehat{\kappa }$ so that $2\widehat{\kappa }=\nabla \cdot \widehat{\alpha }={\nabla }^2\widehat{\psi }$. However, even in this case there is still no guarantee that $\widehat{\kappa }$ corresponds to a physical mass distribution. Depending on the SPT, the modified mass profile may become nonmonotic or even non-positive definite in particular regions of the lens plane.

Provided the curl component of $\widehat{\alpha }$ is sufficiently small, it was shown in USS17 that one can define a curl-free deflection law $\tilde{\alpha }$ which is very similar to $\widehat{\alpha }$ in the sense that their difference is smaller than the astrometric accuracy ε_acc of current observations: $$|\tilde{\alpha }(\theta )-\widehat{\alpha }(\theta )|\text{ =:}|\text{Δ}\alpha (\theta )|<{\epsilon }_{\text{acc}},$$(10)

in a finite region $\mathcal{U}$ where multiple images occur. Therefore, $\tilde{\alpha }$ can be derived as the gradient of a deflection potential $\tilde{\psi }$, which is caused by a mass distribution $\tilde{\kappa }$ corresponding to a gravitational lens. The central question of the validity of an SPT was addressed in USS17. Whereas $\widehat{\alpha }$ yields exactly the same lensed image positions as the original lens, $\tilde{\alpha }$ does not. Because of observational uncertainties and additional physical reasons such as substructures in the mass distribution, we cannot reproduce observed positions to better than a few milliarcseconds (mas) with a smooth mass model (for a detailed discussion, see SS14). A given SPT $\widehat{\beta }(\beta )$ should be flagged as being valid as long as a corresponding curl-free $\tilde{\alpha }$ leads to lensed image shifts smaller than a few mas. In this sense, the SPT is only an approximate invariance transformation. The condition (10) was chosen in USS17 as the criterion to assess whether an SPT is valid or not. The relevance of this choice is reconsidered in detail in Sect. 4.1.

Because it will be of practical interest for deriving SPT-modified time delays in the nonaxisymmetric case (see Sect. 4), we recall here the explicit expressions for $\tilde{\psi }$ and $\tilde{\alpha }=\nabla \tilde{\psi }$. These can essentially be obtained by formulating the “action” $$S={\displaystyle {\int }_{\mathcal{U}}{\left|\nabla \tilde{\psi }-\widehat{\alpha }\right|}^2}{\text{d}}^2\theta$$(11)

for which finding a minimum leads to the Neumann problem: $$\{\begin{array}{cccc}& {\nabla }^2\tilde{\psi }=\nabla \cdot \widehat{\alpha }\text{ =:}2\widehat{\kappa },& \text{for all\hspace{0.17em}}& \theta \in \mathcal{U},\\ & \nabla \tilde{\psi }\cdot n=\widehat{\alpha }\cdot n,& \text{for all\hspace{0.17em}}& \theta \in \partial \mathcal{U},\end{array}$$(12)

where $\partial \mathcal{U}$ represents the boundary curve of $\mathcal{U}$ and n the outward directed normal vector. From Eq. (11), we see that the condition (10) plays a central role in defining a curl-free counterpart $\tilde{\alpha }$ of the SPT-modified deflection law $\widehat{\alpha }$. We also note that the first relation in Eq. (12) implies $\tilde{\kappa }=\widehat{\kappa }$ for all $\theta \in \mathcal{U}$. The Neumann problem can be solved by means of a Green’s function for which an analytical solution is known when $\mathcal{U}$ is a disk of radius R. Thus, the deflection potential $\tilde{\psi }$ evaluated at the position ϑ in the lens plane explicitly reads (Unruh et al. 2017; Wertz & Orthen 2018) $$\tilde{\psi }(\vartheta )=\langle \tilde{\psi }\rangle +2{\displaystyle {\int }_{\mathcal{U}}{H}_1}(\vartheta ;\theta )\widehat{\kappa }(\theta ){\text{d}}^2\theta -{\displaystyle {\int }_{\partial \mathcal{U}}{H}_2}(\vartheta ;\theta )\widehat{\alpha }\cdot n\text{\hspace{0.17em}d}s,$$(13)

where $\langle \tilde{\psi }\rangle$ is the average of $\tilde{\psi }$ on $\mathcal{U}$, ds the line element of the boundary curve $\partial \mathcal{U}$, $$H_1(\vartheta ;\theta )=\frac{1}{4\pi }\left[\ln \left(\frac{{\left|\vartheta -\theta \right|}^2}{R^2}\right)+\ln \left(1-\frac{2\vartheta \cdot \theta }{R^2}+\frac{\left|\vartheta {|}^2\right|\theta {|}^2}{R^4}\right)-\frac{|\theta {|}^2}{R^2}\right],$$(14)

and $$H_2(\vartheta ;\theta )=\frac{1}{4\pi }\left[2\ln \left(\frac{{\left|\vartheta -\theta \right|}^2}{R^2}\right)-1\right].$$(15)

The corresponding deflection angle $\tilde{\alpha }$ can be derived by obtaining the gradient of H₁ and H₂ with respect to ϑ, which reads $$\begin{array}{lll}\tilde{\alpha }(\vartheta )\hfill & =\hfill & \frac{1}{\pi }{\displaystyle {\int }_{\mathcal{U}}\left(\frac{\vartheta -\theta }{|\vartheta -\theta {|}^2}+\frac{|\theta {|}^2\vartheta -R^2\theta }{R^4-2R^2\vartheta \cdot \theta +\left|\vartheta {|}^2\right|\theta {|}^2}\right)}\widehat{\kappa }(\theta )\text{\hspace{0.17em}d}\theta \hfill \\ \hfill & \hfill & -\frac{1}{\pi }{\displaystyle {\int }_{\partial \mathcal{U}}\frac{\vartheta -\theta }{|\vartheta -\theta {|}^2}}\widehat{\alpha }\cdot n\text{\hspace{0.17em}d}s.\hfill \end{array}$$(16)

Of course, $\tilde{\psi }$ and $\tilde{\alpha }$ depend on the radius R of the circular region $\mathcal{U}$ and must satisfies the two conditions R > |ϑ| and R not too large to ensure the criterion (10) to be satisfied.

To quantitatively study the impact of the SPT on time delays, it will be necessary to explicitly define a mapping $\widehat{\beta }(\beta )$. We will focus most of this work on an isotropic SPT described by the radial stretching of the form (9). In particular, we will consider the special case where the deformation function f(|β|) is the lowest-order expansion of more general functions: $$f(|\beta |)=f_0+\frac{f_2}{2{\theta }_{\text{E}}^2}|\beta {|}^2,$$(17)

where f₀f(0), := $f_2\text{ :=}{\theta }_{\text{E}}^2{f}^{″}(0)$, and θ_E is the Einstein angular radius. When f₂ = 0, Eq. (17) reduces to f₀ and the radial stretching (9) simplifies to a pure MST with λ = 1 + f₀. Such as in SS14 and USS17, we only consider SPT parameters which yield to physically meaningful modified mass profiles.

3 Time delays: the axisymmetric case

Owing to its simplicity, the study of how an SPT affects time delays between lensed images produced by an axisymmetric lens provides a valuable insight on the general nonaxisymmetric case. Since α and θ are collinear, the original lens mapping becomes one-dimensional and reads β = θ − α(θ). We set β > 0 and only consider the two outer² lensed images θ_A and θ_B located on opposite sides of the lens center, that is θ_B < 0 < |θ_B| < θ_A. From Eq. (1), we readily deduce the one-dimensional form of the original time delay Δ t_AB between the image pair (θ_A, θ_B): $$\text{Δ}{t}_{\text{AB}}=\frac{D_{\text{Δ}t}}{c}[\tau ({\theta }_{\text{A}})-\tau ({\theta }_{\text{B}})]\text{ =:}\frac{D_{\text{Δ}t}}{c}\text{Δ}{\tau }_{\text{AB}}.$$(18)

The one-dimensional radial stretching (9) simply reads $$\widehat{\beta }(\beta )=[1+f(\beta )]\beta ,$$(19)

where f(−β) = f(β) to preserve the symmetry. With no loss of generality, 1 + f(β) + β df(β)∕dβ > 0 ensures the SPT to be one-to-one. For the axisymmetric case, the SPT is an exact invariance transformation. Thus, the deflection law $\widehat{\alpha }$ is a curl-free field, $\tilde{\alpha }=\widehat{\alpha }$, and there exists a deflection potential $\widehat{\psi }$ such as $$\frac{\text{d}\widehat{\psi }(\theta )}{\text{d}\theta }=\widehat{\alpha }(\theta )=\alpha (\theta )-f(\beta (\theta ))\beta (\theta ),$$(20)

where in the second step we used the one-dimensional form of Eq. (7). From Eq. (1), we deduce that the SPT-modified extra light travel time $\widehat{T}$ reads $$\widehat{T}(\theta )=\frac{D_{\text{Δ}t}}{c}\left[\frac{1}{2}{\left(\theta -\widehat{\beta }[\beta (\theta )]\right)}^2-\widehat{\psi }(\theta )\right]\text{ =:}\frac{D_{\text{Δ}t}}{c}\widehat{\tau }(\theta ).$$(21)

From Eqs. (3) and (21), the SPT-modified time delay between image pair (θ_A, θ_B) of the same source thus becomes $$\text{Δ}{\widehat{t}}_{\text{AB}}=\widehat{T}({\theta }_{\text{A}})-\widehat{T}({\theta }_{\text{B}})=\frac{D_{\text{Δ}t}}{c}(\widehat{\tau }({\theta }_{\text{A}})-\widehat{\tau }({\theta }_{\text{B}}))\text{ =:}\frac{D_{\text{Δ}t}}{c}\text{Δ}{\widehat{\tau }}_{\text{AB}}.$$(22)

With Eqs. (18) and (22), we show in Sect. 3.1 that the time delay ratios $\text{Δ}{\widehat{t}}_{\text{AB}}/\text{Δ}{t}_{\text{AB}}(\equiv \text{Δ}{\widehat{\tau }}_{\text{AB}}/\text{Δ}{\tau }_{\text{AB}})$ can be highly simplified, revealing an elegant expression in terms of β and $\widehat{\beta }(\beta )$, and valid for any axisymmetric lens and deformation function f(β). We also propose an equivalent form of this relation in terms of the original and SPT-modified mean surface mass densities. We illustrate the analytical results with some examples in Sect. 3.2.

3.1 The SPT-modified time delays

After substituting the one-dimensional form of Eqs. (6) and (20) into Eq. (21), the SPT-modified extra light travel time reads $\widehat{T}({\theta }_i)=\widehat{\tau }({\theta }_i)D_{\text{Δ}t}/c$ with $$\widehat{\tau }({\theta }_i)=\frac{1}{2}{[\alpha ({\theta }_i)-f(\beta ({\theta }_i))\beta ({\theta }_i)]}^2-\psi ({\theta }_i)+{\displaystyle {\int }_0^{{\theta }_i}f}(\beta (\theta ))\beta (\theta )\text{\hspace{0.17em}d}\theta ,$$(23)

up to a constant independent of θ, keeping in mind that β(θ) = θ − α(θ). Because of β(θ_A) = β(θ_B) =: β_s, we have f(β(θ_A)) = f(β(θ_B)) = f(β_s), and the SPT-modified time delays between the images θ_A and θ_B is given by $\text{Δ}{\widehat{t}}_{\text{AB}}=\text{Δ}{\widehat{\tau }}_{\text{AB}}{D}_{\text{Δ}t}/c$ with $$\text{Δ}{\widehat{\tau }}_{\text{AB}}=\text{Δ}{\tau }_{\text{AB}}-f({\beta }_{\text{s}}){\beta }_{\text{s}}({\theta }_{\text{A}}-{\theta }_{\text{B}})+{\displaystyle {\int }_{\left|{\theta }_{\text{B}}\right|}^{{\theta }_{\text{A}}}f}(\beta (\theta ))\beta (\theta )\text{\hspace{0.17em}d}\theta .$$(24)

Due to the lens symmetry, the integral over [θ_B, |θ_B|] does not contribute to $\text{Δ}{\widehat{\tau }}_{\text{AB}}$. With no loss of generality, we thus integrate from |θ_B| instead of θ_B in Eq. (24). To go a step further, the difference Δτ_AB between the original Fermat potentials can also be written as $$\begin{array}{lll}\text{Δ}{\tau }_{\text{AB}}\hfill & =\hfill & -{\beta }_{\text{s}}({\theta }_{\text{A}}-{\theta }_{\text{B}})+\frac{1}{2}\left({\theta }_{\text{A}}^2-{\theta }_{\text{B}}^2\right)-(\psi ({\theta }_{\text{A}})-\psi ({\theta }_{\text{B}})),\hfill \\ \hfill & =\hfill & -{\beta }_{\text{s}}({\theta }_{\text{A}}-{\theta }_{\text{B}})+{\displaystyle {\int }_{\left|{\theta }_{\text{B}}\right|}^{{\theta }_{\text{A}}}\beta }(\theta )\text{\hspace{0.17em}d}\theta ,\hfill \end{array}$$

where in the first step we used the original lens equations α(θ_A) = θ_A − β_s and α(θ_B) = θ_B − β_s, and in the last step we used dψ(θ)∕dθ = θ − β(θ). Combining Eqs. (24) and (26), we then obtain from Eq. (22) the SPT-modified time delay: $$\text{Δ}{\widehat{t}}_{\text{AB}}=\text{Δ}{t}_{\text{AB}}\left[1+f({\beta }_{\text{s}})\right]+\frac{D_{\text{Δ}t}}{c}{\epsilon }_{\text{AB}},$$(27)

where we define ε_AB as $${\epsilon }_{\text{AB}}={\displaystyle {\int }_{\left|{\theta }_{\text{B}}\right|}^{{\theta }_{\text{A}}}\beta }(\theta )\left[f(\beta (\theta ))-f({\beta }_{\text{s}})\right]\text{\hspace{0.17em}d}\theta .$$(28)

For the special case of a pure MST, the deformation function f is independent of β, namely f(β(θ)) = f(β_s) ≡ λ − 1 with $\lambda \in ℝ$. Therefore, ε_AB = 0 and we find $\text{Δ}{\widehat{t}}_{\text{AB}}=\lambda \text{Δ}{t}_{\text{AB}}$ for all axisymmetric lenses, as expected. Considering the radial stretching (19) and a singular isothermal sphere (SIS) lens model, we show explicitly in Appendix A that ε_AB = 0 also holds for all image pairs (θ_A, θ_B), that is for 0 ≤ β < θ_E. In fact, simple analytical arguments reveal that, in general, ε_AB remains very small compared to the other terms in Eq. (27) and can be neglected. The demonstration is explained in detail in Appendix B. As a result, the time delay ratios $\text{Δ}\widehat{t}/\text{Δ}t$ given in Eq. (27) can be simply approximated by $$\frac{\text{Δ}\widehat{t}}{\text{Δ}t}\approx 1+f(\beta )\equiv \frac{\widehat{\beta }(\beta )}{\beta },$$(29)

where we have dropped the subscript AB keeping in mind that the equation holds only for time delay ratios between the same pair of lensed images corresponding to the source $\widehat{\beta }$ and β. For a given radial stretching, Eq. (29) shows that the ratios between SPT-modified and original time delays scale basically like $\widehat{\beta }/\beta$, implying that they depend explicitly on the deformation function f(β), as it is the case for the MST.

As written, Eq. (29) misleadingly suggests that the time delay ratio is insensitive to the original lens profile κ. Consider two original radial mass profiles κ⁽¹⁾ and κ⁽²⁾, which are not related under an SPT, and consider a source position β. We locate the corresponding pairs of brighter lensed images by $({\theta }_{\text{A}}^{(1)},{\theta }_{\text{B}}^{(1)})$ and $({\theta }_{\text{A}}^{(2)},{\theta }_{\text{B}}^{(2)})$. For a given deformation function f(β), Eq. (29) says that $\text{Δ}{\widehat{t}}_{\text{AB}}^{(1)}/\text{Δ}{t}_{\text{AB}}^{(1)}\approx \text{Δ}{\widehat{t}}_{\text{AB}}^{(2)}/\text{Δ}{t}_{\text{AB}}^{(2)}$, but the two time delay ratios are evaluated at two different pairs of positions which depend on the lens models, that is ${\theta }_{\text{A}}^{(1)}\ne {\theta }_{\text{A}}^{(2)}$ and ${\theta }_{\text{B}}^{(1)}\ne {\theta }_{\text{B}}^{(2)}$. When κ⁽²⁾ corresponds to a modified version of κ⁽¹⁾ under the SPT $\widehat{\beta }(\beta )=[1+g(\beta )]\beta$ (with g(β) satisfying the conditions given after Eq. (19)), we have ${\theta }_{\text{A}}^{(1)}={\theta }_{\text{A}}^{(2)}$ and ${\theta }_{\text{B}}^{(1)}={\theta }_{\text{B}}^{(2)}$. However, this case can be reduced to an original radial mass profile κ⁽¹⁾ deformed by an SPT that is defined as the composition of two other SPTs such as $\widehat{\beta }(\beta )=[1+h(\beta )]\beta$ with h(β) = [1 + f(β)][1 + g(β)] − 1. Thus, this leads to $\text{Δ}\widehat{t}/\text{Δ}t\approx 1+h(\beta )$, in agreement with Eq. (29).

The SPT-modified mass profile $\widehat{\kappa }$ of a radial profile κ is also radial (SS14). Therefore, time delays Δt_AB and $\text{Δ}{\widehat{t}}_{\text{AB}}$ should depend only on the image positions and the corresponding surface mass densities in the annulus defined between the images. In particular for Δ t_AB, the major contribution comes from the mean surface mass density ${\langle \kappa \rangle }_{\text{AB}}$ in the annulus $\left|{\theta }_{\text{B}}\right|<\theta <{\theta }_{\text{A}}$ (Gorenstein et al. 1988; Kochanek 2002, 2006). We will show next that the time delay ratios (29) can be expressed only in terms of ${\langle \kappa \rangle }_{\text{AB}}$ and the corresponding SPT-modified ${\langle \widehat{\kappa }\rangle }_{\text{AB}}$. First, we can easily show that $$\begin{array}{lll}{\langle \kappa \rangle }_{\text{AB}}\hfill & \text{:=}\hfill & \frac{2}{{\theta }_{\text{A}}^2-{\theta }_{\text{B}}^2}{\displaystyle {\int }_{\left|{\theta }_{\text{B}}\right|}^{{\theta }_{\text{A}}}\theta }\kappa (\theta )\text{\hspace{0.17em}d}\theta =\frac{m({\theta }_{\text{A}})-m(|{\theta }_{\text{B}}|)}{{\theta }_{\text{A}}^2-{\theta }_{\text{B}}^2}\hfill \\ \hfill & =\hfill & 1-\frac{{\beta }_{\text{s}}}{{\theta }_{\text{A}}-\left|{\theta }_{\text{B}}\right|},\hfill \end{array}$$(30)

where in the last step we used m(θ) = θ α(θ), α(θ_A) = θ_A − β_s, and α(θ_B) = θ_B − β_s. Similarly, we can easily deduce that $${\langle \widehat{\kappa }\rangle }_{\text{AB}}=1-\frac{\widehat{\beta }({\beta }_{\text{s}})}{{\theta }_{\text{A}}-\left|{\theta }_{\text{B}}\right|},$$(31)

where we first used $\widehat{m}(\theta )=\theta \widehat{\alpha }(\theta )$ and Eq. (20), then α(θ_A) = θ_A − β_s and α(θ_B) = θ_B − β_s. Combining Eqs. (29) and (31), we thus obtain for the time delay ratio in terms of mean surface mass densities: $$\frac{\text{Δ}\widehat{t}}{\text{Δ}t}\approx \frac{\widehat{\beta }(\beta )}{\beta }=\frac{1-\langle \widehat{\kappa }\rangle }{1-\langle \kappa \rangle },$$(32)

where we have once again dropped the AB keeping in mind that the mean surface mass densities are evaluated in the annulus defined by the inner and outer radii |θ_B| and θ_A, respectively. As expected, Eq. (32) shows that the ratio between SPT-modified and original time delays depends essentially on mean surface mass densities in the annulus $\left|{\theta }_{\text{B}}\right|<\theta <{\theta }_{\text{A}}$. Finally, we note that the second equality in Eq. (32) is exact.

3.2 Some illustrative examples

To illustrate the results obtained in the previous section, we first consider the deformation function (17) with f₂ = 0.5 and f₀ = 0 to separate theimpact of the MST from that of the SPT. As original lens model, we choose a non-singular isothermal sphere (NIS) characterized by the deflection law: $$\alpha (\theta )=\frac{{\theta }_{\text{E}}\theta }{\sqrt{{\theta }_{\text{c}}^2+{\theta }^2}},$$(33)

where the core radius θ_c is defined such as θ_c = ν θ_E with 0 < ν < 1. To derive time delays in the axisymmetric case, we only consider the three lensed image configurations where the fainter central image is omitted. Thus, we need to sample the source positions inside the radial caustic of angular radius β_r = β(θ_r) where θ_r represents the angular radius of the corresponding radial critical curve. For an NIS, β_r is simply given by ${\beta }_{\text{r}}={\theta }_{\text{E}}{(1-{\nu }^{2/3})}^{3/2}\approx 0.695{\theta }_{\text{E}}$ for ν = 0.1. Using this simple lens model and pySPT, we create pairs of mock images for a uniform set of 34 sources covering the range β = 0.02 θ_E to β = 0.68 θ_E < β_r.

The top panel in Fig. 1 shows $\text{Δ}\widehat{t}/\text{Δ}t$ as a function of β for the corresponding pairs of lensed images. We see that the time delay ratios scale remarkably well like the function $1+f(\beta )\equiv \widehat{\beta }(\beta )/\beta$, as predicted by Eq. (29). According to Eq. (29), the strongest effect of the SPT on time delays arises for a source as close as possible to the radial caustic, that is for β → β_r. Thus, in our first example, the theoretical maximum time delay ratio equals $\text{Δ}\widehat{t}/\text{Δ}t\approx 1.12$ for f₂ = 0.5 and ν = 0.1, leading to an impact of 12% on H₀. As long as it leads to a physical meaningful $\widehat{\kappa }$, larger (resp. smaller) values of |f₂| lead to larger (resp. smaller) time delay ratios. To quantitatively evaluate the accuracy of Eq. (29), we compare the numerically evaluated quantity $1-(\text{Δ}\widehat{t}/\text{Δ}t)/(\widehat{\beta }/\beta )$ to unity, as shown in the bottom panel in Fig. 1. From Eq. (27), it follows that $1-(\text{Δ}\widehat{t}/\text{Δ}t)/(\widehat{\beta }/\beta )=\left|{\eta }_{\text{AB}}\right|$ with ${\eta }_{\text{AB}}\text{ :=}({\epsilon }_{\text{AB}}/\text{Δ}{\tau }_{\text{AB}})/(\widehat{\beta }/{\beta }_{\text{s}})$. The quantity |η_AB| is smaller than 10⁻⁴ for β ≤ 0.5 θ_E, reaching a maximum of |η_AB|≈ 6 × 10⁻⁴ ≪ 1 for β = 0.68 θ_E, confirming that ε_AB can be neglected in Eq. (29) in such a case. For an NIS deformed by a radial stretching characterized by Eq. (17), it is possible to derive an analytical solution for ε_AB, hence for |η_AB|, by solving Eq. (28). This analytical solution is represented in the bottom panel in Fig. 1 and fits perfectly the numerical evaluations of $1-(\text{Δ}\widehat{t}/\text{Δ}t)/(\widehat{\beta }/\beta )$ at each source position, as expected.

We have successfully tested the relation (29) for various axisymmetric lens profiles deformed by different deformation functions. As additional examples, we consider the two deformation functions: $$f(\beta )=\frac{2f_0}{\cosh (\beta /{\beta }_0)-f_0},$$(34)

with ${\beta }_0={\theta }_{\text{E}}\sqrt{3(1-f_0)/(1+f_0)}$ and f₀ = −0.32, and $$f(\beta )=f_0+{\beta }_0^2{f}_2{\beta }^2\left[2\left({\beta }_0^2+{\beta }^2\right)\right],$$(35)

with β₀ = 0.8 θ_E, f₀ = −1∕3, and f₂ = 1∕9. The choice for the two deformation functions (34) and (35) is justified by the fact that the resulting SPT-modified mass profiles $\widehat{\kappa }$ are approximately power laws near the tangential critical curve, that is $\widehat{\kappa }(\theta )\approx \widehat{\kappa }({\theta }_{\text{E}}){(\theta /{\theta }_{\text{E}})}^{-\upsilon }$ (SS14). In both cases, we adopt an NIS with θ_c = 0.1 θ_E as original lens model and the same source sample as in the first example. Figure 2 shows the time delay ratios between the lensed images for each source. As expected, $\text{Δ}\widehat{t}/\text{Δ}t$ fits almost perfectly the function $1+f(\beta )\equiv \widehat{\beta }(\beta )/\beta$. For β = 0, the two deformation functions simplify to f(β) = f₀ and the corresponding SPTs reduce to pure MSTs leading to $\text{Δ}\widehat{t}/\text{Δ}t=0.68$ and $\text{Δ}\widehat{t}/\text{Δ}t=2/3$, respectively. Therefore, any changes from these values reflect the impact of the SPT. For β = 0.68 θ_E, Fig. 2 shows an impact of around 3.6% and 2.2% on H₀, respectively, which is significantly smaller than what we have obtained for the first example.

Not all combinations of SPT deformation parameters and original mass profiles κ yield a physically meaningful SPT-modified mass profile, namely $\widehat{\kappa }$ monotonically decreasing and positive definite (Schneider & Sluse 2014). In addition, the maximum time delay ratio also depends on κ since the latter directly defines the size of the radial caustic β(θ_r), namely the region in the source plane that produces multiple images. In summary, the way the SPT affects the time delays is very sensitive to the choice of the deformation function f, the associated deformation parameters, the original mass profile κ and lensed image positions. For these reasons, we restrain ourselves to draw generalized quantitative conclusions in the axisymmetric case. However, our numerical tests suggest an effect of a few percent in general. We will show in the next section that the simple connection between the time delay ratios and the source position ratios may still be very strong in the nonaxisymmetric case.

Fig. 1 Impact of an SPT described by the radial stretching $\widehat{\beta }(\beta )=1+f_2{\beta }^2/(2{\theta }_{\text{E}}^2)$, with f2= 0.5, on time delays image pairs generated by a NIS, with core θc = 0.1 θE. Top panel: ratio between SPT-modified and original time delays (black inverted triangles) for each source position. The time delay ratios scale almost perfectly like $\widehat{\beta }/\beta$ (red curve), as predicted by Eq. (29). Bottom panel: numerical confirmation of the validity of the approximation performed in Eq. (29) which consists in neglecting the term εAB. The solid black curve illustrates |ηAB| as an analytical function of β and perfectly fits the quantity $1-(\text{Δ}\widehat{t}/\text{Δ}t)/(\widehat{\beta }/\beta )$ numerically evaluated for each source position.

Fig. 2 Impact on time delays of two different SPTs defined such that the corresponding $\widehat{\kappa }$ is approximately a power law near the tangential critical curve. As predicted by Eq. (29), the time delay ratios (black dots and squares) scale almost perfectly like $\widehat{\beta }/\beta$ (red and green curves).

4 Time delays: the nonaxisymmetric case

In this section, we drop the axisymmetry assumption for the original lens model. The SPT-modified deflection law $\widehat{\alpha }$ is thus not a curl-free field in general and there exists no deflection potential $\widehat{\psi }$ which satisfies $\nabla \widehat{\psi }=\widehat{\alpha }$. To define a physically meaningful modified extra light travel time, we consider the deflection law $\tilde{\alpha }$, the closest curl-free approximation to $\widehat{\alpha }$ which satisfies the criterion (10) for all θ over a region $\mathcal{U}$ (see Eq. (16)), and the associated deflection potential $\tilde{\psi }$ (see Eq. (13)). For the rest of this section, we follow USS17 and consider ε_acc ≈ 5 × 10⁻³ θ_E over the circular region |θ|≤ 2 θ_E, where the approximation stems from the typical positional accuracy of the Hubble Space Telescope.

Within the region $\mathcal{U}$, the lensed images $\tilde{\theta }$ of the source $\widehat{\beta }$ satisfying the lens mapping $\widehat{\beta }=\tilde{\theta }-\tilde{\alpha }(\tilde{\theta })$ are expected to be sufficiently close to the corresponding original images $\theta =\beta +\alpha (\theta )=\widehat{\beta }+\widehat{\alpha }(\theta )$ to not be distinguished observationally. However, we show in Sect. 4.1 that the criterion (10) defined in USS17 cannot guarantee the difference $\left|\text{Δ}\theta \right|\text{ :=}|\tilde{\theta }-\theta |$ between the SPT-modified image position $\tilde{\theta }$ of the source $\widehat{\beta }$ and the image position θ of the source β to be smaller thanε_acc over the whole region $\mathcal{U}$. However, for specific pairs of original and SPT-modified sources leading to indistinguishable image configurations, we illustrate in Sect. 4.2 the typical behavior of the time delay ratios. Finally, based solely on a numerical optimization, we slightly modify the source mapping $\widehat{\beta }(\beta )$ by relaxing the isotropic condition of the SPT. It follows that the region where |Δ θ | < ε_acc can be substantially extended. From this ad hoc source mapping, $\tilde{\alpha }$ and $\tilde{\psi }$, we illustrate the corresponding alternative time delay ratios in Sect. 4.3.

4.1 Criterion for the validity of an SPT

To illustrate the limit of the criterion |Δα(θ)| < ε_acc, we consider a situation similar to SS14 and USS17, namely a quadrupole lens composed of an NIS plus external shear γ_p (NISg) for which the deflection law is defined by $$\alpha (\theta )=\frac{{\theta }_{\text{E}}\theta }{\sqrt{{\theta }_{\text{c}}^2+|\theta {|}^2}}-{\gamma }_{\text{p}}\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)\theta ,$$(36)

where the core radius is set to θ_c = 0.1 θ_E. The originalsource mapping is transformed by a radial stretching (9) with a deformation function f of the form (17). The adopted SPT is thus defined by $$\widehat{\beta }(\beta )=\left(1+f_0+\frac{f_2}{2{\theta }_{\text{E}}^2}|\beta {|}^2\right)\beta .$$(37)

For the rest of this section, we set the deformation parameter f₂ to be f₂ = 0.4 and the external shear magnitude to be γ_p = 0.1. Furthermore, we exclude the effect of a pure MST by simply choosing f₀ = 0. According to USS17 (see their Fig. 4), this specific pair (f₂, γ_p) constitutes an allowed pair of parameters in a sense it fulfills the criterion |Δ α (θ)| < ε_acc over the circular region |θ | ≤ 2 θ_E. Figure 3 shows the map |Δα(θ)| over a circular grid |θ | ≤ 2 θ_E in the lens plane. This figure is similar to the map |Δα(θ)| illustrated in Fig. 7 in USS17, although they used f₂ = 0.55. It turns out that the ratio |Δα(θ_i)|∕|Δα(θ_j)| remains unaffected when f₂ varies, but is sensitive to variations of R or γ_p. Actually, γ_p is the only parameter that explicitly characterizes the degree of asymmetry of the original lens model.

The next step consists in determining how well the deflection law $\tilde{\alpha }$ allows us to reproduce the original lensed images of a source. To this aim, we create a set of mock images θ of a sample of sources β that cover the first quadrant of the source plane. We restrict the grid of sources to 0 ≤ |β|≤ 0.66 θ_E where multiple images are produced (see top panel in Fig. 4). Then, we produce the images $\tilde{\theta }$ of the corresponding SPT-modified sources $\widehat{\beta }(\beta )$. The bottom panel in Fig. 4 shows the image positions θ and the color-coding represents |Δθ| in units of θ_E. The same color-coding is applied to the sources (top panel) where only the largest offset, denoted by |Δ θ |_max, are shown.

An “unexpected” conclusion can be drawn from Fig. 4; even though |Δ α (θ)| < ε_acc over the region |θ|≤ 2 θ_E (as shown in Fig. 3), many of image configurations are characterized by |Δ θ |≫ ε_acc for at least one lensed image. This implies that these image configurations can be distinguished from the original ones and the corresponding SPT can no longer be flagged as valid. Furthermore, the largest offsets |Δ θ | occur near the tangential critical curve. It comes with no surprise that the corresponding regions in the source plane are thus located near the tangential caustic curve. To address this behavior, we first consider what the quantity |Δ α (θ)| really represents. As defined in Eq. (10), both $\widehat{\alpha }$ and $\tilde{\alpha }$ are evaluated at the same position θ in the lens plane. Therefore, we have $\widehat{\alpha }(\theta )=\theta -\widehat{\beta }$ and $\tilde{\alpha }(\theta )=\theta -\tilde{\beta }$ where $\tilde{\beta }$ is the source position of the image θ under the deflection law $\tilde{\alpha }$. Combining the two latter equations leads to $$\text{Δ}\alpha (\theta )=\widehat{\beta }(\theta )-\tilde{\beta }(\theta )\text{ =:}\text{Δ}\beta .$$(38)

Equation (38) shows that minimizing |Δα(θ)| is equivalent to minimizing |Δβ| with no guarantee on |Δθ|. Indeed, let us consider a position θ close to a critical line for which |Δα(θ)| < ε_acc, for example (θ_x∕θ_E, θ_y∕θ_E) = (0.5, 1.0) (see Fig. 3 and bottom panel in Fig. 4). The corresponding source β (θ) is necessarily close to a caustic, so is $\widehat{\beta }(\theta )$. Thus, the source $\tilde{\beta }(\theta )$ lies in a region of the source plane where even small shifts |Δβ| can lead to significantly different image positions. This explains why regions where |Δ θ |≫ ε_acc are those which surround the critical curves. Furthermore, whereas ${\theta }_i-\widehat{\alpha }({\theta }_i)={\theta }_j-\widehat{\alpha }({\theta }_j)$ is satisfied for all i ≤ j, we have ${\theta }_i-\tilde{\alpha }({\theta }_i)\ne {\theta }_j-\tilde{\alpha }({\theta }_j)$, meaning that the θ_i are not lensed images of a unique source under the deflection law $\tilde{\alpha }$. Thus, the criterion |Δα(θ_i)| < ε_acc for a lensed image configuration θ_i is based upon positions that are not linked under the deflection law $\tilde{\alpha }$. These few simple arguments suggest with no loss of generality that the choice of Eq. (10) as a validity criterion may not be the most appropriate one.

Let us now evaluate $\tilde{\alpha }$ at the position $\tilde{\theta }=\widehat{\beta }+\tilde{\alpha }(\tilde{\theta })$ instead of θ and consider the difference $\tilde{\alpha }(\tilde{\theta })-\widehat{\alpha }(\theta )$. We readily find that $$\left|\tilde{\alpha }\left(\tilde{\theta }\right)-\widehat{\alpha }(\theta )\right|=\left|\tilde{\theta }-\theta \right|=\left|\text{Δ}\theta \right|,$$(39)

which corresponds exactly to the image shift induced by the SPT that we expect to be smaller than ε_acc. Assuming that $\left|\text{Δ}\theta \right|$ is small, we can show to first order that $$\begin{array}{lll}\widehat{\beta }\hfill & =\hfill & \tilde{\theta }-\tilde{\alpha }(\tilde{\theta }),\hfill \\ \hfill & =\hfill & \theta +\text{Δ}\theta -\widehat{\alpha }(\theta +\text{Δ}\theta )-\text{Δ}\alpha (\theta +\text{Δ}\theta ),\hfill \\ \hfill & \approx \hfill & \theta -\widehat{\alpha }(\theta )+\left(1-\frac{\partial \widehat{\alpha }}{\partial \theta }\right)\text{Δ}\theta -\text{Δ}\alpha (\theta ).\hfill \end{array}$$(40)

Thus, for all positions |θ|≤ 2 θ_E not located on a critical curve, Eq. (40) leads to $$\left|\text{Δ}\theta \right|\approx \left|{\widehat{\mathcal{A}}}^{-1}(\theta )\text{Δ}\alpha (\theta )\right|.$$(41)

Equation (41) clearly shows that the offsets |Δθ| are related to |Δα(θ)| through the SPT-modified Jacobi matrix $\widehat{\mathcal{A}}(\theta )$ whose impact become larger as we get closer to the critical curves. Figure 5 illustrates the quantity $\left|{\widehat{\mathcal{A}}}^{-1}(\theta )\text{Δ}\alpha (\theta )\right|$ using two different color-coding. The upper panel shows the same color-coding as used in Fig. 3 for comparison. It shows that a significant part of the region|θ| < 2 θ_E is characterized by |Δθ| > ε_acc even though|Δα(θ)| < ε_acc. The lower panel adopts a color-coding which allows us to better visualize regions characterized by the largest offsets. These regions surround the two critical curves represented by the two green lines. We confirm the validity of the first order Eq. (41) by comparing Fig. 5 with Fig. 6, which represents explicitly the quantity $\left|\tilde{\theta }-\theta \right|$. As expected, small differences can be observed very close to the critical curves where higher order terms in Eq. (40) become significant and cannot be ignored. In addition, Fig. 6 is much more time consuming to obtain than Fig. 5. Indeed, a single $\left|\tilde{\theta }-\theta \right|$ evaluation requires $\tilde{\theta }$ to be calculated first, that is solving the lens equation $\widehat{\beta }=\tilde{\theta }-\tilde{\alpha }(\tilde{\theta })$ that implies numerous $\tilde{\alpha }$ evaluations. In contrast, a single Eq. (41) evaluation requires only one $\tilde{\alpha }$ evaluation. For this reason, the grid density in bottom panel in Fig. 5 is 20 times higher than in Fig. 6.

Equation (41) confirms that the criterion |Δ α (θ)| < ε_acc for the validity of an SPT cannot guarantee the angular separation |Δθ| to be smaller than the astrometric accuracy of current observations, at least in regions nearby critical curves. To construct the curl-free deflection field $\tilde{\alpha }$, USS17 have considered the action defined in Eq. (11) for which they found a minimum. This approach is based explicitly on the validity criterion (10), which is not satisfactory and should be reconsidered. A new appropriate criterion would of course imply the definition of a new action to be minimized, leading to a new definition for $\tilde{\alpha }$ and $\tilde{\psi }$. Such a new approach is beyond the scope of this paper and will not be addressed here. Nevertheless, it remains possible to quantitatively estimate the impact of the SPT on time delays with the means available. In the next section, we first focus on the subset of source positions β depicted in the top panel in Fig. 4 that yields |Δθ| < ε_acc.

Fig. 3 Map of |Δα(θ)| over the circular region |θ|≤ 2 θE for f2 = 0.4, θc = 0.1 θE, and γp = 0.1. This figure issimilar to Fig. 7 in Unruh et al. (2017) with f2 = 0.55, even thoughit is based on a different approach (see the text for more details).

Fig. 4 Top panel: grid of source positions β covering the radial range 0 ≤|β|≤ 0.66 θE in the first quadrant of the source plane. Each source produces a set of lensed images θ (shown in the bottom panel) under the original deflection law α and the corresponding set of $\tilde{\theta }$ under the curl-free deflection law $\tilde{\alpha }$. The color-coding refers to the largest offset |Δθ|max associated with each source. The solid black curves locate the caustic curves (top panel) and the critical curves (bottom panel), respectively, for the NIS plus external shear. Bottom panel: set of mock lensed images θ produced by the source positions shown in the top panel and lensed by an NIS plus external shear characterized by θc = 0.1 θE and γp = 0.1. The color-coding represents the offsets |Δθ| in units of θE between θ and the images $\tilde{\theta }$.

Fig. 5 Maps of $\left|{\widehat{\mathcal{A}}}^{-1}(\theta )\text{Δ}\alpha (\theta )\right|\approx \left|\text{Δ}\theta \right|$ in units of θE for f2 = 0.4 and γp = 0.1. Top panel: we use the same color-coding as in Fig. 3 to explicitly show that |Δ α (θ)| < 5 × 10−3 θE over a region of the lens plane does not guarantee the image offsets |Δθ| to be smallerthan 5 × 10−3 θE over the same region. Bottom panel: we adjust the color-coding to bring out regions where the image offsets |Δ θ | are the largest, namely the critical curves.

Fig. 6 Map of $|\tilde{\theta }-\theta |\text{ =:}\left|\text{Δ}\theta \right|$ in units of θE for f2 = 0.4 and γp = 0.1. Clear differences with Fig. 5 are observed for positions θ located almost on the radial critical curve, where the approximation adopted in Eq. (41) is expected to become not valid. This map is much more time consuming to obtain than the ones represented in Fig. 5.

4.2 The SPT-modified time delays for valid configurations

In the previous section, we have shown that the criterion defined in Eq. (10) does not guarantee |Δ θ | < ε_acc for all |θ | ≤ 2 θ_E. However, the top panel in Fig. 4 also shows sources (mainly outside the tangential caustic curve) for which the corresponding largest offsets |Δθ|_max between original and SPT-modified image configurations are smaller than ε_acc. Adopting the same original lens model as in the previous section, Fig. 7 shows the quantity $$|\text{Δ}\theta (\beta ){|}_{\text{max}}\text{ := max}\left(\left|{\widehat{\mathcal{A}}}^{-1}(\theta (\beta ))\text{Δ}\alpha (\theta (\beta ))\right|\right),$$(42)

in units of θ_E over the first quadrant in the source plane. The region outside the radial caustic curve is irrelevant in our case because it does not lead to multiple image configurations. The two disjointed hatched regions denoted by B₁ and B₂ are defined such as all sources inside B₁ ∪B₂ lead to image configurations characterized by |Δθ|_max ≤ ε_acc. For f₂ = 0.4 and γ_p = 0.1, the region B₁ ∪B₂ covers around 56% of the area enclosed by the radial caustic curve. Smaller values for f₂ yield larger B₁ ∪B₂ areas, up to 100% when f₂ = 0 (SPT reduced to an MST) or γ_p = 0 (axisymmetric lens). The regions B₁ and B₂ are situated on both sides of the tangential caustic curve. The very high area ratio between these two regions (1–560 in this case) indicates that most of the valid image configurations are composed of two images (the fainter third central one is always omitted). Moreover, the few valid four component configurations are very symmetric, suggesting comparable time delays between opposite image pairs.

Provided that $\widehat{\kappa }$ is physically meaningful, the curl-free deflection field $\tilde{\alpha }$ yields indistinguishable image configurations for sources β ∈B₁ ∪B₂ as compared to the original α. Although these valid image configurations are of limited interest for time delay cosmography³, the resulting model ambiguities may still prevent us from performing a robust lens modeling. Thus, even though the adopted SPT is not valid over all the region |θ|≤ 2 θ_E, we propose in this section to analyze the time delay ratios of these particular image configurations between the original and SPT-modified models. To this aim, we consider an original nonaxisymmetric mass distribution which produces n lensed images θ_i of a source β ∈B₁ ∪B₂. The time delay Δ t_ij between a pair of lensed images θ_i and θ_j is defined in Eq. (3). The corresponding SPT-modified time delay $\text{Δ}{\tilde{t}}_{ij}$ have to be evaluated at image positions ${\tilde{\theta }}_i$ and ${\tilde{\theta }}_j$, respectively, leading to $$\text{Δ}{\tilde{t}}_{ij}=\tilde{T}\left({\tilde{\theta }}_i\right)-\tilde{T}\left({\tilde{\theta }}_j\right)=\frac{D_{\text{Δ}t}}{c}\left[\tilde{\tau }\left({\tilde{\theta }}_i\right)-\tilde{\tau }\left({\tilde{\theta }}_j\right)\right]\text{ =:}\frac{D_{\text{Δ}t}}{c}\text{Δ}{\tilde{\tau }}_{ij},$$(43)

where the SPT-modified Fermat potential is defined by $$\tilde{\tau }\left(\tilde{\theta }\right)=\frac{1}{2}{\left[\tilde{\theta }-\widehat{\beta }\left(\tilde{\theta }\right)\right]}^2-\tilde{\psi }\left(\tilde{\theta }\right)=\frac{1}{2}{\left|\tilde{\alpha }\left(\tilde{\theta }\right)\right|}^2-\tilde{\psi }\left(\tilde{\theta }\right).$$(44)

We present here the representative results obtained for two classes of models: the quadrupole NISg as defined in Eq. (36) and a non-singular isothermal elliptical lens (NIE). The NIE surface mass density κ is intrinsically nonaxisymmetric and is defined by (see e.g., Keeton 2001) $$\kappa (\theta )=\frac{{\theta }_{\text{E}}}{2\sqrt{{\theta }_{\text{c}}^2+{\rho }^2}},$$(45)

where the variable ρ, constant on ellipses with axis ratio $q=\sqrt{(1-ϵ)/(1+ϵ)}$, is characterized by $$\rho =\sqrt{\frac{{\theta }_x^2}{1-ϵ}+\frac{{\theta }_y^2}{1+ϵ}}.$$(46)

For the rest of this section, the quadrupole model parameters are set to (θ_c, γ_p) = (0.1 θ_E, 0.1) and the NIE model parameters to (θ_c, ϵ) = (0.1 θ_E, 0.15). These models are deformed by an SPT corresponding to a radial stretching defined in Eq. (37) with f₀ = 0 and f₂ = 0.4. In both cases, we used pySPT to create mock images θ of two separated grids of sources β, which cover together the corresponding regions B₁ ∪B₂. The size and shape of B₁ ∪B₂ are defined by both the model and SPT parameters, hence differ from the quadrupole to the NIE (see Fig. 7 for the quadrupole and Fig. 8 for the NIE). Making use of Eq. (16), we compute the corresponding images $\tilde{\theta }$ of the SPT-modifed sources $\widehat{\beta }$. We finally derive the time delay Δt_ij and $\text{Δ}{\tilde{t}}_{ij}$ and represent their ratios in Fig. 9 (for the NISg) and Fig. 10 (for the NIE). Because of their similarities, we discuss the NISg and the NIE simultaneously.

In the double image regime, the time delay ratios scale almost like the ratios $\left|\widehat{\beta }\right|/\left|\beta \right|$. The color-coding in Fig. 10 refers to the azimuth angle of β, denoted as ∠β. Even though the dispersion is slightly larger for the NIE, the deviations from $\left|\widehat{\beta }\right|/\left|\beta \right|$ do not exceed 0.5% in all cases. The inset in Fig. 10 clearly shows that the dispersion of the time delay ratios is the effect of the relative direction of β with respect to the orientation of the axis of the elliptical iso-density contours (here equal to 0°). The deviations from $\left|\widehat{\beta }\right|/\left|\beta \right|$ are maximum for ∠β = 0° and ∠β = 90°, and minimum for ∠β ≈ 45°. A similar behavior is observed for the NISg, but with respect to the orientation of the external shear (also equal to 0°). We suggest that Eq. (29), valid for the axisymmetric case (see Sect. 3.1), may also be valid in the nonaxisymmetric case for sufficiently large values of |β|, $$\frac{\text{Δ}\tilde{t}}{\text{Δ}t}\approx 1+f(|\beta |)=\frac{\left|\widehat{\beta }(\beta )\right|}{\left|\beta \right|}.$$(47)

Actually, even for the most unfavorable cases, Eq. (47) provides at least a fairly good estimate of $\text{Δ}\tilde{t}/\text{Δ}t$. It turns out that these two particular examples are representative of the time delay ratio behavior for double image configurations produced by a nonaxisymmetric lens. Thus, the impact of the SPT in the double image regime comes mainly from the ratios $\left|\widehat{\beta }\right|/\left|\beta \right|$, in the same way as for the axisymmetric case. In particular, the largest time delay ratios ${(\text{Δ}\tilde{t}/\text{Δ}t)}_{\text{max}}$ (considering only the double image configurations for now) is obtained for the source position β ∈B₂ characterized by the largest radial coordinate |β| and denoted as βmax. Therefore, ${(\text{Δ}\tilde{t}/\text{Δ}t)}_{\text{max}}$ depends on βmax and the latter depends on both the deformation function and the original lens model parameters, which define the size of B₂. For the NISg model depicted in Fig. 9, we find ${\beta }_{\text{max}}^{\text{NISg}}\approx (0.562,0.413){\theta }_{\text{E}}$, $\left|{\beta }_{\text{max}}^{\text{NISg}}\right|\approx 0.697{\theta }_{\text{E}}$, leading to ${(\text{Δ}\tilde{t}/\text{Δ}t)}_{\text{max}}^{\text{NISg}}\approx 1.121$, that is an impact of around 12% on H₀. For the NIE model depicted in Fig. 10, we find ${\beta }_{\text{max}}^{\text{NIE}}\approx (0.335,0.0){\theta }_{\text{E}}$, $\left|{\beta }_{\text{max}}^{\text{NIE}}\right|\approx 0.335{\theta }_{\text{E}}$, leading to ${(\text{Δ}\tilde{t}/\text{Δ}t)}_{\text{max}}^{\text{NIE}}\approx 1.028$, that is an impact of around 3% on H₀. Similar to the axisymmetric case, the impact of the SPT on time delays may substantially vary according to the nature of the original lensmodel.

A different behavior is observed for the case of the quadruple image regime. As first pointed out in SS13 from a pure empirical case, the time delay ratios of image pairs between the original and SPT-modified models are not conserved, that is ${(\text{Δ}\tilde{t}/\text{Δ}t)}_{ij}\ne {(\text{Δ}\tilde{t}/\text{Δ}t)}_{ik}$ with i < k ≤ 4. For this reason, even though only three independent time delays can be obtained from a quadruple image configurations, we represent in Figs. 9 and 10 the time delay ratios for all six image permutations (i, j) ∈ [(1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4)]. We note that the criterion chosen for ordering the images (1–4) is the extra light travel time, from smallest to largest. The pair of opposite images, namely (θ₁, θ₂) and (θ₃, θ₄), leads to $\text{Δ}\tilde{t}/\text{Δ}t$ close to 1, regardless of the adopted original lens model we have tested (see the pairs of green triangles close to 1 in bottom panels in Figs. 9 and 10). Owing to the symmetry of the image configurations (B₁∋ β ~ 0), the time delays Δt₁₂ and Δt₃₄ are smaller than the time delays between other image combinations, tending towards 0 when β approaches 0. The same holds true for the SPT-modified time delays, while we note that $\text{Δ}{\tilde{t}}_{12}\approx \text{Δ}{t}_{12}$ and $\text{Δ}{\tilde{t}}_{34}\approx \text{Δ}{t}_{34}$. For sources β ∈ B₁, the mean impact of the SPT, denoted as $\langle \text{Δ}\tilde{t}/\text{Δ}t\rangle$, is around of a few percent for both the NISg and the NIE. In contrast to the NISg, the impact of the SPT for the case of the NIE is larger in the quadruple image regime than in the double image regime. This only reflects that $1+f\left(\left|{\beta }_{\text{max}}^{\text{NIE}}\right|\right)\lesssim {\langle \text{Δ}\tilde{t}/\text{Δ}t\rangle }_{}^{\text{NIE}}$ while $1+f\left(\left|{\beta }_{\text{max}}^{\text{NISg}}\right|\right)\gg {\langle \text{Δ}\tilde{t}/\text{Δ}t\rangle }_{}^{\text{NISg}}$.

Fig. 7 Map of |Δθ(β)|max in units of θE for f2 = 0.4 and γp = 0.1 (see Eq. (42)). The inner (resp. outer) green line represents the tangential (resp. radial) caustic curve. The two hatched regions (B1 and B2) delimited by blue curves demarcate parts of the source plane where |Δθ(β)|max ≤ 5 × 10−3 θE. The region B1 lies inside the tangential caustic curve while B2 lies outside. The inset highlight the region around the position β = 0.

Fig. 8 Grid of sources β located inside the region B1 ∪B2 for the NIE with (θc, ϵ) = (0.1 θE, 0.15). The inner (resp. outer) green line represents the tangential (resp. radial) caustic curve. The color-coding refers to the offsets |Δ θ |∕θE between the lensed images θ of the sources β and the lensed images $\tilde{\theta }$ of the SPT-modified sources $\widehat{\beta }$.

Fig. 9 Time delay ratios of image pairs between the NISg and the corresponding SPT-modified model. The model parameters are (θc, γp) = (0.1 θE, 0.1) and the radial stretching is characterized by f2 = 0.4. Top panel: $\text{Δ}\tilde{t}/\text{Δ}t$ normalized by the ratio $\left|\widehat{\beta }\right|/\left|\beta \right|$ is close to 1 in the double image regime and for opposite images in the quadruple image regime. Bottom panel: the impact of the SPT on the time delays is around a few percent, reaching a maximum of 12% for the particular case of a source located almost on the radial caustic curve but still inside B2.

Fig. 10 Time delay ratios of image pairs between the NIE and the corresponding SPT-modified model. The model parameters are (θc, ϵ) = (0.1 θE, 0.15) and the radial stretching is characterized by f2 = 0.4. The time delay ratios in the quadruple (resp. double) image regime are depicted with triangles (resp. inverted triangles). The red line shows the source position ratios $\left|\widehat{\beta }\right|/\left|\beta \right|$ and the color-coding refers to the azimuth angle ∠β. Top panel: even though the dispersion is slightly larger compared to the quadrupole model, $\text{Δ}\tilde{t}/\text{Δ}t$ normalized by the ratio $\left|\widehat{\beta }\right|/\left|\beta \right|$ is still close to 1 in the double image regime and for opposite images in the quadruple image regime. Bottom panel: the time delay ratio dispersion is clearly related to ∠β (see the text for more details). The impact of the SPT on the time delays is also around a few percent, reaching a maximum of around 5% for aparticular combination of non-opposite images of a source located inside B1.

4.3 The alternative SPT-modified time delays

In Sect. 4.1, we have shown that the capability of the deflection law $\tilde{\alpha }$ to predict the same multiple images as predicted by α (with an accuracy of ε_acc) is very limited (see Eq. (41)). Based upon a representative example, Fig. 7 shows that only a very small region (B₁) in the source plane leads to indistinguishable quadruple image configurations. In this section, we investigate a method to extend the region in the lens plane where the offsets |Δθ| are smaller than ε_acc. The idea consists in finding source positions β^⋆ in the vicinity of $\widehat{\beta }(\beta )$ that lead to new image positions ${\theta }^{\star }={\beta }^{\star }+\tilde{\alpha }\left({\theta }^{\star }\right)$ in such a way that the offsets |Δ θ ^⋆||θ^⋆ −θ| are as small as possible. The search for each β^⋆ is based on the numerical minimization of the cost function h(β^⋆) defined by $$h({\beta }^{\star })=\text{max}\left(\left|{\theta }^{\star }{}_i\left({\beta }^{\star }\right)-{\theta }_i\left(\widehat{\beta }\right)\right|\right),$$(48)

using the Levenberg–Marquardt algorithm (Levenberg 1944). Because we expect β ^⋆ to be close to the corresponding $\widehat{\beta }$, we always choose the latter as first guesses while we do not restrict β^⋆ to share the same direction as $\widehat{\beta }$. Thus, the resulting source mapping ${\beta }^{\star }(\widehat{\beta })$, hence β^⋆(β), may not be isotropic as for the radial stretching.

As we shall see, this approach may drastically increase the region B₁ ∪B₂ (in particular B₁) while benefiting from a simple implementation. However, we must point out that this method suffers several flaws. First, the source mapping β ^⋆ (β) lacks a solid analytical basis. While $\widehat{\beta }(\beta )$ gives rise to $\widehat{\alpha }$ which is analytically connected to $\tilde{\alpha }$ by definition⁴, there is no apparent link between $\tilde{\alpha }$ and β^⋆. Moreover, the way it is obtained precludes any further analytical investigation. Second, there is no definitive guarantee for the source mapping β ^⋆ (β) to be one-to-one over $\mathcal{U}$. Finally, a successful minimization of the cost function h for a given β^⋆ does not guarantee |Δθ^⋆| to be smaller than ε_acc. Indeed, the solution β ^⋆ only corresponds to the one for which the cost function h(β^⋆) is the smallest in the vicinity of $\widehat{\beta }$, being potentially larger than ε_acc. This is particularly true for sources which are located very close to the caustic curves. At least, we have |Δ θ ^⋆|≤|Δθ| where the equality holds when ${\beta }^{\star }=\widehat{\beta }$. For these reasons, we point out that this numerical approach cannot supplant the analytical reconsideration of how the curl-free deflection law $\tilde{\alpha }$ is defined. However, the combination of β^⋆(β) and $\tilde{\alpha }$ constitutes a physically meaningful alternative to β and α, and deserves to be considered.

To illustrate the method, we adopt the same lens model and SPT as in Sect. 4.1. We also consider the same grid of sources β covering the first quadrant of the source planeand restricted to 0 ≤|β|≤ 0.66 θ_E. We illustrate the results of the numerical optimization in Figs. 11 and 12. The ratio |β ^⋆ |∕|β| plotted against |β|∕θ_E in Fig. 11 clearly shows the slight anisotropy of the source mapping β ^⋆ (β) resulting from the numerical optimization. Indeed, for original sources located on a quarter circle with a radius |β |, the corresponding |β ^⋆ | depend on the azimuth angle ∠β; |β ^⋆ | is larger than $\left|\widehat{\beta }\right|$ for ∠β > 45° and smaller for ∠β < 45°. For sources leading to double image configurations, the |β^⋆| scatter is around 1%. Furthermore, discontinuities in the mapping β^⋆(β) appear when a source crosses the tangential caustic curve. Two particular jumps are highlighted for sources passing through the two cusps. The ones located on the β_x-axis and β_y-axis are explicitely depicted in Fig. 11. Due to these discontinuities, an extended source which crosses the tangential caustic curve is not mapped smoothly onto an SPT-modified extended source. This effect propagates to the image plane, but the impact on the corresponding lensed image is not observable, as shown in Fig. 12. Similar to Fig. 4, the bottom panel in Fig. 12 represents the image positions θ and the color-coding shows the offsets |Δθ^⋆| in unit of θ_E. The same color-coding is also applied to the sources where only the largest offset between the corresponding pairs of lensed images are shown (see top panel in Fig. 12). Figure 12 shows that almost all the offsets |Δ θ ^⋆| are now smaller than ε_acc, even for sources located inside the tangential caustic curve. It is worth stating that a finer source grid would have led to a larger number of sources located very close to the caustic curves, for which the optimized cost function may be larger than ε_acc. Compared to Fig. 4, some image positions depicted in Fig. 12 show an offset |Δ θ ^⋆| (after the optimization process) larger than the offset |Δθ| (before the optimization process). For example, the image position θ_k = (1.242, 0) θ_E of the source β = (0.117, 0) θ_E is characterized by |Δθ_k| = 0.001 θ_E while $\left|\text{Δ}{\theta }^{\star }{}_k\right|=0.005{\theta }_{\text{E}}$. This behavior stems from the fact that, for a given n-image configuration, the optimization process minimizes only the largest offset but not all the n offsets simultaneously (because of the max(⋅) function in Eq. (48)). Thus, while the largest offset becomes smaller, the offset $\left|\text{Δ}{\theta }^{\star }{}_k\right|$ also varies during the optimization process, leading at the end to $\left|\text{Δ}{\theta }_k\right|<\left|\text{Δ}{\theta }^{\star }{}_k\right|\le h({\beta }^{\star })$.