© ESO 2019

1. Introduction

In strong gravitational lensing systems, in which a galaxy or a galaxy cluster causes multiple images or strong image distortions of background sources, one often neglects the inhomogeneities of the gravitational field between the observer and the lens, and between the lens and the source (see, e.g., Kochanek 2006; Treu 2010; Bartelmann 2010; Kneib & Natarajan 2011, for reviews on strong lensing systems). This usually provides a very good approximation, since the lensing strength of the main lens over the region where strong lensing effects occur is much larger than the typical distortion effects of matter along the line-of-sight. The latter is comparable to the typical strength of cosmic shear effects (see, e.g., Bartelmann & Schneider 2001; Schneider 2006; Munshi et al. 2008; Hoekstra 2013), and amounts to about 1 or 2% of the distortion in the strong-lens region.

Whereas the propagation effects are thus small, the interest of these weak distortions has been renewed, for at least two different reasons. The first is that strong lens systems may be biased toward showing relatively strong line-of-sight structures, thereby increasing their lensing probability. This effect is most likely affecting the lensing cross sections for the formation of giant arcs in clusters (see Bartelmann et al. 1998 for stating the “arc statistic problem”, and Meneghetti et al. 2013 for a recent update of the issue). Cosmological simulations (Puchwein & Hilbert 2009, and references therein) indicate that line-of-sight structure can indeed substantially modify the lensing efficiency of clusters. More recently, Bayliss et al. (2014) found a significant overdensity of galaxy groups along the line-of-sight toward strong-lensing clusters, observationally supporting the presence of this bias.

The second reason for the renewed interest in intervening distortions of strong-lens systems is their use for a precise determination of parameters, most noticible the Hubble constant from measured time delays (Refsdal 1964) in multiply imaged quasi-stellar objects (see, e.g., Kochanek 2003, 2006; Treu 2010; Suyu et al. 2010, 2013). The quality of modern imaging data and the accuracy of time delay estimates allows one to derive estimates of the Hubble constant with a formal error of ∼6% (Suyu et al. 2013). At this level of precision, line-of-sight effects may become highly relevant in these strong-lensing systems (e.g., Wong et al. 2011; Collett et al. 2013; Greene et al. 2013, and references therein).

With additional deflectors along the line-of-sight, the lens equation, which relates the true source position to the observed positions of images, needs to be generalized to include deflection at more than one distance from the observer. Blandford & Narayan (1986) established the theory of multi-plane gravitational lensing (see also Chap. 9 of Schneider et al. 1992), where the mapping between images and their source is affected by the action of several deflectors at different redshifts. This multi-plane lensing is employed in ray-tracing simulations – see Refsdal (1970) for the earliest ray-tracing simulations, Jain et al. (2000), Hilbert et al. (2009), and references therein – where the three-dimensional mass distribution between observer and source is partitioned into separate slices, and the mass distribution in each slice is treated as a gravitational lens plane.

The full theory of multiple deflection lensing is required if there are two or more strong deflectors along the same line-of-sight toward a background source. For galaxy-scale lensing, such systems are rare, since such an alignment is not very probable. However, in the sample of several hundred galaxy-scale strong lens systems currently known, there are examples of multiple lenses at two different redshifts (Chae et al. 2001; Gavazzi et al. 2008). These systems may be of particular interest, since they may be used in principle to determine cosmological distance ratios (Collett & Auger 2014), and therefore to constrain the cosmic expansion history, although in practice the mass-sheet transformation renders this a difficult task (Schneider 2014).

Far more common are situations where there is a single strong lens in the line-of-sight to a distant source, and several other deflectors at different redshifts situated sufficiently far away from the strong-lensing region of the main deflector such that their deflection angle can be linearized across this strong-lensing region. Kovner (1987) considered this case of one main lens, combined with linearized deflections at different redshifts. The effects of the linear deflectors can be summarized into a set of matrices which describe the mapping between angles at the vertices of light cones to the separation vectors along the light cone at the redshifts of lens and source. Schneider (1997) reconsidered this situation, using the general formalism of light propagation in an inhomogeneous universe (Seitz et al. 1994; hereafter SSE), and related this generalized gravitational lens equation to the one where the deflection occurs in only a single plane, with a tidal deflection matrix added (see also Keeton et al. 1997). McCully et al. (2014) further generalized this theory to consider several main lenses, together with linearized deflections between the main lens planes. In particular, McCully et al. (2014) emphasized the advantages of this hybrid framework for the modeling of strong lens systems with multiple main deflectors along the line-of-sight (see also McCully et al. 2017).

One of the prime motivations for the work of McCully et al. (2014) was the derivation of the time-delay function for such generalized lensing systems, as this is required for relating the time delay to the scale-length of the Universe, i.e., the Hubble radius. Using the Millennium Simulation, Jaroszynski & Kostrzewa-Rutkowska (2014) investigated the impact of the line-of-sight matter distribution on strong lensing properties, including the time delay; they concluded that the intermediate mass distribution leads to a spread of ∼6% in the product of Hubble constant and time delay, for strong lensing systems with source redshifts z_s ∼ 2. Arguably, the largest obstacle for the time-delay method to obtain accurate estimates for the Hubble constant is the degeneracy of the mass model due to the mass-sheet transformation (MST; Falco et al. 1985; see Schneider & Sluse 2013 for a recent discussion) or the more general source position transformation (Schneider & Sluse 2014). Recently, Schneider (2014) has shown that an MST also exists for the case of lenses at two different distances from us. In particular, this MST leads to a scaling of all time delays (from sources on both source planes) by the same factor, thus precluding that the degeneracy due to the MST can be broken by time-delay ratios.

In this paper, these results are further generalized to the case of arbitrarily many main lens planes, together with linear deflections between the main lens (and source) planes. After a short summary of the propagation equations in an inhomogeneous universe, as they apply to the case under consideration here, we derive the time-delay function for the case of several main lens planes. Whereas this has been done before by McCully et al. (2014), our expression for the light travel time function is expressed in a form which allows us to apply Fermat’s theorem in gravitational lensing. Thus, the lens equation can be obtained from requiring that the light travel time function is stationary with respect to all impact vectors in the main lens planes. With this procedure, we derive an alternative form of the lens equation, which involves the impact vectors of three consecutive lens planes and the deflection angle in the middle one of them. We show that this iterative form of the lens equation is equivalent to the more standard form; in order to do so, we first obtain a very general relation between distance matrices.

We then turn to the MST in this case, and find a curious property of its behavior: The uniform, isotropic scaling factor which characterizes the MST alternates between a free parameter λ and unity from one lens/source plane to the next. The corresponding modification of the deflection angle in the main lens planes corresponds to a scaling of the deflection, plus the addition of a tidal deflection matrix (which in the absence of linear deflections between lens planes reduces to the addition of a uniform mass sheet). Finally we show that all time delays, for sources located on any plane, scale by the same factor under an MST, precluding the possibility of breaking the degeneracy from the MST by measured time-delay ratios.

2. Generalized multi-plane lens equation

2.1. Optical tidal equation

The propagation of light rays follows from the geodesic equation, specialized to a perturbed Robertson–Walker metric (see Schneider et al. 1992; Seitz et al. 1994; Bartelmann 2010, and references therein). In particular, for infinitesimally small light bundles, the separation vector ξ of a ray from the reference (or “central”) ray of the bundle is given by the optical tidal equation

$$\begin{array}{c}\hfill {\mathit{\xi }}^{″}(\lambda )=\mathcal{T}(\lambda )\mathit{\xi }(\lambda ),\end{array}$$(1)

where 𝒯 is the optical tidal matrix, evaluated at the affine parameter λ along the central ray. Here and throughout this paper, a prime denotes differentiation with respect to the affine parameter λ, which is related to the redshift z, the cosmic scale factor a = 1/(1 + z), or cosmic time t through

$$\begin{array}{c}\hfill \mathrm{d}\lambda =-\frac{c\mathrm{d}t}{(1+z)}=-\frac{c\mathrm{d}a}{H}=\frac{c\mathrm{d}z}{H{(1+z)}^2},\end{array}$$(2)

where we have assumed that λ increases with redshift, i.e., decreases with cosmic time,

$$\begin{array}{c}\hfill H=\dot{a}/a=H_0\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}{(1+z)}^3+(1-{\mathrm{\Omega }}_{\mathrm{m}}-{\mathrm{\Omega }}_{\mathrm{\Lambda }}){(1+z)}^2+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}\end{array}$$(3)

