Gravitational lensing theory
Theory of multi-plane lensing
Gravity is able to perform lensing analyses of lens systems composed of multiple lensing planes at different distances (redshifts). This capability is implemented in terms of distance coefficients and is based on the multi-plane lensing equations described below.
As discussed in the literature (see, e.g. Seitz & Schneider 1992 and Hilbert et al. 2009), multi-plane lensing equations can be written in an iterative way, essentially following backwards the light path. A light ray, during his path, will be subject to many deflections at its intersection with the various lensing planes. These deflections are often described in the literature in terms of the reduced deflection angle $\alpha_i = d_{is} \hat \alpha_i / d_s$, where $d_{is}$ (respectively $d_s$) is the distance between the $i$-th lens plane and the source (respectively, the observer and the source), and $\hat \alpha_i$ is the true deflection angle at the $i$-th lensing plane. The last quantity is actually the physical deflection, and depends solely on the mass distribution at the lens plane, with no dependence on any distance.
Since Gravity needs to deal with many lens planes and many sources at different distances, it is much more efficient to write multi-plane lens equations in terms of $\hat \alpha$ (rather than $\alpha$, as is done usually). To this purpose, let us call $\delta_i$ the vector denoting the intersection of the extension of the $i$ light path segment with the plane containing the observer. By definition, we have therefore $\delta_1 = 0$. Calling $x_i$ the angular position of the intersection of the light path with the $i$-th lensing plane, as seen from the observer, we have
\[\delta_{i+1} = \delta_i + d_i \hat \alpha_i(x_i) \; , \qquad x_{i+1} = x_i - \rho_i \delta_{i+1} \; .\]
Here $\rho_i = d_{i,i+1} / (d_i d_{i+1})$ is the ratio of the distance between the $i$-th plane and the next $i+1$-th plane and the distances of the $i$-th and $i+1$-th plane to the observer. Note that all distances here and below are taken to be comoving transverse cosmological distances (rather than angular-diameter distances as usually done: the difference is important for the computation of the time-delay below). For a flat cosmological model, the comoving transverse distances are identical to the comoving distance (line of sight), and therefore $d_{i,i+1} = d_{i+1} - d_i$.
The differential of previous equations gives the Jacobian matrix of the lensing equation (note that in case of multi-plane lensing the Jacobian matrix is often not symmetric). We define
\[\Delta_1 = \begin{pmatrix} 0 & 0 \\ 0 & 0 \end{pmatrix} \; , \qquad A_1 = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \; .\]
as the null and identity matrices. We then have
\[\Delta_{i+1} = \Delta_i + d_i \frac{\partial \hat \alpha_i(x_i)}{\partial x_i} A_i \; , \qquad A_{i+1} = A_i - \rho_i \Delta_{i+1} \; .\]
Finally, the time-delay expression can be computed as
\[T = \frac{1}{c} \sum_i \frac{1}{2} \rho_i \delta_{i+1}^2 - d_i \psi_i(x_i) \; ,\]
where $\psi_i(x_i)$ is the Fermat's potential of the $i$-th lens plane computed at the angular position $x_i$.