{'fraction_non_alphanumeric': 0.04323544941310264, 'fraction_numerical': 0.012875187847772227, 'mean_word_length': 4.211450947190179, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gravitational lenses in arbitrary space-times

The precision reached by current and forthcoming strong-lensing observations requires to accurately model various perturbations to the main deflector. Hitherto, theoretical models have been developed to account for either cosmological line-of-sight perturbations, or isolated secondary lenses via the multi-plane lensing framework. This article proposes a general formalism to describe multiple lenses within an arbitrary space-time background. The lens equation, and the expressions of the amplification and time delays, are rigorously derived in that framework. Our results may be applied to a wide range of set-ups, from strong lensing in anisotropic cosmologies, to line-of-sight perturbations beyond the tidal regime.


Introduction
Gravitational lensing (Schneider et al., 1992) has become a key probe of the content, structure, and physical laws of our Universe. While weak lensing teaches us about the distribution of matter on cosmic scales (e.g. Asgari et al., 2021;Gatti et al., 2020), strong lensing lies amongst the best tools to measure today's cosmic expansion rate H 0 (Refsdal, 1964;Wong et al., 2020), and encapsulates valuable information on the small-scale distribution and nature of the dark matter (Despali et al., 2018;Diego et al., 2018;Harvey et al., 2020;Blum et al., 2020).
The quality of current imaging and data analysis requires an equally elaborate theoretical modelling of the strong lenses. That means complex models for the mass distribution of the main deflector, responsible for, e.g., the formation of multiple images of the same source, but also of the secondary deflectors which may perturb the effect of the main lens (Keeton et al., 1997;Tihhonova et al., 2018;Çağan Şengül et al., 2020;Li et al., 2020). Such perturbations are referred to as line-of-sight effects.
Two distinct frameworks have been developed to model line-of-sight perturbations in strong gravitational lensing. The first approach (Kovner, 1987;Seitz and Schneider, 1994;Bar-Kana, 1996;Schneider, 1997;Birrer et al., 2017) consists in treating secondary deflectors in the tidal regime, i.e., adding external convergence and shear to the main lens model. This technique is well suited to describe cosmological perturbations which would add to the effect of a lens; it has been successfully applied to measuring cosmic shear with Einstein rings by Kuhn et al. (2020). On the contrary, the second approach consists in modelling all the secondary deflectors as thin lenses. This formalism, called multi-plane lensing (Blandford and Narayan, 1986), is thus adapted to the description of isolated perturbers within an otherwise empty Universe, or within an ideal Friedmann-Lemaître-Robertson-Walker (FLRW) cosmology. The applicability of the original multi-plane formalism was then extended by the introduction of tidal planes (McCully et al., 2014), voids (McCully et al., 2017, or cosmological perturbations (Schneider, 2019); however, it has never been considered in full generality.
The present article proposes to fill this gap with an ultimate multi-plane formalism, where one or several lenses with arbitrary velocities may be placed in any smooth space-time background. This formalism encapsulates all the key lensing observables into a single versatile language. Our results may be applied to a wide range of set-ups; three specific examples will be provided in this article: (1) one lens with cosmological perturbations; (2) one lens in anisotropic cosmology, which recently attracted renewed attention (Akrami et al., 2020;Migkas et al., 2020;Secrest et al., 2021;Fosalba and Gaztañaga, 2020); and (3) multiple lenses in an under-dense Universe. Furthermore, this work establishes the necessary tools to accurately describe line-of-sight corrections in strong lensing beyond the standard convergence and shear (Fleury et al., 2019b(Fleury et al., , 2021.
The article is organised as follows. We start in section 2 with a discussion on the nature of the gravitational fields encountered by light beams, thereby defining the dichotomy between reference space-time and lenses. In section 3 we consider the case where a single lens is placed within an arbitrary reference space-time, before moving to the general N -lens case in section 4. We conclude in section 5.
We adopt units for which the speed of light is unity. Bold symbols (α, β, . . .) indicate two-dimensional Euclidean vectors, i.e., components of two-dimensional vectors over an orthonormal basis. Bold calligraphic symbols (A, D, . . .) refer to matrices, and cursive letters (P, F , L , . . .) to time-like or null geodesics.

Rays and beams of light
Let a source of light and an observer be placed in an arbitrary space-time. We call light beam the set of light rays that connect the source to the observer; a beam may have several components if the source is multiply imaged. We assume that the wavelength of light is much smaller than any relevant length scale of the problem (eikonal approximation) so that light rays are null geodesics.

Smooth and rough fields
As a beam of light propagates from the source to the observer, it may be deflected, distorted, and somehow split, by the gravitational field that it experiences. The formalism proposed in this article relies on the broad distinction between two categories of gravitational fields: smooth fields on the one hand, and rough fields on the other hand. In smooth-field regions, the light beam is continuously (and possibly strongly) distorted, but its integrity is preserved. In rough-field regions, the light beam experiences sudden deflections, which may give rise to multiple images. Thus, rough-field regions correspond to what is usually referred to as lenses, while smooth-field regions constitute a reference space-time, where light propagates from one lens to the next one. In the standard multi-plane lensing framework, the reference space-time is either Minkowski or FLRW. We shall not make this assumption here. 1 These considerations are illustrated in fig. 1. Specifically, we shall say that a gravitational field is smooth if the light beam can be considered infinitesimal within that field. In other words, a field is smooth if the corresponding space-time curvature is slowly varying and homogeneous within the light beam's cross section (Fleury et al., 2017(Fleury et al., , 2019a. Otherwise, we shall say that the field is rough. Let us illustrate this terminology with the example of a point mass M . The curvature that it produces at a distance r reads R µνρσ R µνρσ = 12(2GM/r 3 ) 2 ; hence the typical scale over which it changes appreciably is r. As a light beam with cross-sectional area A passes next to the mass, it experiences a smooth field if A r 2 , and a rough field otherwise. The picture gets slightly more complicated if we recall that the beam's cross-sectional area A(r) to depend on r, due to light focusing. Denoting A 0 the unlensed area of the beam, its lensed counterpart reads A(r) = µ(r)A 0 , where µ(r) ≡ (1 − r 4 E /r 4 ) −1 is the point-lens magnification, and r E its Einstein radius. The smooth-field condition then becomes A 0 r 2 (1 − r 4 E /r 4 ). Another example is light propagating through a diffuse distribution of matter, such as a gas cloud or a dark-matter halo. In that situation, space-time curvature is effectively dominated by its Ricci component (Fleury et al., 2017), which is mostly controlled by the matter density field ρ. Consequently, the corresponding complete space-time smooth reference lens = + rough-field region gravitational field is smooth if ρ is almost homogeneous on the scale of the beam's cross-section, and rough otherwise. Finally, we note that gravitational fields are not always directly generated by matter distributions in a Newtonian-like manner. The most immediate example is gravitational waves, which are nothing but propagating curvature. In our terminology, a gravitational wave is a smooth field if its wavelength is much larger than the beam's diameter, and rough otherwise. Perhaps even more relevant to lensing are the infinite-wavelength gravitational waves that characterise anisotropic cosmologies. For instance, the anisotropic expansion of Bianchi I models produces a homogeneous Weyl curvature, i.e. a smooth tidal field, which continuously shears and rotates light beams as they propagate .

Embedding lenses in the reference space-time
In the remainder of this article, we aim to concretely implement the distinction between smooth and rough fields into a generalised multi-plane lensing framework. But before doing so, let us briefly explain how one may technically treat the embedding of a lens within an arbitrary reference space-time.
Consider some rough-field region that is traversed by the light beam. In astrophysically relevant situations, the roughness of the field is due to a localised lumpy matter distribution (the lens), as depicted in fig. 2. Let L be the world-line of the lens's centre of mass. If we neglect the lens's self force, then L is a time-like geodesic of the space-time metricḡ µν generated by the rest of the universe, i.e., the reference space-time. Thus, we may introduce Fermi normal coordinates X α along L , which materialise the rest frame of the lens. With such coordinates, the reference metric is essentially Minkowskian across the rough-field region,ḡ αβ = η αβ , up to small corrections due to the local curvature of the reference space-time. These corrections are negligible in any astrophysically relevant situation, because the size of the region is on the order of the beam's cross section. Using the coordinates X α , one may then compute the gravitational field generated by the lens as if it were isolated, modulo the aforementioned small curvature corrections.
In the remainder of this article, we shall assume that the lenses are non-compact matter distributions with non-relativistic velocity dispersion. This notably excludes black-hole or neutron-star systems, and relativistic gases, but it provides an accurate description of any other astronomically relevant lens. In that context, we can treat rough gravitational fields as linear Newtonian perturbations, and hence apply the standard description of thin lenses (Schneider et al., 1992). Relaxing this assumption would affect the lens modelling and complicate the computation of time delays, but it would not change the essence of the framework that we propose in sections 3 and 4.

One lens in an arbitrary reference space-time
In this section, we tackle the case of a light beam that only propagates in smooth-field regions, except one localised rough-field region modelled as a thin lens. This problem has been studied by several authors in order to evaluate the impact of cosmological perturbations on the properties of a strong lens. Kovner (1987) refers to it as a thick lens in the telescope approximation; Schneider et al. (1992); Seitz and Schneider (1994) calls it generalised quadrupole lens; Bar-Kana (1996) writes about lensing with large-scale structure; and Schneider (1997), whose presentation is the closest to ours, simply calls it the cosmological lens. The corresponding formalism has been recently applied to the weak lensing of Einstein rings by Birrer et al. (2017), with the perspective of novel synergies between weak-lensing and strong-lensing observations (Kuhn et al., 2020).
The results derived in this section include, or are formally equivalent to, the aforementioned works'. However, we insert them in a broader context and extend their range of application -a novel example is proposed in section 3.6. Furthermore, to our knowledge we propose in section 3.4 the first rigorous proof of the expression of the time delay for lensing within an arbitrary reference space-time.

Geometric set-up
Let us describe in detail the geometry of the problem. We shall start with a four-dimensional picture, which forms the rigorous basis of the subsequent threedimensional picture, which in turn is more adapted to practical calculations. The four-dimensional discussion may be skipped by any reader who is not particularly interested in fully accurate definitions. (fig. 3). Let S , L , O be the world-lines of the source, lens and observer respectively. Let S ∈ S and O ∈ O be the emission and observation events; we call physical ray P a null geodesic connecting S to O. As we model the rough-field region as a thin lens, we will treat the physical ray as a set of two geodesics of the reference space-time, which undergoes a sudden deflection and a pause in the vicinity of the lens. The unlensed ray U is the null geodesic of the reference space-time that connects S to O. The intersection event O is earlier than O; their separation is the lensing time delay ∆t.

Four-dimensional picture
We introduce a fiducial ray F defined as the geodesic of the reference space-time that connects O to the lens's world-line L . 2 We denote with L ≡ F ∩ L the intersection event. Importantly, we assume that the rays P, U , F are all very The solid thick line indicates the physical ray P, the dashed line is the unlensed ray U , the dotted line is the fiducial ray F , and the thin solid line is the continuation C of the physical ray without deflection beyond the lens plane. The two grey shaded regions represent the lens and source world-sheets Σ d , Σ s ; their sections (thick grey lines) orthogonal to L or S indicate the lens and source planes at different times.
close to each other, so that any of their respective separations can be treated as infinitesimal in the reference space-time.
At L we define the lens plane as the two-dimensional space that is orthogonal to both L and F . In other words, the lens plane is strictly spatial (made of simultaneous events) in the lens's rest frame, and it is orthogonal to the spatial direction of the fiducial ray in that frame. The lens plane is well defined in the immediate vicinity of L, i.e., up to distances that are much smaller than the curvature radii of the reference space-time. So far we have defined the lens plane at the time of the event L; we may then generalise it to other times through parallel transport along L . Physically speaking, it means that the lens plane is non rotating. The three-dimensional time-like space that is swept by the lens plane as time goes on will be referred to as the lens world-sheet Σ d .
As the physical ray passes near the lens, it is effectively slowed down and deflected by its gravitational potential. In the thin-lens model adopted here, the deflection and delay are instantaneous. In other words, everything happens as if the photon were pausing and suddenly changing direction as it intersects Σ d . This pause is of course an idealised modelling of an otherwise continuous process, but it is an integral part of the thin lens approximation. We call P and P the beginning and end of the pause, respectively. Since the pause happens at a fixed position in the lens's frame, the separation between P and P is parallel to L . The duration of the pause represents the so-called potential time delay; 3 it generally differs from the final time delay ∆t measured along O, not only because these intervals are expressed in different rest frames, but also because ∆t contains additional, geometrical contributions.
We call x the (common) position of P and P in the lens plane, with respect to the origin set by L . 4 Similarly, we call r the position, in the lens plane, of the intersection R ≡ U ∩ Σ d between the unlensed ray and the lens world-sheet.
Let us now describe what happens in the vicinity of the source. Just like we defined the lens world-sheet Σ d from F and L , we define the source world-sheet Σ s from F and S . 5 Let us call F ≡ F ∩ Σ s the intersection of the fiducial ray and the source world-sheet. The line parallel to S and passing through F will be taken as the origin of the source plane. We call s the position of S in the source plane with respect to that fiducial origin.
Finally, let C be the continuation of the physical ray without deflection nor delay in the lens plane. We call I ≡ C ∩ Σ s its intersection with the source world-sheet. This event may be though of as the image event, i.e. the event that would be observed at O in the same direction as the physical ray in the absence of the lens. We denote with i the position of I in the source plane.
The above discussion shows that, as long as the various rays can be considered infinitesimally close to each other in the reference space-time, we can univocally define the notions of lens plane and source plane, and the various position vectors x, r, s, i in those planes. We shall now safely proceed with a three-dimensional representation of the problem, which will allow us to represent angles more easily. Figure 4, which is a spatial representation of the space-time diagram of fig. 3, corresponds to the more traditional way of picturing gravitational lensing by a thin deflector. The fiducial ray F , which is a null geodesic of the reference space-time, was chosen to go through the lens's position for simplicity, but any other nearby ray would be eligible. The lens plane (resp. source plane) is perpendicular to the fiducial ray in the lens's (resp. source's) rest frame. The position of the source S with respect to the fiducial origin F of the source plane is s.

Three-dimensional picture.
The unlensed ray U is the geodesic of the reference space-time connecting O to the S; we call β its angular separation with respect to the fiducial ray. Hence, β represents the direction in which the source would be observed in the reference space-time, without the lens. The unlensed ray intersects the lens plane at R, whose position with respect to L is r.
The physical ray P is made of two portions of geodesics of the reference spacetime, which intersect at P in the lens plane. The position x of P with respect to L is where the photon pauses and is deflected. We denote withα the deflection angle of the physical ray at P . Importantly,α is defined in the rest frame of the lens; it is thus subject to aberration effects when evaluated in another frame.
We denote with θ the separation between the physical ray and the fiducial ray at O, i.e. the observed direction of the image in the presence of the lens. The difference α ≡ θ − β, not to be confused withα, may be called displacement angle. The analogue of α in the source plane, i.e. the difference between the emission directions of the physical and unlensed rays, is denoted σ.
Finally, the continued ray C is the null geodesic of the reference space-time that coincides with the physical ray between O and P . As such, it is not deflected at P and intersects the source plane at a different point I, called image position. It represents the position of a source that would be observed in the direction θ in the absence of the lens. We call i the vector connecting F to I in the source plane.

Lens equation for one lens
Now that the geometric set-up has been fully described, we are ready to derive the lens equation for one thin deflector embedded in the arbitrary reference space-time.

Light propagation in the reference space-time.
In the reference space-time, by definition, the rays F , U , P, C are considered infinitesimally close to each other. Thus, the relative behaviour of any two of these rays can be described by the Sachs framework (Sachs, 1961), which is based on geodesic deviation. In what follows, we shall use a number of results of this formalism without deriving them; we refer the interested reader to, e.g., Fleury (2015, Chapt. 2) for further details. Let λ be a past-oriented affine parameter along the fiducial ray F , such that λ = 0 at O. Just like we defined the lens plane and source plane, we may introduce a local screen space at any point λ of F . Let us denote with x(λ) the screenspace separation between the physical and fiducial rays at λ. This vector field thus interpolates between 0 at O, x in the lens plane and s in the source plane.
In the reference space-time, x(λ) satisfies the differential equation where R is a particular screen-space projection of the Riemann curvature tensor of the reference space-time, called the optical tidal matrix. In fact, because eq. (1) is linear, we immediately conclude that it would equally apply to the separation of any two rays in the reference space-time. (1) is a linear second-order differential equation, so any of its solution is linearly related to its initial conditions. If this initial condition is considered at O where x(0) = 0, then there exists then there exists a 2 × 2 matrix D, called Jacobi matrix, such that x(λ) = D(λ)ẋ(0), where a dot denotes a derivative with respect to λ. The affine-parameter derivativeẋ(0) of x at O is related to the angle θ measured in the observer's rest frame asẋ(0) = ω 0 θ, with ω 0 the cyclic frequency of light in that frame. 6 More generally, if any two rays of the reference space-time emerge from, or converge to, a point (a) with angular separation θ a , then their transverse separation x b at another point (b) reads

Jacobi matrices. Equation
where ω a is the fiducial photon's cyclic frequency as measured in the same rest frame where θ a was defined. The Jacobi matrix D a b and the product ω a θ a are frame-independent. The presence of ω a is thus essential to account for aberration effects in θ a . We stress that a,b are not indices; they represent the positions where x, θ are evaluated. The non-standard notation "a b" in D a b is designed to clearly indicate that the two rays merge at (a). If the roles of (a) and (b) were swapped, i.e. if we considered two rays merging at (b) instead of (a), then their separation at (a) would read where the T superscript indicates the matrix-transpose operation. Equation (3) is known as Etherington's reciprocity law (Etherington, 1933); see Fleury (2015, § 2.2.3) for a derivation using the same conventions as this article.

Lens equation.
From the definition of the Jacobi matrix, we can express the position of the image i with respect to the source s in two different ways, where o, d, s respectively refer to the observer, deflector (or lens), and source positions. The deflection angleα(x) depends on the position x where the physical ray pierces the lens plane. For a thin lens made of non-relativistic and non-compact matter,α is dictated by the surface mass density Σ(x) in the lens plane (Schneider et al., 1992), where | . . . | denotes the Euclidean norm. The deflection angle can also be expressed as the gradient ofψ, which is twice the projected gravitational potential of the lens, The latter indeed satisfies the projected Poisson equation ∆ψ = 8πGΣ(x), where ∆ denotes the two-dimensional Laplacian and G is Newton's constant. Since x = D o d ω o θ, we conclude from eq. (4) that the lens equation reads with z d the observed redshift of the lens, 1 + z d = ω d /ω o . (7) is fully general, provided that light indeed encounters only one rough-field region, which can be modelled as a thin lens. No assumption is made on the reference space-time apart from the smoothness of its curvature. Hence, the deflector's redshift z d must not be understood as a cosmological redshift in general. Sinceα represents the deflection angle in the rest frame of the lens, it is by definition independent of the lens' motion. In that context, the redshift term 1 + z d encodes aberration effects. For instance, if the lens recedes from the observer, then its redshift increases, and the net observed deflection (1 + z d )α increases as well.

Important remarks. Equation
Let us finally stress that β is fundamentally linked to the reference space-time. Indeed, β represents the direction in which one would observe the source without the lens, i.e., if light were only propagating in smooth-field regions. In particular, it does not represent the direction where the source would be seen in an empty Universe, because the smooth-field regions do affect light propagation. If one wishes to work with another reference space-time, which does not correspond to the geometry where light propagates from one lens to another (e.g. Minkowski or FLRW), then one must explicitly allow for the difference between β and the new notion of unlensed direction. This will be illustrated with two concrete examples in sections 3.5 and 3.6 source plane σ dt s ds wave-fronts Figure 5. Two signals, emitted simultaneously with an angle σ from s + ds, are equivalent to two signals emitted with a relative delay dt s = σ · ds from s. Here, everything happens as if the physical signal (solid) was emitted slightly before the unlensed signal (dashed), so that dt s < 0 consistently with the opposite orientations of σ and ds.

Amplification matrix for one lens
From the lens equation (7), we immediately derive the amplification matrix of the system, i.e. the Jacobian matrix of the lens mapping θ → β(θ), whereĤ ij ≡ ∂ 2ψ /∂x i ∂x j is the Hessian matrix of the deflection potential. Note that, contrary to the latter, A is generally not symmetric, due to the coupling between the lens and the reference space-time.

Time delay for one lens
The time delay between the images of strong-lensing systems is a key observable in astronomy and cosmology (Wong et al., 2020). Its expression in the presence of cosmological perturbations can be found in Schneider et al. (1992); Bar-Kana (1996); Schneider (1997), but without a direct proof. Here we propose a rigorous derivation of the time-delay formula, inspired from the wave-front method introduced by Schneider et al. (1992, § 5.3). Let ∆t denote the time delay between the physical signal and the unlensed signal. In other terms, if a signal emitted by the source reached the observer at t 0 in the absence of the lens, then it would reach it at t = t 0 + ∆t in the presence of the lens. Note that ∆t is not directly observable, because it involves two signals that propagate in different space-times. It is, however, a convenient theoretical intermediate.
The time delay is conveniently parameterised as ∆t(θ, β), because it generally depends on the source position s = ω o D o s β, and on the point x = ω o D o d θ where the physical ray pierces the lens plane. In terms of that parameterisation, the observable time delay between two images A and B of the same source reads (θ B , β).
3.4.1 Differential time delay. The time delay between two signals emitted simultaneously in different directions depends on the source's position s. Indeed, as depicted in fig. 5, if the source lies at s + ds, then everything happens as if the two signals were emitted from s but with a slight relative delay dt s = σ · ds. Thus, if two signals emitted from s are observed with a time delay ∆t, then shifting the source by ds results in an additional delay where the redshift factor 1 + z s = dt o /dt s allows for the fact that dt s and d∆t were defined in distinct frames.
3.4.2 Time-delay formula. From eq. (9), we see that the expression of the time delay ∆t may be obtained by a line integral of the vector field σ(s), between an arbitrary reference point and s. In general, the line integral of a vector field depends on the path along which the integral is performed; except if the vector field is a gradient, which turns out to be the case. Namely, the emission angle σ between the lensed and unlensed rays depends on the source position s as where the scalar function T reads In eq. (10), the derivative dT /ds must be understood as a total derivative, in the sense that it accounts for the variation of both natural variables θ, β of T under a variation of s. A detailed proof of eq. (10) is provided in appendix A. Combining it with eq. (9), we immediately conclude that d∆t = dT , that is 7 ∆t(θ, β) = T (θ, β) + cst .
Therefore, the time delay between two different images A, B of the same source reads ∆t AB (β) = T (θ A , β) − T (θ B , β). Note that any function of β could be added to the expression of T (θ, β) without changing the observable ∆t AB . The time-delay matrix T generalises the more common notion of time-delay distance to an arbitrary reference space-time. Indeed, if the latter is chosen as the FLRW space-time, then T =τ 1, whereτ ≡ (1 + z d )D o dDo s /D d s is usually called the time-delay distance. Contrary to the Jacobi matrices that compose it, the time-delay matrix is symmetric, The time-delay matrix is related to, but different from, the telescope matrix introduced by Kovner (1987), and then used by Schneider et al. (1992); Seitz and Schneider (1994); Schneider (1997) within the generalised quadrupole lens model. Kovner (1987); Seitz and Schneider (1994); Schneider (1997) proved its symmetry with purely algebraic arguments. In appendix B, instead, we propose a geometric proof relying on Etherington's reciprocity law (3).

Fermat's principle.
Just like time-like and space-like geodesics can be defined from a stationary-time and stationary-length principle, null geodesics may be defined from Fermat's principle; see e.g. Schneider et al. (1992). In the present context, Fermat's principle states that the arrival time of a physical ray is stationary with respect to small variations of the position x where it pierces the lens plane. In other words, all things being fixed (notably s, r), physical rays must satisfy ∂T /∂x = 0. This property can be checked explicitly as follows. We first rewrite the first term of T as α · T α = (1 + z d )α ·Tα, whereT = ω d D o d D −1 o s D d s is also a symmetric matrix. Then, using the identityTα = x − r we immediately find so that physical light rays are indeed those whose deflection angleα is dictated by the physics of the lens plane.

Example: one lens in a perturbed cosmological background
Suppose that the reference space-time can be treated as a weakly perturbed homogeneousisotropic FLRW model. The associated Jacobi matrix takes the form In eq. (16),D a b = (1 + z b ) −1 f K (χ b − χ a ) denotes the angular diameter distance of (b) measured from (a) in the FLRW space-time, χ being the radial comoving distance, and f K (χ) ≡ sin( √ Kχ)/ √ K, with K is the spatial-curvature parameter. The second quantity, is the amplification matrix due to cosmological perturbations about FLRW, i.e. caused by the large-scale matter inhomogeneities in the Universe. Its key components are the convergence κ a b and complex shear γ a b ; we have neglected the anti-symmetric part of A a b , which represents the rotation of light beams with respect to parallel transport, because it is of the order of γ 2 if γ 1 (Fleury, 2015, § 2.3.2). At linear order in the matter density contrast δ(η, χ, x) ≡ (ρ −ρ)/ρ, the convergence and shear read (Fleury et al., 2019a) where η 0 denotes today's conformal time, and ϕ is the polar angle of x = x(cos ϕ, sin ϕ).

Lens equation.
In the lens equation (7), β stands for the direction in which a source would be observed without the lens, but still in the presence of the weak cosmological perturbations. In the lensing literature, however, it is understandably more common to write the lens equation in terms of the completely unlensed source position. That direction, which we shall denoteβ = A o s β, would be the observed position of the source without both the strong lens and the cosmological perturbations, i.e. in an ideal FLRW Universe. In terms ofβ, eq. (7) reads Equation (20) had already been obtained in the literature under various forms; for instance, it is equivalent to eq. (6.7) of Kovner (1987), eq. (14) of Bar-Kana (1996), eq. (35) of McCully et al. (2014); it also coincides with eq. (3.3) of Birrer et al. (2017) in the critical-mass-sheet approximation defined therein.

Equivalent lens.
In the absence of cosmological perturbations, i.e. for A a b = 1, the lens equation would take the form In the presence of perturbations, the lens equation (20) is formally equivalent to eq. (21), if one replacesᾱ(θ) by which may be called the equivalent lens. This shows that, in full generality, a lens modelᾱ must be supplemented with 9 real parameters (3 κ a b and 3 complex γ a b ) to properly account for cosmological perturbations. Degeneracies between these parameters might occur ifᾱ enjoys symmetries. In general, α eq (θ) cannot be written as a gradient, which means that it does not derive from a potential. An alternative approach (Schneider, 1997) which circumvents this issue consists in first applying a transformation β →β ≡ A o d A −1 d s β to the source plane. The resulting equivalent lens then derives from a potential.

Amplification matrix.
Following the discussion of section 3.5.1, we shall consider amplifications with respect to the homogeneous Universe, rather than the amplification due to the sole lens. The corresponding "total" amplification matrix is defined as A tot = dβ/dθ = A o s A, and reads where A ≡ 1 − dᾱ/dθ is the amplification matrix that would characterise the lens in the absence of cosmological perturbations, i.e. in an ideal FLRW reference space-time.
Equation (25) confirms that line-of-sight perturbations do not only add to the effect of a lens, but they also modify the effect of the lens itself.

Time delays.
For a perturbed FLRW reference space-time, the general expression (11) of the time delay becomes at first order in the cosmological perturbations, where, on the one hand where we denoted δA a b ≡ A a b − 1 for short, gathers all the corrections due to cosmological perturbations. Taken at face value, the correction δT thereby induced is quite complex. However, for practical analyses of time-delay observations, these may be reduced to a single external convergence and shear. First, since the source positionβ is unknown and hence a free parameter in such analyses, it does not make any difference whether one considers β = A −1 o sβ instead. Second, if the lens modelψ(θ) is general enough, 9 then it may effectively account for the corrections due to A o d in the argument ofψ. In that context, the time-delay model that must be used reads with an "external" amplification matrix

Example: lensing in an anisotropic Universe
As a second illustration, consider the case of a lens placed in a homogeneous but anisotropic Universe. If its homogeneity hyper-surfaces have no intrinsic curvature, ..
x N +1 = s p l a n e l = 1 p l a n e l p l a n e l + 1 s o u r c e p l a n e P U F Figure 6. Same as fig. 4, but with N lenses labelled by l. The transverse vector x l represents the position where the physical ray intersects the lth lens plane. The observer would correspond to l = 0 (x 0 = 0) and the source to l = N + 1 (x N +1 = s).
then it may be described by the Bianchi I geometry (Ellis and MacCallum, 1969). In the Bianchi I space-time, cosmic expansion is the same everywhere, but it is faster in some directions than in others. In comoving coordinates, its metric reads where a(η) is the volume scale factor, and the three β i (η), which sum to zero, encode the anisotropy of expansion. The propagation of light in Bianchi I cosmologies has been thoroughly investigated in Fleury et al. (2015Fleury et al. ( , 2016, thereby extending previous works on the same topic (Saunders, 1969). Let us restrict for simplicity to the weak-anisotropy limit (β i 1), where the optical properties of Bianchi I reduce to those of a perturbed FLRW whose amplification matrix reads We used the matrix where s o 1 , s o 2 denote the Sachs basis at the observer. Note that in A BI o s the source (s) and observer (o) are defined from their conformal time; this does not account for the change in redshift due to the anisotropic expansion, which reads 1 + z = (1 +z) All the results of section 3.5 can be used to describe lensing in a weakly anisotropic Universe, if one uses the A BI o s as amplification matrix all along.

N lenses in an arbitrary reference space-time
Let us now turn to the more involved case where light travels through an arbitrary number N of rough-field regions. A three-dimensional representation of the set-up is depicted in fig. 6, where the various quantities are defined in the same way as in the one-lens case. To our knowledge, such a situation had never been considered in full generality, although the hybrid framework proposed by McCully et al. (2014) may allow one to treat the most relevant cases in practice.

Lens equation for N lenses
For N lenses, the most direct derivation of the lens equation differs from the singlelens case in its structure. Let us introduce again a past-oriented affine parameter λ along F , such that λ = 0 at the observer and increases towards the source. We still denote with x(λ) the transverse separation between the physical and fiducial rays. Between two successive lens planes, x(λ) evolves smoothly according to the Sachs equation (1) of the reference space-time. When light reaches a lens plane, its sudden deflection implies a discontinuity of the derivativeẋ ≡ dx/dλ.

Deflection by a lens.
Let us denote ∆ẋ l ≡ẋ(λ + l ) −ẋ(λ − l ) the discontinuity oḟ x in the lth plane. This quantity is not exactly the deflection angleα l , because the latter represents the discontinuity in the way x changes with proper distance in the lens' rest frame. As shown in footnote 6, proper distance and affine parameter are related by d = ωdλ, so that where the minus sign comes from the conventional orientation ofα.

Between two lenses.
Contrary to section 3, the propagation from one lens plane to the next one does not necessarily involve merging rays; as such, it cannot be expressed in terms of Jacobi matrices only. However, we may still exploit the fact that x(λ) is the solution of a linear second-order differential equation (1). In full generality, the state of that solution at λ is encoded in the "phase space" vector Cauchy's theorem then implies the existence of a 4 × 4 Wronski matrix W such that for any λ 1 , λ 2 , X(λ 2 ) = W(λ 2 ← λ 1 )X(λ 1 ) .
Although the Wronski matrix will not appear in the final lens equation, it is a very convenient tool for its derivation. Note that the top-right 2 × 2 block of W(λ 2 ← λ 1 ) is actually the Jacobi matrix D 1 2 . Indeed, if two rays cross at λ 1 , then X 1 = (0,ẋ 1 ), and X 2 = (D 1 2ẋ1 ,ẋ 2 ). Another important property of W, which trivially follows from its definition (36), is the product law In what follows, we denote by W(l + 1 ← l) the Wronski matrix that relates X − l+1 ≡ X(λ − l+1 ) before its deflection in the (l + 1)th plane, to X + l ≡ X(λ + l ) after its deflection in the lth plane,

Recursion relation and lens equation. Since x(λ)
is continuous at λ l , the discontinuity of the phase-space vector X reads Denoting X l ≡ X − l (just before deflection at λ l ) for short, eqs. (38) and (39) yield the recursion relation X l+1 = W(l + 1 ← l)(X l + ∆X l ), which is solved as Finally, isolating the first two components, x l , of X l in the above, noting that X 0 = (0, ω o θ), and expressing the result in terms of β l ≡ (ω o D o l ) −1 x l , we find for any l, which only involves Jacobi matrices. The angle β l represents the direction in which a source at x l in the lth plane would be observed in the absence of foreground lenses m < l. The case l = N + 1, with β N +1 = β yields the explicit lens equation Note that the single-lens equation (7) is recovered for N = 1. Equation (41) matches eq. (10) of Schneider (2019), derived for a perturbed FLRW reference space-time.

Amplification matrix for N lenses
Just like the lens equation is a recursion relation, the amplification matrix for N lenses takes a recursive form. From β l we shall define the intermediate amplification matrix A l ≡ dβ l /dθ, which characterises the distortions of an infinitesimal source in the lth plane due to the foreground lenses m < l. By differentiating eq. (41), we find the recursion relation whereĤ m ij ≡ ∂ψ m /∂x i ∂x j , and with initial condition A 1 = 1 since β 1 = θ. The complete amplification matrix, accounting for all the lenses, is A ≡ dβ/dθ = A N +1 . single-lens time-delay formula (11) in each of these triangles. Precisely, in the lth triangle, two signals emitted simultaneously at x l+1 , the first one being undeflected (observed in the direction β l+1 ) and the other one deflected in the lth plane (observed in the direction β l ), would be received with a delay up to a constant, with for any l, m such that 0 < l < m, We note that the time-delay matrix satisfies the following unfolding relation; 10 ; for any three planes l, m, n such that 0 < l < m < n, See appendix C for a proof. Although we shall not use it here, eq. (46) is very useful to derive the expression of time delays in the context of multi-plane lensing with a dominant lens (Fleury et al., 2021). The delay between the actual ray (observed in the direction θ) and the undeflected ray (observed in the direction β) is then the sum of all these partial delays:

Rigorous derivation.
Although it yields the correct result, the above intuitive derivation of eq. (47) is actually incomplete. Indeed, we applied the single-lens time delay formula (11) to two non-physical rays. In particular, the "deflection angle" of the lth triangle, which is formed by the rays observed in the directions β l , β l+1 , and which we may denoteα l , is not the physical deflection angleα l = ∂ψ l /∂x l (see fig. 7). This is not a mere detail, because the equality between the deflection angle and the gradient of the lensing potential was a key of the derivation of eq. (11) proposed in appendix A. Therefore, nothing guarantees in principle that eq. (11) can be applied to rays that do not exhibit the correct deflection angle.
Fortunately, despite this weakness in the way that eq. (47) was obtained, the formula itself turns out to be correct. Let us prove this point. First of all, we note that the differential expression (9), i.e. d∆t = (1 + z s )σ · ds, still applies here because it relies on strictly local arguments in the source plane. Therefore, if we can show that (1 + z s )σ = dT /ds, then it would imply d∆t = dT just like in the single-lens case, and hence eq. (47) would follow.
Under a small variation of the source position s, all the intermediate positions x l (and thus β l ) change accordingly, so that each contribution T l of T varies as A similar calculation has already been performed in appendix A, except that (i) what was denoted s there is now x l+1 ; and (ii) the deflection angle must be replaced by the geometrical angleα l . In other words, eq. (73) applied to the configuration of the lth triangle reads from which it is easy to get the derivative with respect to s by multiplying both sides with (dx l+1 /ds) T ; the variable with respect to which one differentiates is just a matter of parameterisation here. The last step consists in noticing the identitỹ which clearly appears in fig. 7, and from which we finally get When summing the dT l /ds, all the terms cancel two by two, except the first one proportional to σ 0 ≡ 0, and the last one proportional to σ N ≡ σ. Therefore, which concludes our proof.

Fermat's principle.
Like in the single-lens case, Fermat's principle states that a light ray is physical if and only if the time-delay function is stationary for this ray.
Considering T as a function of x 1 , . . . , x N instead of β 1 , . . . , β N , one gets with similar calculations as in appendix A. Therefore, the function T is stationary with respect to changes of x 1 , . . . , x N for and only for the physical ray.

Example: lenses in an under-dense Universe
As an illustration of the framework developed in this section, consider the situation where a fraction f ∈ [0, 1] of the matter in the Universe is homogeneously distributed, while the rest is under the form of lenses. Such a scenario is comparable to the Einstein-Straus Swiss-cheese model (Einstein and Straus, 1945;Kantowski, 1969;Fleury et al., 2013), although in the latter the point-masses are surrounded by empty regions, while here we rather consider lenses placed within a homogeneous but under-dense cosmos. This under-dense Universe stands for our reference space-time, for which the Jacobi matrix is scalar, where D a b is given by the Kantowski-Dyer-Roeder distance (Kantowski, 1969;Dyer and Roeder, 1973;Dyer, 1973;Dyer and Roeder, 1974;Fleury, 2014) with smoothness parameter f . For z < 2, it may be approximated up to a few percent precision by the standard FLRW distance corrected by a negative convergence (Kainulainen and Marra, 2009) which was obtained from eq. (18) with δ = f − 1.
In that setup, the lens recursion (41) becomes which is identical to the original multi-plane recursion (Blandford and Narayan, 1986), up to the expression of the distances. Such an approach efficiently meets the somehow heavier formalism developed by McCully et al. (2017, § 2.2.2) to describe lenses separated by cosmic voids. Time delays are also affected by the under-density of the reference space-time; the associated function is identical to the one of the standard multi-plane framework, with except that the involved distances are changed, in particular

Conclusion
In this article, we have proposed a comprehensive and efficient framework to model gravitational lensing by one or several deflectors placed in an arbitrary reference spacetime. Our formalism relies on the dichotomy between smooth-field regions, which form our reference space-time where light beams can be considered infinitesimal, and rough-field regions which can be described as isolated thin lenses. In that context, we have derived the lens equations, and the expressions of the amplification matrix and time delays. We illustrated our general results to: a single lens with cosmological perturbations; a single lens in an anisotropic Universe; and to N lenses in an under-dense Universe. In our derivations of the lens equation, amplification matrix and time delays, we have assumed that the lenses could be individually described by their Newtonian gravitational potential. This assumption is well motivated in most astrophysically relevant situations, but it does not apply to very hot systems, or in the vicinity of compact objects. Dealing with such relativistic lenses would require to change that specific part of the modelling, but it would not affect how lenses are embedded in the smooth reference space-time.
Lensing by multiple deflectors had already been actively studied in the literature. The specific additions of the present work can be summarised as follows. In section 3, we extended the description of a single lens with cosmological perturbations (e.g. Schneider, 1997) to a lens within any reference space-time. In particular, we proposed in section 3.4 the first rigorous derivation of the expression of time delays in that general context. Section 4 further extended the results of section 3 to an arbitrary number of lenses, thereby generalising the standard multi-plane lensing formalism, which was hitherto limited to the Minkowski or FLRW reference space-times.
Our work establishes a firm basis for the description of gravitational lensing by several deflectors; in particular, it shall be applied to the accurate treatment of line-of-sight effects in strong gravitational lensing beyond the tidal approximation (Fleury et al., 2021, in prep.).  fig. 4, except that we added the vector y, which is obtained by continuing the physical ray without deflection from the source plane to the observer plane.
The set-up is reminded in fig. 8, where we added the useful quantity y = −ω d D d oα = ω s D s o σ in the observer plane.
The two natural variables of T being θ, β, we decompose the total derivative as dT ds = dθ ds Let us first compute ∂T /∂θ and ∂T /∂β. Thanks to the symmetry of T , derivatives of α · T α can be considered as hitting only the first α, while cancelling the pre-factor 1/2, so that To express T α, we note from fig. 8 = (1 + z s )σ , where the last two steps may be followed with fig. 8. We finally obtain dT ds = (1 + z d ) dx ds T α − dψ dx + (1 + z s )σ = (1 + z s )σ , which concludes the proof. We chose to explicitly keep the term proportional tô α − dψ/dx in the result in order to apply it more easily to the N -lens case.

B Symmetry of the time-delay matrix
Let us prove that the time-delay matrix is symmetric, i.e., T = T T . For that purpose, we shall express the vector y in the observer's plane (see fig. 8) in two different ways. The proof will only rely on Etherington's reciprocity law, D a b = −D T a b . A first option consists in following the chain y →α → i − s → α, invoking the corresponding Jacobi matrices. The computation goes as follows, Comparing eqs. (78) and (83), which are valid for any α, concludes the proof.

C Unfolding relation (46) for the time-delay matrix
Let us prove that, for any three lens planes l, m, n such that 0 < l < m < n, the time-delay matrix satisfies the unfolding relation This relation is illustrated in fig. 9, which is constructed as follows: (i) consider an arbitrary position y in the observer plane; (ii) connect that position to the origin N of the nth lens plane with a ray of the reference space-time, and call L, M the It follows that (T −1 ln − T −1 lm − T −1 mn )y = 0, and therefore eq. (84) because the above holds for any y in the observer plane.