is the Hubble function, H₀ is the Hubble constant, and Ω_m and Ω_Λ are the cosmic density parameters in matter and vacuum energy.

We consider 𝒯 to consist of three separate components, 𝒯 = 𝒯_bg + 𝒯_sm + 𝒯_cl. The first is the optical tidal matrix of the homogeneous background universe, and is given by (see SSE)

$$\begin{array}{c}\hfill {\mathcal{T}}_{\mathrm{bg}}=-\frac{3}{2}{\left(\frac{H_0}{c}\right)}^2{\mathrm{\Omega }}_{\mathrm{m}}{(1+z)}^5\mathcal{I},\end{array}$$(4)

where ℐ is the two-dimensional unit matrix. If we consider a light bundle with vertex at redshift z = 0 and affine parameter λ = 0, which is only subject to 𝒯_bg, then the solution of Eq. (1) is ξ(λ)=D(λ)θ, where θ is the angle that the light ray under consideration encloses with the fiducial ray at the vertex, and D(λ) is the solution of the differential equation

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^{2}D(\lambda )}{\mathrm{d}{\lambda }^2}=-\frac{3}{2}{\left(\frac{H_0}{c}\right)}^2{\mathrm{\Omega }}_{\mathrm{m}}{(1+z)}^{5}D(\lambda ),\end{array}$$(5)

with initial condition D(0)=0 and D′=1. Here, D is the angular-diameter distance of the homogeneous universe, as a function of affine parameter (or redshift).

We consider two different kinds of inhomogeneities here. The first of them is related to small-scale density inhomogeneities, such as galaxies and their dark matter halos. For those, the optical tidal matrix 𝒯_cl is a strong function of position, and there is no longer a linear relation between enclosed angle θ and separation ξ of a light ray, for finite θ. In strong-lensing applications, we are typically interested only in that region of these small-scale inhomogeneities where multiple images can be formed, which corresponds to a few tens of kiloparsecs for galaxy lenses. These small-scale inhomogeneities will be considered explicitly as “main lenses” in the following. The second kind of inhomogeneities is due to the large-scale mass distribution of the universe. We assume that the corresponding gravitational field is sufficiently smooth, so that the tidal effects can be considered approximately constant across the region where strong-lensing effects occur. Hence we assume that over such a region, 𝒯_sm can be considered to depend only on λ, not on the actual position of the light ray. According to SSE, this contribution of the optical tidal matrix is

$$\begin{array}{c}\hfill {\left({\mathcal{T}}_{\mathrm{sm}}\right)}_{\mathit{ij}}=-\frac{{(1+z)}^2}{c^2}(2\frac{{\partial }^2\varphi }{\partial {\xi }_i\partial {\xi }_j}+{\delta }_{\mathit{ij}}\frac{{\partial }^2\varphi }{\partial {\xi }_3^2}),\end{array}$$(6)

where ϕ is the Newtonian potential sourced by the density inhomogeneity, i.e., satisfying the (three-dimensional) Poisson equation ${\mathrm{\nabla }}_{\xi }^2\varphi =4\pi G(\rho -\overline{\rho })$, where $\overline{\rho }(z)$ is the mean matter density in the Universe, and we assumed that the light ray propagates in the ξ₃-direction.

2.2. Generalized multi-plane lens equation

We consider a direction in the sky with has a set of main lenses along the line-of-sight to a distant source, located at redshifts z_i or affine parameters λ_i. As shown in SSE (see also Schneider 1997), the separation vector ξ(λ) then becomes

$$\begin{array}{c}\hfill \mathit{\xi }(\lambda )=\mathcal{D}(\lambda )\mathit{\theta }-{\displaystyle {\sum }_i}{\mathcal{D}}_i(\lambda )[\widehat{{\mathit{\alpha }}_i}({\mathit{\xi }}_i)-{\widehat{\mathit{\alpha }}}_i^0]\mathrm{H}(\lambda -{\lambda }_i),\end{array}$$(7)

where θ is the angle the light encloses with the fiducial ray at the observer. The distance matrix 𝒟(λ) solves the optical tidal equation

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\mathcal{D}(\lambda )}{\mathrm{d}{\lambda }^2}=[{\mathcal{T}}_{\mathrm{bg}}(\lambda )+{\mathcal{T}}_{\mathrm{sm}}(\lambda )]\mathcal{D}(\lambda ),\end{array}$$(8)

with initial conditions 𝒟(0) = 0 and 𝒟′(0) = ℐ, and the distance matrices 𝒟_i(λ) solve the same differential equation, but with initial conditions 𝒟_i(λ_i) = 0 and (d𝒟_i/dλ)(λ_i) = (1+z_i)ℐ. They are the distance matrices which apply for light rays having their vertex at λ_i. In Eq. (7), H(λ − λ_i) is the Heaviside step function. The deflection angle ${\widehat{\mathit{\alpha }}}_i({\mathit{\xi }}_i)$ is given in terms of the surface mass density Σ_i(ξ_i) as

$$\begin{array}{c}\hfill {\widehat{\mathit{\alpha }}}_i({\mathit{\xi }}_i)=\frac{4G}{c^2}{\displaystyle \int {\mathrm{d}}^2{\xi }^{\prime }{\mathrm{\Sigma }}_i({\mathit{\xi }}^{\prime })\frac{{\mathit{\xi }}_i-{\mathit{\xi }}^{\prime }}{|{\mathit{\xi }}_i-{\mathit{\xi }}^{\prime }{|}^2},}\end{array}$$(9)

and ${\widehat{\mathit{\alpha }}}_i^0$ denotes the deflection angle of the fiducial ray in the ith lens plane. If we now perform the translation $\mathit{\xi }(\lambda )=\stackrel{\sim }{\mathit{\xi }}(\lambda )+\mathit{\eta }(\lambda )$, with

$$\mathit{\eta }(\lambda )={\displaystyle {\sum }_i}{\mathcal{D}}_i(\lambda ){\widehat{\mathit{\alpha }}}_i^0\mathrm{H}(\lambda -{\lambda }_i),$$

then $\stackrel{\sim }{\mathit{\xi }}$ satisfies Eq. (7) with the ${\widehat{\mathit{\alpha }}}_i^0$ set to zero. In the following, we always assume this (unobservable) translation, and drop the tilde on $\stackrel{\sim }{\mathit{\xi }}$ henceforth. For the impact vectors ξ_j in the jth plane, we then obtain

$$\begin{array}{c}\hfill {\mathit{\xi }}_j=\mathcal{D}({\lambda }_j)\mathit{\theta }-{\displaystyle {\sum }_{i=1}^{j-1}{\mathcal{D}}_i({\lambda }_j){\widehat{\mathit{\alpha }}}_i({\mathit{\xi }}_i)\equiv {\mathcal{D}}_j\mathit{\theta }-{\sum }_{i=1}^{j-1}{\mathcal{D}}_{\mathit{ij}}{\widehat{\mathit{\alpha }}}_i({\mathit{\xi }}_i).}\end{array}$$(10)

2.3. Calculation of the distance matrices

The distance matrices 𝒟i(λ) depend on the large-scale matter distribution around the line-of-sight, which is dominated by dark matter and thus difficult to determine observationally. One possibility to estimate 𝒯_sm from observations is to assume that galaxies provide a good tracer of the total matter distribution on large scales. From an observed distribution of galaxies (or a particular kind of galaxies, like luminous red galaxies) around the line-of-sight, an estimate of the tidal field can be obtained; this is the strategy proposed in Collett et al. (2013) and Smith et al. (2014).

The propagation of light through the large-scale structure is the subject of cosmological weak lensing, or cosmic shear (see, e.g., Schneider 2006). In cosmic shear, one usually describes the propagation matrices in comoving coordinates. In order to connect the formulation given here, where the 𝒟_i relate angles to proper transverse separation vectors (which is appropriate, as the matter distribution of the main lenses – galaxies or clusters – are most conveniently described in physical scales), to that used in weak lensing, we show in the appendix that

$$\begin{array}{c}\hfill {\mathcal{D}}_i(\chi )=D_i(\chi )\mathcal{I}-\frac{2}{c^2}{\displaystyle {\int }_{{\chi }_i}^{\chi }\mathrm{d}{\chi }^{\prime }\frac{a(\chi )}{a({\chi }^{\prime })}{f}_k(\chi -{\chi }^{\prime })\mathsf{H}(\varphi ({\chi }^{\prime })){\mathcal{D}}_i({\chi }^{\prime }),}\end{array}$$(11)

where χ is the comoving distance, f_k(χ) is the comoving angular diameter distance, D_i(χ)=a(χ)f_k(χ − χ_i), which satisfies the differential equation

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{f}_k(\chi )}{\mathrm{d}{\chi }^2}=-K{f}_k(\chi )\end{array}$$(12)

with initial conditions f_k(0)  =  0 and df_k(0)/dχ  =  1. Furthermore, $K=({\mathrm{\Omega }}_{\text{m}}+{\mathrm{\Omega }}_{\mathrm{\Lambda }}-1)H_0^2/c^2$ is the spatial curvature of the universe, and H(ϕ) is the two-dimensional Hessian of the gravitational potential ϕ, evaluated in comoving transverse coordinates. Provided the perturbations are small, so that 𝒟_i deviates only slightly from D_iℐ, we can replace 𝒟_i(χ′) by D_i(χ′) in the integrand,

$$\begin{array}{c}\hfill {\mathcal{D}}_i(\chi )=D_i(\chi )\mathcal{I}-\frac{2}{c^2}{\displaystyle {\int }_{{\chi }_i}^{\chi }\mathrm{d}{\chi }^{\prime }\frac{a(\chi )}{a({\chi }^{\prime })}{f}_k(\chi -{\chi }^{\prime })\mathsf{H}(\varphi ({\chi }^{\prime }))D_i({\chi }^{\prime }).}\end{array}$$(13)

In weak lensing, this approximation is called “the neglect of lens-lens coupling”, often also termed as “Born-approximation”. This approximation is very accurate and certainly sufficient for the purposes discussed in the current context¹. Thus, the deviation of 𝒟_i from D_iℐ is given as a line-of-sight integral over the tidal force field, with a distance-dependent weighting. We also note that in the approximation (13), the distance matrices 𝒟_i are symmetric, which is generally not the case if the exact expression (11) is used.

For statistical studies, instead of trying to obtain the tidal field along the line-of-sight toward the sources from observations of the galaxy distribution in those directions, one can also derive the probability distribution for the distance matrices from cosmological simulations, as has been done in Suyu et al. (2013).

3. Time-delay function and the iterative lens equation

In this section, we first derive the light travel time function (LTTF) corresponding to the generalized multi-plane lens equation considered in the previous section. An alternative form of the lens equation is then derived from Fermat’s principle, which in the current context states that the lens equation is equivalent to setting the gradient of the LTTF with respect to all impact vectors equal to zero (i.e., real light rays correspond to stationary points of the LTTF). Then we show that this iterative form of the lens equation is equivalent to Eq. (10).

3.1. Time-delay function

We first consider a single lens plane at z₁ with a source at z₂, in which case the lens equation reads

$$\begin{array}{c}\hfill {\mathit{\xi }}_2={\mathcal{D}}_2\mathit{\theta }-{\mathcal{D}}_{12}\widehat{\mathit{\alpha }}({\mathit{\xi }}_1)={\mathcal{D}}_2{\mathcal{D}}_1^{-1}{\mathit{\xi }}_1-{\mathcal{D}}_{12}\widehat{\mathit{\alpha }}({\mathit{\xi }}_1),\end{array}$$(14)

where we set $\widehat{\mathit{\alpha }}({\mathit{\xi }}_1)\equiv {\widehat{\mathit{\alpha }}}_1({\mathit{\xi }}_1)$. We define the LTTF τ(ξ₁, ξ₂) to be the excess light travel time from a source at ξ₂ to the observer caused by the deflection in the lens plane at ξ₁. This excess travel time has two components, a geometrical one (a bent ray is longer than an unbent one), and a potential one caused by the retardation of photons in the gravitational potential of the deflector.

In standard lens theory, with unperturbed angular-diameter distances D_ij, the potential part of the time delay function takes the form cτ_pot = −(1 + z₁)(D₁D₂/D₁₂)ψ(θ), where ψ(θ) is the deflection potential, which satisfies ${\mathrm{\nabla }}_{\mathit{\theta }}\psi =\mathit{\alpha }(\mathit{\theta })\equiv (D_{12}/D_2)\widehat{\mathit{\alpha }}(D_1\mathit{\theta })$. Now we define the potential $\widehat{\psi }({\mathit{\xi }}_1)$ such that it satisfies ${\mathrm{\nabla }}_{{\mathit{\xi }}_1}\widehat{\psi }=\widehat{\mathit{\alpha }}$. This potential is a multiple of ψ, i.e., $\widehat{\psi }({\mathit{\xi }}_1)=k\psi ({\mathit{\xi }}_1/D_1)$, up to an irrelevant additive constant. To find k, we consider

$$\widehat{\mathit{\alpha }}={\mathrm{\nabla }}_{{\mathit{\xi }}_1}\widehat{\psi }=\frac{k}{D_1}{\mathrm{\nabla }}_{\mathit{\theta }}\psi =\frac{k}{D_1}\mathit{\alpha }=\frac{k{D}_{12}}{D_1{D}_2}\widehat{\mathit{\alpha }},$$

yielding

$$\begin{array}{c}\hfill \widehat{\psi }({\mathit{\xi }}_1)=\frac{D_1{D}_2}{D_{12}}\psi ({\mathit{\xi }}_1/D_1).\end{array}$$(15)

Hence, the potential part of the LTTF takes the form $c{\tau }_{\mathrm{pot}}=-(1+z_1)\widehat{\psi }({\mathit{\xi }}_1)$. As expected, it does not depend on the cosmological distances, since it is caused solely by the local effect of propagating through a gravitational field, and the corresponding time interval is then redshifted by a factor (1 + z₁) to the observer. Accordingly, also in the case of perturbed light propagation between the lens planes, the potential part of the time delay must have the same form, since it is unaffected by propagation effects.

We find an explicit expression for τ(ξ₁, ξ₂) by requiring that ∇_ξ1τ(ξ₁, ξ₂)=0 is equivalent to the lens Eq. (14). This fixes τ up to a multiplicative constant and terms which depend solely on ξ₂. The multiplicative constant is fixed by the explicit expression given above for the potential part of τ, and the additive constant (which is irrelevant for time delay measurements) is fixed by requiring that the geometrical part of τ should vanish if the light ray is undeflected by the main lens. This then yields

$$\begin{array}{cc}\hfill c\tau ({\mathit{\xi }}_1,{\mathit{\xi }}_2)=& \frac{1+z_1}{2}{({\mathsf{C}}_{12}{\mathit{\xi }}_1-{\mathcal{D}}_{12}^{-1}{\mathit{\xi }}_2)}^{\mathrm{t}}{\mathsf{C}}_{12}^{-1}({\mathsf{C}}_{12}{\mathit{\xi }}_1-{\mathcal{D}}_{12}^{-1}{\mathit{\xi }}_2)\hfill \\ \hfill & -(1+z_1)\widehat{\psi }({\mathit{\xi }}_1),\hfill \end{array}$$(16)

where we defined

$$\begin{array}{c}\hfill {\mathsf{C}}_{\mathit{ij}}={\mathcal{D}}_{\mathit{ij}}^{-1}{\mathcal{D}}_j{\mathcal{D}}_i^{-1}.\end{array}$$(17)

As shown in Schneider (1997) – see also Kovner (1987) – the matrix C_ij is symmetric. Indeed, we see that

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{{\mathit{\xi }}_1}\tau ({\mathit{\xi }}_1,{\mathit{\xi }}_2)=\frac{1+z_1}{c}({\mathsf{C}}_{12}{\mathit{\xi }}_1-{\mathcal{D}}_{12}^{-1}{\mathit{\xi }}_2-\widehat{\mathit{\alpha }})=0\end{array}$$(18)

is equivalent to the lens Eq. (14), as is easily verified by multiplying the foregoing expression by 𝒟₁₂ from the left. We also note that τ = 0 if $\widehat{\psi }=0$ and if the ray is unbent, i.e., if ${\mathit{\xi }}_2={\mathcal{D}}_2{\mathcal{D}}_1^{-1}{\mathit{\xi }}_1$, as we required for τ. McCully et al. (2014) obtained a somewhat different form for τ, which they showed to be equivalent to the expression given here. However, it must be pointed out that this equivalence applies only to physical light rays, i.e., those which satisfy the gravitational lensing Eq. (14). For those rays, we could write the light travel time as

$$c\tau =\frac{1+z_1}{2}{\widehat{\mathit{\alpha }}}_1^{\mathrm{t}}{\mathsf{C}}_{12}^{-1}{\widehat{\mathit{\alpha }}}_1-(1+z_1)\widehat{\psi }(\mathit{\xi }),$$

where the lens Eq. (14) was used to eliminate ξ₂. However, Eq. (16) is more general, as it yields the light travel time for all kinematically possible rays, not only for those for which the bend by the main lenses equals the actual deflection angle as calculated as the gradient of $\widehat{\psi }$. This more general form of the τ is needed if the lens equation is to be derived from Fermat’s principle.

In case of several main lens planes, the LTTF is obtained by considering the replacement of the actual light ray by successively straighter rays, i.e., by removing the bends of the ray and the gravitational potentials they traverse (see Sect. 9.2 of Schneider et al. 1992). Removing the bend and deflection potential in the first plane leads to the contribution (16) of the LTTF. Subsequent removal of the bend and potential in the second lens plane yields a similar contribution, with the indices (1,2) replaced by (2,3). Iterating this consideration, we obtain for the general case

$$\begin{array}{cc}\hfill c\tau ({\mathit{\xi }}_1,{\mathit{\xi }}_2,\cdots ,{\mathit{\xi }}_{N+1})=& {\displaystyle {\sum }_{i=1}^N(1+z_i)[\frac{1}{2}{({\mathsf{C}}_{i,i+1}{\mathit{\xi }}_i-{\mathcal{D}}_{i,i+1}^{-1}{\mathit{\xi }}_{i+1})}^{\mathrm{t}}}\hfill \\ \hfill & \times {\mathsf{C}}_{i,i+1}^{-1}({\mathsf{C}}_{i,i+1}{\mathit{\xi }}_i-{\mathcal{D}}_{i,i+1}^{-1}{\mathit{\xi }}_{i+1})-{\widehat{\psi }}_i({\mathit{\xi }}_i)]\hfill \end{array}$$(19)

as the sum over terms of the form (16) for the individual planes.

Any ray connecting the source at ξ_N + 1 and the observer is fully characterized by the impact vectors ξ_i, 1 ≤ i ≤ N in the lens planes, since between the planes, it follows the propagation Eq. (1) whose solution is uniquely determined by the two impact vectors at consecutive planes. Therefore, the actual light rays are singled out as those for which the LTTF is stationary, with respect to variations of the impact vectors in the N lens planes. Thus, the lens equation is obtained by setting the derivative of τ with respect to the ξ_j equal to zero. For each j ≥ 2, two terms of the above sum contribute, namely the terms i = j and i = j − 1. We obtain

$$\begin{array}{cc}\hfill {\mathrm{\nabla }}_{{\mathit{\xi }}_j}(c\tau )=& (1+z_j)[{\mathsf{C}}_{j,j+1}{\mathit{\xi }}_j-{\mathcal{D}}_{j,j+1}^{-1}{\mathit{\xi }}_{j+1}-{\widehat{\mathit{\alpha }}}_j({\mathit{\xi }}_j)]\hfill \\ \hfill & +(1+z_{j-1}){\left({\mathcal{D}}_{j-1,j}^{-1}\right)}^{\mathrm{t}}({\mathsf{C}}_{j-1,j}^{-1}{\mathcal{D}}_{j-1,j}^{-1}{\mathit{\xi }}_j-{\mathit{\xi }}_{j-1})\hfill \\ \hfill =& 0,\hfill \end{array}$$(20)

or

$$\begin{array}{cc}\hfill {\mathit{\xi }}_{j+1}=& {\mathcal{D}}_{j,j+1}[{\mathsf{C}}_{j,j+1}+\frac{1+z_{j-1}}{1+z_j}{\left({\mathcal{D}}_{j-1,j}^{-1}\right)}^{\mathrm{t}}{\mathsf{C}}_{j-1,j}^{-1}{\mathcal{D}}_{j-1,j}^{-1}]{\mathit{\xi }}_j\hfill \\ \hfill & -\frac{1+z_{j-1}}{1+z_j}{\mathcal{D}}_{j,j+1}{\left({\mathcal{D}}_{j-1,j}^{-1}\right)}^{\mathrm{t}}{\mathit{\xi }}_{j-1}-{\mathcal{D}}_{j,j+1}{\widehat{\mathit{\alpha }}}_j({\mathit{\xi }}_j).\hfill \end{array}$$(21)

This equation relates the position vectors in three consecutive planes to the deflection angle in the middle plane, quite in contrast to the lens Eq. (10) which contains all impact vectors ξ_i for a given ξ_j, 1 ≤ i ≤ j − 1. Hence, this new lens equation is more “local” than the original one.

In the following we explicitly show that these two forms of the lens equation are equivalent. For this, we first need to derive a general relation between distance matrices.

3.2. A relation between distance matrices

Consider the pairs of light rays sketched in Fig. 1, where the first has a vertex at λ_q and encloses an angle θ with the fiducial ray. At the affine parameter λ_s, its separation vector from the fiducial ray is ξ_s. The second light ray has its vertex at λ_r and intersects the first ray at λ_s; this then specifies its direction ϑ relative to the fiducial ray. At the intersection point, the two rays enclose an angle φ. From the geometry of the figure, we find

$$\begin{array}{cc}\hfill \mathit{\theta }& ={\mathcal{D}}_{\mathit{qr}}^{-1}{\mathit{\xi }}_r={\mathcal{D}}_{\mathit{qs}}^{-1}{\mathit{\xi }}_s={\mathcal{D}}_{\mathit{qt}}^{-1}{\mathit{\xi }}_t;\hfill \\ \hfill \mathit{\vartheta }& ={\mathcal{D}}_{\mathit{rt}}^{-1}({\mathit{\xi }}_t+\mathrm{\Delta }{\mathit{\xi }}_t)={\mathcal{D}}_{\mathit{rs}}^{-1}{\mathit{\xi }}_s.\hfill \end{array}$$(22)

Fig. 1. Sketch of two light rays through four consecutive planes, with 0 ≤ λq <  λr <  λs <  λt. The rays are not deflected in the lens planes. The first ray has its vertex at λq and encloses an angle θ with the fiducial ray; the second ray with vertex at λr encloses an angle ϑ with the fiducial ray. Both rays intersect at λs. The geometry of this figure yields the relation (25) between distance matrices.

We use the latter equation to derive a relation between the 𝒟’s, by expressing all vectors in terms of ξ_s. Using the first of (22), we get ${\mathit{\xi }}_t={\mathcal{D}}_{\mathit{qt}}{\mathcal{D}}_{\mathit{qs}}^{-1}{\mathit{\xi }}_s$. Furthermore, the figure shows that Δξ_t = 𝒟_stφ. On the other hand, ξ_r = −𝒟_s(λ_r)φ.

We now have to relate 𝒟_s(λ_r), the backward extension of the solution 𝒟_s(λ) of Eq. (8), to 𝒟_rs. For that, we consider two solutions of Eq. (8), 𝒟_r(λ) and 𝒟_s(λ), with their appropriate initial conditions at λ_r and λ_s, respectively, and define the matrix

$$\begin{array}{c}\hfill \mathsf{W}(\lambda )=\frac{\mathrm{d}{\mathcal{D}}_r^{\mathrm{t}}}{\mathrm{d}\lambda }(\lambda ){\mathcal{D}}_s(\lambda )-{\mathcal{D}}_r^{\mathrm{t}}(\lambda )\frac{\mathrm{d}{\mathcal{D}}_s}{\mathrm{d}\lambda }(\lambda ),\end{array}$$(23)

which is the Wronskian of the differential Eq. (8); here, the superscript “t” denotes the transpose of a matrix. The derivative of W vanishes, due to Eq. (8); hence, W is a constant. Evaluating Eq. (23) at λ = λ_r and making use of the initial conditions of 𝒟_r(λ) yields 𝖶 = (1 + z_r) 𝒟_s(λr). Similarly, at λ = λ_s we find $\mathsf{W}=-(1+z_s)\left({\mathcal{D}}_r^{\mathrm{t}}\right)({\lambda }_s)=-(1+z_s){\mathcal{D}}_{\mathit{rs}}^{\mathrm{t}}$. Thus, we obtain

$$\begin{array}{c}\hfill {\mathcal{D}}_s({\lambda }_r)=-\frac{1+z_s}{1+z_r}{\mathcal{D}}_{\mathit{rs}}^{\mathrm{t}},\end{array}$$(24)

which is Etherington’s theorem in matrix form (Etherington 1933). With this relation, we then find that

$$\mathrm{\Delta }{\mathit{\xi }}_t={\mathcal{D}}_{\mathit{st}}\mathit{\phi }=-{\mathcal{D}}_{\mathit{st}}{\mathcal{D}}_s^{-1}({\lambda }_r){\mathit{\xi }}_r=\frac{1+z_r}{1+z_s}{\mathcal{D}}_{\mathit{st}}{\left({\mathcal{D}}_{\mathit{rs}}^{\mathrm{t}}\right)}^{-1}{\mathit{\xi }}_r.$$

Using Eq. (22) and collecting terms,

$${\mathcal{D}}_{\mathit{rt}}{\mathcal{D}}_{\mathit{rs}}^{-1}{\mathit{\xi }}_s={\mathit{\xi }}_t+\mathrm{\Delta }{\mathit{\xi }}_t={\mathcal{D}}_{\mathit{qt}}{\mathcal{D}}_{\mathit{qs}}^{-1}{\mathit{\xi }}_s+\frac{1+z_r}{1+z_s}{\mathcal{D}}_{\mathit{st}}{\left({\mathcal{D}}_{\mathit{rs}}^{\mathrm{t}}\right)}^{-1}{\mathcal{D}}_{\mathit{qr}}{\mathcal{D}}_{\mathit{qs}}^{-1}{\mathit{\xi }}_s$$

follows. Since this relation is valid for all ξ_s, a general relation between distance matrices is obtained:

$$\begin{array}{c}\hfill \frac{1+z_r}{1+z_s}{\mathcal{D}}_{\mathit{st}}{\left({\mathcal{D}}_{\mathit{rs}}^{\mathrm{t}}\right)}^{-1}{\mathcal{D}}_{\mathit{qr}}={\mathcal{D}}_{\mathit{rt}}{\mathcal{D}}_{\mathit{rs}}^{-1}{\mathcal{D}}_{\mathit{qs}}-{\mathcal{D}}_{\mathit{qt}},\end{array}$$(25)

where we multiplied the resulting equation by 𝒟_qs from the right². This result generalizes the corresponding relation in Schneider (2016) derived for the case that the distance matrices 𝒟 reduce to scalars.

Indeed, a relation of this kind is expected to hold: Consider λ_t ≡ λ as a variable. The two matrix-valued functions 𝒟_q(λ) and 𝒟_r(λ) are linearly independent solutions of the transport Eq. (8), provided λ_r ≠ λ_q. Therefore, the solution 𝒟_s(λ) can be written as a linear combination of the other two. This combination should be of the form

$$\begin{array}{c}\hfill {\mathcal{D}}_s(\lambda )=[{\mathcal{D}}_r(\lambda ){\mathcal{D}}_r^{-1}({\lambda }_s)-{\mathcal{D}}_q(\lambda ){\mathcal{D}}_q^{-1}({\lambda }_s)]\mathsf{X},\end{array}$$(26)

which satisfies one of the initial conditions, 𝒟_s(λ_s) = 0. The matrix X is determined from the second initial condition; our result (25) shows that

$$\begin{array}{c}\hfill {\mathcal{D}}_s(\lambda )=[{\mathcal{D}}_r(\lambda ){\mathcal{D}}_r^{-1}({\lambda }_s)-{\mathcal{D}}_q(\lambda ){\mathcal{D}}_q^{-1}({\lambda }_s)]\frac{1+z_s}{1+z_r}{\mathcal{D}}_{\mathit{qs}}{\mathcal{D}}_{\mathit{qr}}^{-1}{\mathcal{D}}_{\mathit{rs}}^{\mathrm{t}}.\end{array}$$(27)

3.3. Equivalence of Eqs. (10) and (21)

We shall now show that the two forms (10) and (21) of the lens equation are equivalent. As a first step, we rewrite Eq. (21) in a form that admits a simple geometrical interpretation. Specializing Eq. (25) to q = 0, r = j − 1, s = j, t = j + 1 yields

$$\begin{array}{c}\hfill \frac{1+z_{j-1}}{1+z_j}{\mathcal{D}}_{j,j+1}{\left({\mathcal{D}}_{j-1,j}^{\mathrm{t}}\right)}^{-1}{\mathcal{D}}_{j-1}={\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}{\mathcal{D}}_j-{\mathcal{D}}_{j+1}.\end{array}$$(28)

We next consider the prefactor of ξ_j in Eq. (21). Using ${\mathsf{C}}_{j-1,j}^{-1}{\mathcal{D}}_{j-1,j}^{-1}={\mathcal{D}}_{j-1}{\mathcal{D}}_j^{-1}$, which is obtained from the definition (17) of C, we find that this prefactor becomes

$$\begin{array}{c}\hfill {\mathcal{D}}_{j+1}{\mathcal{D}}_j^{-1}+\frac{1+z_{j-1}}{1+z_j}{\mathcal{D}}_{j,j+1}{\left({\mathcal{D}}_{j-1,j}^{-1}\right)}^{\mathrm{t}}{\mathcal{D}}_{j-1}{\mathcal{D}}_j^{-1}={\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1},\end{array}$$(29)

where in the last step we made use of Eq. (28) and the fact that the inversion and transposition operations on a matrix commute. Thus, we can rewrite Eq. (21) in the form

$$\begin{array}{cc}\hfill {\mathit{\xi }}_{j+1}=& {\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}{\mathit{\xi }}_j\hfill \\ \hfill & -\frac{1+z_{j-1}}{1+z_j}{\mathcal{D}}_{j,j+1}{\left({\mathcal{D}}_{j-1,j}^{-1}\right)}^{\mathrm{t}}{\mathit{\xi }}_{j-1}-{\mathcal{D}}_{j,j+1}{\widehat{\mathit{\alpha }}}_j({\mathit{\xi }}_j).\hfill \end{array}$$(30)

Fig. 2. Propagation of a light ray (thick bent line) between three consecutive planes. The vertical line is the optical axis, with respect to which the separation vectors ξ are measured. The geometry of this figure yields the lens Eq. (30) – see text.

We note that this generalizes equation (4.47) of SSE to the case of general distance matrices between main lens planes. This form of the lens equation can be immediately interpreted geometrically. For this, we consider Fig. 2, from which we read off

$$\mathit{\eta }+{\mathit{\xi }}_{j+1}={\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}{\mathit{\xi }}_j-{\mathcal{D}}_{j,j+1}{\widehat{\mathit{\alpha }}}_j.$$

Furthermore, η = 𝒟_j,j+1φ; on the other hand, ξ_j−1 = −𝒟_j,j−1φ. Eliminating φ from these two relations and making use of Eq. (24), we find

$$\mathit{\eta }=\frac{1+z_{j-1}}{1+z_j}{\mathcal{D}}_{j,j+1}{\left({\mathcal{D}}_{j-1,j}^{-1}\right)}^{\mathrm{t}}{\mathit{\xi }}_{j-1}.$$

Together, these two equations reproduce Eq. (30). In this form, the equation not only is confined to three consecutive lens planes, but all distance matrices occurring here are those between these three planes.

In order to show the equivalence of Eqs. (10) and (21), it is useful to rewrite the prefactor of ξ_j − 1 in Eq. (30) in a different form. Making use again of Eq. (28), we obtain

$$\begin{array}{cc}\hfill {\mathit{\xi }}_{j+1}=& {\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}{\mathit{\xi }}_j\hfill \\ \hfill & +({\mathcal{D}}_{j+1}{D}_{j-1}^{-1}-{\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}{\mathcal{D}}_j{\mathcal{D}}_{j-1}^{-1}){\mathit{\xi }}_{j-1}\hfill \\ \hfill & -{\mathcal{D}}_{j,j+1}{\widehat{\mathit{\alpha }}}_j({\mathit{\xi }}_j).\hfill \end{array}$$(31)

We prove the equivalence by induction; for j = 1, this equivalence is seen from Eq. (18). Hence we assume that it is true for all planes up to j. Then, taking the difference between Eq. (10) for ξ_j + 1 and Eq. (31),

$$\begin{array}{cc}\hfill \mathrm{\Delta }=& {\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}{\mathit{\xi }}_j\hfill \\ \hfill & +({\mathcal{D}}_{j+1}{D}_{j-1}^{-1}-{\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}{\mathcal{D}}_j{\mathcal{D}}_{j-1}^{-1}){\mathit{\xi }}_{j-1}-{\mathcal{D}}_{j,j+1}{\widehat{\mathit{\alpha }}}_j({\mathit{\xi }}_j)\hfill \\ \hfill & -{\mathcal{D}}_{j+1}\mathit{\theta }+{\displaystyle {\sum }_{i=1}^j{\mathcal{D}}_{i,j+1}{\widehat{\mathit{\alpha }}}_i({\mathit{\xi }}_i),}\hfill \end{array}$$(32)

we need to show that Δ = 0. We first replace ξ_j − 1 and ξ_j by their expressions from Eq. (10), which holds because of the induction assumption,

$$\begin{array}{cc}\hfill \mathrm{\Delta }=& {\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}({\mathcal{D}}_j\mathit{\theta }-{\displaystyle {\sum }_{i=1}^{j-1}{\mathcal{D}}_{i,j}{\widehat{\mathit{\alpha }}}_i({\mathit{\xi }}_i)})\hfill \\ \hfill & +({\mathcal{D}}_{j+1}-{\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}{\mathcal{D}}_j)(\mathit{\theta }-{\displaystyle {\sum }_{i=1}^{j-2}{\mathcal{D}}_{j-1}^{-1}{\mathcal{D}}_{i,j-1}{\widehat{\mathit{\alpha }}}_i({\mathit{\xi }}_i)})\hfill \\ \hfill & -{\mathcal{D}}_{j,j+1}{\widehat{\mathit{\alpha }}}_j({\mathit{\xi }}_j)-{\mathcal{D}}_{j+1}\mathit{\theta }+{\displaystyle {\sum }_{i=1}^j{\mathcal{D}}_{i,j+1}{\widehat{\mathit{\alpha }}}_i({\mathit{\xi }}_i).}\hfill \end{array}$$(33)

From this equation, one sees immediately that the terms ∝θ cancel each other. Second, the two terms $\propto {\widehat{\mathit{\alpha }}}_j$ add up to zero. Third, also the sum of the two terms $\propto {\widehat{\mathit{\alpha }}}_{j-1}$ is zero. Thus, what remains to be shown is that the prefactor of the terms $\propto {\widehat{\mathit{\alpha }}}_i$,

$$\begin{array}{cc}\hfill K_i=& {\mathcal{D}}_{i,j+1}-{\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}{\mathcal{D}}_{i,j}\hfill \\ \hfill & -({\mathcal{D}}_{j+1}-{\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}{\mathcal{D}}_j){\mathcal{D}}_{j-1}^{-1}{\mathcal{D}}_{i,j-1},\hfill \end{array}$$(34)

for i ≤ j − 2 vanish. For this, we consider again Eq. (25), setting r = j − 1, s = j, t = j + 1, once with q = 0, and once with q = i. This then yields

$$\begin{array}{cc}\hfill \frac{1+z_{j-1}}{1+z_j}{\mathcal{D}}_{j,j+1}{\left({\mathcal{D}}_{j-1,j}^{\mathrm{t}}\right)}^{-1}& =({\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}{\mathcal{D}}_j-{\mathcal{D}}_{j+1}){\mathcal{D}}_{j-1}^{-1}\hfill \\ \hfill & =({\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}{\mathcal{D}}_{i,j}-{\mathcal{D}}_{i,j+1}){\mathcal{D}}_{i,j-1}^{-1}.\hfill \end{array}$$(35)

After multiplying by 𝒟_i,j−1, we see that the final equality shows that K_i = 0, which proves that Δ = 0, and thus the equivalence of the two forms of the lens equation. The equation K_i = 0 itself provides an interesting relation between distance matrices.

4. Mass-sheet transformation

In standard gravitational lensing, with a single deflector between the source and observer, there is a transformation of the mass distribution of the lens which keeps most observables invariant, the mass-sheet transformation (MST, see Falco et al. 1985). Since this transformation is accompanied by a uniform isotropic scaling in the source plane, all magnifications are scaled by the same factor, so that magnification (and thus observable flux) ratios are unchanged. The MST changes the product of time delay and Hubble constant, though, and the corresponding degeneracy can thus be broken by measuring the time delay in lens systems, assuming the Hubble constant to be known from other cosmological observations (see Schneider & Sluse 2013, and references therein). S14 has recently shown that a MST also exists in the case of two lens planes and two source planes. In this section, we show that also for perturbed gravitational lens systems, as considered in this paper, such a MST does exist.

4.1. Single main lens plane

We start with the case of a single lens plane, using the lens Eq. (14), and modify the deflection angle ${\widehat{\mathit{\alpha }}}_1({\mathit{\xi }}_1)$ to the new form

$$\begin{array}{c}\hfill {\widehat{\mathit{\alpha }}}_1^{\prime }({\mathit{\xi }}_1)=\lambda {\widehat{\mathit{\alpha }}}_1({\mathit{\xi }}_1)+{\mathsf{G}}_1{\mathit{\xi }}_1,\end{array}$$(36)

where λ is a real number, and G₁ is a matrix³. Throughout this section, a prime denotes a mass-sheet transformed quantity. Thus, the modified deflection angle is a scaled version of the original one, plus a term linear in the impact vector. If G₁ is symmetric, this linear term corresponds to a tidal matrix, i.e., adding a uniform mass sheet to the scaled lens mass distribution, plus an external shear. The modified lens equation then becomes

$$\begin{array}{c}\hfill {\mathit{\xi }}_2^{\prime }={\mathcal{D}}_2\mathit{\theta }-{\mathcal{D}}_{12}{\widehat{\mathit{\alpha }}}_1^{\prime }({\mathit{\xi }}_1)={\mathcal{D}}_2\mathit{\theta }-{\mathcal{D}}_{12}[\lambda {\widehat{\mathit{\alpha }}}_1({\mathit{\xi }}_1)+{\mathsf{G}}_1{\mathcal{D}}_1\mathit{\theta }].\end{array}$$(37)

As for the orignal MST, we require that the modified impact vector ${\xi }_2^{\prime }$ is related to the original one by a uniform, isotropic scaling, ${\xi }_2^{\prime }={\nu }_2{\xi }_2$, where ν₂ is the scaling factor. Thus we require

$$\begin{array}{c}\hfill {\mathcal{D}}_2\mathit{\theta }-{\mathcal{D}}_{12}[\lambda {\widehat{\mathit{\alpha }}}_1({\mathit{\xi }}_1)+{\mathsf{G}}_1{\mathcal{D}}_1\mathit{\theta }]={\nu }_2[{\mathcal{D}}_2\mathit{\theta }-{\mathcal{D}}_{12}{\widehat{\mathit{\alpha }}}_1({\mathit{\xi }}_1)].\end{array}$$(38)

In order to have the terms $\propto {\widehat{\mathit{\alpha }}}_1$ equal on both sides of Eq. (38), we need to set ν₂ = λ, as is also the case for the MST in standard lensing – the scaling of the source plane (here plane number 2) is the same as that of the deflection angle. The remaining terms are all ∝θ, and setting them equal on both side leads to 𝒟₂ − 𝒟₁₂𝖦₁𝒟₁ = λ𝒟₂, or

$$\begin{array}{c}\hfill {\mathsf{G}}_1=(1-\lambda ){\mathcal{D}}_{12}^{-1}{\mathcal{D}}_2{\mathcal{D}}_1^{-1}=(1-\lambda ){\mathsf{C}}_{12}.\end{array}$$(39)

Since C₁₂ is symmetric (see Schneider 1997), G₁ is indeed a tidal matrix. This single-main plane MST was also derived by McCully et al. (2014). Thus, in a generalized gravitational lens situation, the MST requires a shear in addition to a uniform mass sheet⁴.

4.2. Two main lens planes

We now consider a second lens plane at z₂, with the source plane being located at z₃. The modified lens equation then reads

$$\begin{array}{c}\hfill {\mathit{\xi }}_3^{\prime }={\mathcal{D}}_3\mathit{\theta }-{\mathcal{D}}_{13}{\widehat{\mathit{\alpha }}}_1^{\prime }({\mathit{\xi }}_1^{\prime })-{\mathcal{D}}_{23}{\widehat{\mathit{\alpha }}}_2^{\prime }({\mathit{\xi }}_2^{\prime }),\end{array}$$(40)

and we require the modified deflection angle ${\widehat{\mathit{\alpha }}}_2\prime$ to be chosen such that the 3-plane is just uniformly scaled relative to the original one, i.e., ${\xi }_3^{\prime }={\nu }_3{\xi }_3$. This condition then yields

$$\begin{array}{cc}\hfill {\mathcal{D}}_3\mathit{\theta }& -{\mathcal{D}}_{13}[\lambda {\widehat{\mathit{\alpha }}}_1({\mathit{\xi }}_1)+{\mathsf{G}}_1{\mathit{\xi }}_1]-{\mathcal{D}}_{23}{\widehat{\mathit{\alpha }}}_2^{\prime }({\mathit{\xi }}_2^{\prime })\hfill \\ \hfill & ={\nu }_3[{\mathcal{D}}_3\mathit{\theta }-{\mathcal{D}}_{13}{\widehat{\mathit{\alpha }}}_1({\mathit{\xi }}_1)-{\mathcal{D}}_{23}{\widehat{\mathit{\alpha }}}_2({\mathit{\xi }}_2)].\hfill \end{array}$$(41)

In order to account for the term $\propto {\widehat{\mathit{\alpha }}}_2$ on the r.h.s. of Eq. (41), the modified deflection has to be of the form⁵

$$\begin{array}{cc}\hfill {\widehat{\mathit{\alpha }}}_2^{\prime }({\mathit{\xi }}_2^{\prime })& ={\nu }_3{\widehat{\mathit{\alpha }}}_2({\mathit{\xi }}_2^{\prime }/\lambda )+{\mathsf{G}}_2{\mathit{\xi }}_2^{\prime }\hfill \\ \hfill & ={\nu }_3{\widehat{\mathit{\alpha }}}_2({\mathit{\xi }}_2^{\prime }/\lambda )+\lambda {\mathsf{G}}_2[{\mathcal{D}}_2\mathit{\theta }-{\mathcal{D}}_{12}{\widehat{\mathit{\alpha }}}_1({\mathit{\xi }}_1)].\hfill \end{array}$$(42)

This choice then yields equal terms $\propto {\widehat{\mathit{\alpha }}}_2$ on both sides. Equating the terms $\propto {\widehat{\mathit{\alpha }}}_1$ leads to the condition −λ𝒟₁₃ + λ𝒟₂₃𝖦₂𝒟₁₂ = −ν₃𝒟₁₃, or

$$\begin{array}{c}\hfill {\mathsf{G}}_2=(1-{\nu }_3/\lambda ){\mathcal{D}}_{23}^{-1}{\mathcal{D}}_{13}{\mathcal{D}}_{12}^{-1}.\end{array}$$(43)

Using the same arguments as in the Appendix of Schneider (1997), it is straightforward to show that any combination of distance matrices of the form

$$\begin{array}{c}\hfill {\mathcal{D}}_{\mathit{st}}^{-1}{\mathcal{D}}_{\mathit{rt}}{\mathcal{D}}_{\mathit{rs}}^{-1}\text{is symmetric,}\end{array}$$(44)

for 0 ≤ z_r <  z_s <  z_t. Hence, G₂ is symmetric, and thus corresponds to a tidal matrix. Equating the terms ∝θ in Eq. (41) then leads to

$$\begin{array}{c}\hfill (1-{\nu }_3){\mathcal{D}}_3=(1-\lambda ){\mathcal{D}}_{13}{\mathcal{D}}_{12}^{-1}{\mathcal{D}}_2+\lambda (1-{\nu }_3/\lambda ){\mathcal{D}}_{13}{\mathcal{D}}_{12}^{-1}{\mathcal{D}}_2,\end{array}$$(45)

which has the unique solution ν₃ = 1. Thus, as is the case for the standard multi-plane lens discussed in Schneider (2014), the MST does lead to no scaling in the plane j = 3. Therefore,

$$\begin{array}{c}\hfill {\mathsf{G}}_2=(1-1/\lambda ){\mathcal{D}}_{23}^{-1}{\mathcal{D}}_{13}{\mathcal{D}}_{12}^{-1}.\end{array}$$(46)

The implied scaling of the mass distribution in the plane i = 2 that follows from Eq. (42) is discussed in Schneider (2014); in short, the surface mass density distribution giving rise to ${\widehat{\mathit{\alpha }}}_2({\mathit{\theta }}_2)$ needs to be scaled in amplitude and scale-length to yield a deflection ${\widehat{\mathit{\alpha }}}_2(\lambda {\mathit{\theta }}_2)$.

4.3. Arbitrary number of planes

Here we generalize the MST to an arbitrary number of source/lens planes. It turns out that the lens equation in the form (31) is better suited for that purpose. We write it in the form

$$\begin{array}{c}\hfill {\mathit{\xi }}_{j+1}={\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}{\mathit{\xi }}_j+{\mathsf{B}}_j{\mathit{\xi }}_{j-1}-{\mathcal{D}}_{j,j+1}{\widehat{\mathit{\alpha }}}_j({\mathit{\xi }}_j),\end{array}$$(47)

where B_j is the term in parenthesis in Eq. (31). We now assume a scaling ${\xi }_j^{\prime }={\nu }_j{\xi }_j$ in every plane, and set the scaled deflection angles to be

$$\begin{array}{c}\hfill {\widehat{\mathit{\alpha }}}_j^{\prime }({\mathit{\xi }}_j^{\prime })={\nu }_{j+1}{\widehat{\mathit{\alpha }}}_j({\mathit{\xi }}_j^{\prime }/{\nu }_j)+{\mathsf{G}}_j{\mathit{\xi }}_j^{\prime }={\nu }_{j+1}{\widehat{\mathit{\alpha }}}_j({\mathit{\xi }}_j)+{\nu }_j{\mathsf{G}}_j{\mathit{\xi }}_j.\end{array}$$(48)

Note that Eqs. (36) and (42) are special cases of the relation (48) for j = 1, 2, respectively, and ν₁ = 1, ν₂ = λ, ν₃ = 1. Then, from ${\xi }_j^{\prime }={\nu }_j{\xi }_j$, we obtain from Eq. (47)

$$\begin{array}{cc}\hfill {\mathit{\xi }}_{j+1}^{\prime }=& {\nu }_j{\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}{\mathit{\xi }}_j+{\nu }_{j-1}{\mathsf{B}}_j{\mathit{\xi }}_{j-1}\hfill \\ \hfill & -{\nu }_{j+1}{\mathcal{D}}_{j,j+1}{\widehat{\mathit{\alpha }}}_j({\mathit{\xi }}_j)-{\nu }_j{\mathcal{D}}_{j,j+1}{\mathsf{G}}_j{\mathit{\xi }}_j\hfill \\ \hfill =& {\nu }_{j+1}[{\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}{\mathit{\xi }}_j+{\mathsf{B}}_j{\mathit{\xi }}_{j-1}-{\mathcal{D}}_{j,j+1}{\widehat{\mathit{\alpha }}}_j({\mathit{\xi }}_j)].\hfill \end{array}$$(49)

The terms $\propto {\widehat{\mathit{\alpha }}}_j$ cancel each other. The terms ∝ξ_j − 1 yield the condition ν_j + 1 = ν_j − 1, and equating the terms ∝ξ_j yields

$$\begin{array}{c}\hfill {\mathsf{G}}_j=(1-{\nu }_{j+1}/{\nu }_j){\mathcal{D}}_{j,j+1}^{-1}{\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1},\end{array}$$(50)

which is symmetric according to Eq. (44) and thus represents a tidal matrix. Thus, we obtain ν_j = λ for j even, and ν_j = 1 for j odd. Correspondingly,

$$\begin{array}{cc}\hfill {\mathsf{G}}_j& =(1-1/\lambda ){\mathcal{D}}_{j,j+1}^{-1}{\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}\mathrm{for}j\mathrm{even};\hfill \\ \hfill {\mathsf{G}}_j& =(1-\lambda ){\mathcal{D}}_{j,j+1}^{-1}{\mathcal{D}}_{j-1,j+1}{\mathcal{D}}_{j-1,j}^{-1}\mathrm{for}j\mathrm{odd}.\hfill \end{array}$$(51)

We note that Eqs. (39) and (43) are special cases of Eq. (51).

Hence we find that the MST in multiple (lens and source) plane gravitational lensing exhibits a curious behavior: The scaling factor in every second plane is just unity, whereas it is λ in the other half of the planes. In particular that means that a “standard candle” or “standard rod” in one of the planes with j odd cannot be used to break the degeneracy related to the MST, as the images of these sources are unaffected by the MST. The prefactor in the tidal matrices G_j are positive on every other plane, and negative on the remaining ones. If one disregards the perturbations between lens planes, so that the distance matrices 𝒟_ij reduce to angular diameter distances D_ij, then the G_j become scalars proportional to the density of uniform mass sheets; in this case, positive and negative densities of these sheets alternate.

4.4. Transformation of the time delay

We next consider how the MST affects the time delays. For that, we assume to have a source on plane number N + 1, with its light being deflected in N main lens planes. The corresponding LTTF is given in Eq. (19), where ξ_N + 1 is the position of the source in its source plane.

To obtain the corresponding function cτ′ after the MST, we first need to consider the transformation of the potential time delay. That is, we need to find the transformed deflection potential ${\widehat{\psi }}_i\prime ({\mathit{\xi }}_i\prime )$, which needs to satisfy ${\mathrm{\nabla }}_{{\mathit{\xi }}_i\prime }{\widehat{\psi }}_i\prime ({\mathit{\xi }}_i\prime )={\widehat{\mathit{\alpha }}}_i\prime ({\mathit{\xi }}_i\prime )$. For this, we make the ansatz ${\widehat{\psi }}_i\prime ({\mathit{\xi }}_i\prime )=a{\widehat{\psi }}_i({\mathit{\xi }}_i\prime /{\nu }_i)+{\mathit{\xi }}_i^{\prime \mathrm{t}}{\mathsf{G}}_i{\mathit{\xi }}_i\prime /2$. Taking the gradient yields ${\mathrm{\nabla }}_{{\mathit{\xi }}_i\prime }{\widehat{\psi }}_i\prime ({\mathit{\xi }}_i\prime )=(a/{\nu }_i){\widehat{\mathit{\alpha }}}_i({\mathit{\xi }}_i\prime /{\nu }_i)+{\mathsf{G}}_i{\mathit{\xi }}_i\prime$. This is seen to agree with ${\widehat{\mathit{\alpha }}}_i\prime ({\mathit{\xi }}_i\prime )$ in Eq. (48), provided a = ν_iν_i + 1 = λ. Thus,

$$\begin{array}{c}\hfill {\widehat{\psi }}_i^{\prime }({\mathit{\xi }}_i^{\prime })=\lambda {\widehat{\psi }}_i({\mathit{\xi }}_i^{\prime }/{\nu }_i)+\frac{1}{2}{\mathit{\xi }}_i^{\prime \mathrm{t}}{\mathsf{G}}_i{\mathit{\xi }}_i^{\prime }.\end{array}$$(52)

We then obtain for the transformed LTTF

$$\begin{array}{cc}\hfill c{\tau }^{\prime }=& {\displaystyle {\sum }_{i=1}^N(1+z_i)[\frac{1}{2}{({\nu }_i{\mathsf{C}}_{i,i+1}{\mathit{\xi }}_i-{\nu }_{i+1}{\mathcal{D}}_{i,i+1}^{-1}{\mathit{\xi }}_{i+1})}^{\mathrm{t}}{\mathsf{C}}_{i,i+1}^{-1}}\hfill \\ \hfill & \times ({\nu }_i{\mathsf{C}}_{i,i+1}{\mathit{\xi }}_i-{\nu }_{i+1}{\mathcal{D}}_{i,i+1}^{-1}{\mathit{\xi }}_{i+1})\hfill \\ \hfill & -\lambda {\widehat{\psi }}_i({\mathit{\xi }}_i)-\frac{{\nu }_i^2}{2}{\mathit{\xi }}_i^{\mathrm{t}}{\mathsf{G}}_i{\mathit{\xi }}_i].\hfill \end{array}$$(53)

We now show that

$$\begin{array}{c}\hfill {\tau }^{\prime }=\lambda \tau +F({\mathit{\xi }}_{N+1}),\end{array}$$(54)

which means that the transformed LTTF just scales by a factor λ, independent of the plane on which the source is located, plus a function which only depends on the location of the source, and thus cancels when considering time delays, i.e., differences between τ for pairs of multiple images of the source.

Taking the difference of τ′−λτ, we first note that the terms $\propto {\widehat{\psi }}_i$ drop out. Second, we note that the terms containing products of ξ_i and ξ_i + 1 also cancel, since ν_iν_i + 1 = λ. Thus we find that

$$\begin{array}{cc}\hfill c({\tau }^{\prime }-\lambda \tau )=& {\displaystyle {\sum }_{i=1}^N\frac{(1+z_i)}{2}\{{\mathit{\xi }}_i^{\mathrm{t}}[({\nu }_i^2-\lambda ){\mathsf{C}}_{i,i+1}-{\nu }_i^2{\mathsf{G}}_i]{\mathit{\xi }}_i}\hfill \\ \hfill & +({\nu }_{i+1}^2-\lambda ){\xi }_{i+1}^{\mathrm{t}}{\left({\mathcal{D}}_{i,i+1}^{-1}\right)}^{\mathrm{t}}{\mathsf{C}}_{i,i+1}^{-1}{\mathcal{D}}_{i,i+1}^{-1}{\mathit{\xi }}_{i+1}\}.\hfill \end{array}$$(55)

The term i = 1 of the first sum vanishes, since ν₁ = 1, and Eq. (39) holds. The final term (i = N) of the second sum depends only on the source position, and thus corresponds to the function F(ξ_N + 1) previously mentioned. Since this term is of no interest, we simply drop it from now on. Relabeling the index of the second sum as i → i − 1, we then get

$$\begin{array}{c}\hfill c({\tau }^{\prime }-\lambda \tau )={\mathit{\xi }}_i^{\mathrm{t}}\left[{\displaystyle {\sum }_{i=2}^N\frac{(1+z_i)}{2}({\nu }_i^2-\lambda ){\mathsf{K}}_i}\right]{\mathit{\xi }}_i,\end{array}$$(56)

where the matrices K_i are given as

$$\begin{array}{c}\hfill {\mathsf{K}}_i=\frac{1+z_{i-1}}{1+z_i}{\left({\mathcal{D}}_{i-1,i}^{-1}\right)}^{\mathrm{t}}{\mathsf{C}}_{i-1,i}^{-1}{\mathcal{D}}_{i-1,i}^{-1}+{\mathsf{C}}_{i,i+1}-\left(\frac{{\nu }_i^2}{{\nu }_i^2-\lambda }\right){\mathsf{G}}_i.\end{array}$$(57)

Since, according to the definition (17), ${\mathsf{C}}_{i-1,i}^{-1}{\mathcal{D}}_{i-1,i}^{-1}={\mathcal{D}}_{i-1}{\mathcal{D}}_i^{-1}$, we can rewrite K_i as

$$\begin{array}{c}\hfill {\mathsf{K}}_i=\frac{1+z_{i-1}}{1+z_i}{\left({\mathcal{D}}_{i-1,i}^{-1}\right)}^{\mathrm{t}}{\mathcal{D}}_{i-1}{\mathcal{D}}_i^{-1}+{\mathsf{C}}_{i,i+1}-\left(\frac{{\nu }_i^2}{{\nu }_i^2-\lambda }\right){\mathsf{G}}_i.\end{array}$$(58)

From Eq. (29),

$$\begin{array}{c}\hfill \frac{1+z_{i-1}}{1+z_i}{\left({\mathcal{D}}_{i-1,i}^{-1}\right)}^{\mathrm{t}}{\mathcal{D}}_{i-1}{\mathcal{D}}_i^{-1}={\mathcal{D}}_{i,i+1}^{-1}{\mathcal{D}}_{i-1,i+1}{\mathcal{D}}_{i-1,i}^{-1}-{\mathcal{D}}_{i,i+1}^{-1}{\mathcal{D}}_{i+1}{\mathcal{D}}_i^{-1}\end{array}$$(59)