Recursion relations for gravitational lensing

The weak gravitational lensing formalism can be extended to the strong lensing regime by integrating a nonlinear version of the geodesic deviation equation. The resulting ‘roulette’ expansion generalises the notion of convergence, shear and flexion to arbitrary order. The independent coefficients of this expansion are screen space gradients of the optical tidal tensor which approximates to the usual lensing potential in the weak field limit. From lensed images, knowledge of the roulette coefficients can in principle be inverted to reconstruct the mass distribution of a lens. In this paper, we simplify the roulette expansion and derive a family of recursion relations between the various coefficients, generalising the Kaiser–Squires relations beyond the weak-lensing regime.

for an array of papers describing the implementation of Kaiser and Squires' work to obtain real data. (Note, however, the recent works by Fleury et al. [5][6][7], where the Kaiser-Squires theorem between shear and convergence is found to be violated due to the finiteness of the beam.) The first order correction to the weak-lensing convergence and shear has been calculated, and is called flexion [8][9][10]. More recently, higher order corrections have also been provided to arbitrary order in screen-space derivatives through the so-called 'Roulette' formalism [11,12]. This is derived by integrating the fully non-linear geodesic deviation equation. This formalism allows one to form an extended image and calculate its extended shape from the lensing potential at the centre of the image (or more generally from the optical tidal matrix).
The coefficients of the expansion depend on screen space derivatives of the lensing potential in the limit of a weak gravitational field, and neglecting time delay effects and so on. In this regime, we show here the formalism can be simplified in terms of the standard complex formulation of weak lensing. This then allows us to derive generalisations of the Kaiser-Squires relations to arbitrary order.

Mathematical preliminaries
There are two main routes to calculating the effect of a gravitational lens. One of them is to calculate all the null-geodesics converging at the point of observation. The other route is to start by the geodesic deviation equation (GDE), which typically is calculated to linear order. From a first order Taylor expansion in the Ricci rotation coefficients, the resulting equation contains one power of Riemann; ξ + R a kbk ξ b = O(ξ,ξ) 2 1st order GDE. (1) We follow the notation of [10], and the reader is referred to Fig. 1 therein for an intuitive explanation of the notation: A dot (˙) denotes derivative w.r.t. the null-geodesic 1 and k as an index represents the projection along a congruence of null-geodesics. More accurate information about the lens is obtained by starting from the Bazanski equation [13,14], which takes into account terms to second order in the deviation vector ξ and its derivativeξ . In the context of weak lensing, this was investigated in [10], and is the order that flexion appears. This was generalised to arbitrary order keeping only the maximum number of leading screen-space derivatives in [11,12], 2 where the Roulette formalism is derived. It was shown that if one only keeps leading order screen-space derivatives the equation to solve becomes (1) with the replacement The integral solution of (1) with the replacement (2) is written down by the solution procedure provided in [10]. The resulting integral solution is then expanded into maps for each order m of screen-space derivatives. It is shown how the general solution of the equation may be written as a sum over modes called 'roulettes'. The derivation is rather involved, and we refer the reader to [11] for a comprehensive overview, and [12] for the detailed derivation. The general mth order map is given as Eqs. (11)- (14) in [11]. The completely general roulette amplitudes α m s , β m s (even modes) andᾱ m s ,β m s (odd modes) are given as Eqs. (89)-(94) in [12]. With a symmetric amplification matrix, such as in the flat-sky approximation, all the odd modes vanish. In the present work we restrict attention to the weak-field approximation and neglect all but the highest number of screen space derivatives of the otential. Linearising around Minkowski space we write the perturbations with respect to Poisson gauge as and define a lensing potential We refer the reader to the discussion around Eq. 102 in [12] for further details.

Structure of the paper
The rest of the paper is structured as follows. In Sect. 2 we first write down recursion relations relating the roulette amplitudes to each other. Thereafter, the Roulette formalism (in the above given approximations) is cast into a compact form, using complex notation. We use this to write down a very neat expression for the resulting complex spin roulette amplitudes γ m s = α m s + β m s that invoke no sums or integrals. In Sect. 3 we formulate the recursion relations of the previous section in complex notation. In doing so, it becomes natural to introduce operators, and to write down operator identities. The action of the operator is to take in any one (complex) roulette amplitude and transform it to any other. Then, in Sect. 4 we compare our results with that of earlier works, and confirm agreement. In Sect. 5 we use the recursion relations to find general expressions for the lensing potential derivatives in terms of the spin roulette amplitudes α m s and β m s . Finally, in Sect. 6, we conclude.

Roulette modes in the weak-field approximation
For a general thin lens in the weak-field and flat-sky approximation the Roulette amplitudes [12] take the form where δ is the ordinary Kronecker delta 3 and χ is the distance from the observer to the lens. X , Y are coordinates in the lens-plane and ψ = ψ(X , Y ) is the lensing potential given by Eq. (4). Also the spin s is restricted sucht that 0 ≤ s ≤ m + 1 and the roulette amplitudes α m s , β m s may be non-zero only if m + s is odd. Finally; S m(k) The Eqs. (5) and (6) are horrendous, but it is easy enough to calculate α 0 1 and β 0 1 . We find By the recursion relations introduced below, these two 'lowest modes' of expansion will actually suffice. The recursion relations we find are: and Here (C + + ) and (C + − ) are algebraic coefficients given by The signs in the name of the coefficients (C + + ) and (C − + ) are there to indicate the direction of change in the value m (superscript) and s (subscript). We provide a more compact notation for these relations (and the inverses) in Sect. 3. In order to do so, however, we must first introduce complex variables.

Complex notation
Define first the complex variables By use of Pascal's triangle we now, upon a bit of straight forward algebra, obtain from Eqs. (5) and (6) the expression where we have defined Define similarly the differentiation operator Furthermore, let * denote the complex conjugate, and define also = ∂ 2 To be explicit, the operator inverse −1 is defined such that ( Here f is a function. A particularly useful expression for γ m s may be found by use of the recursion relations written down above (and complexified in the next section). To keep the algebra in the main text at a minimum we only quote the result, and refer the reader to "Appendix E" for a derivation of this relation. Here a − = (m + 1 − s)/2 and Γ m s are numerical coefficients given by Note that the expression for γ m s given by (20) does not involve any sum or integral. Hence it is computationally much more economic then the foregoing expression (16), which involve a sum over integrals. Below we have written down the first few complex roulette amplitudes by use of Eq. (20). Mode γ m s is given in row s, column m, starting with s = 0, m = 0 in the upper left corner.

D-operators: introducing a simpler notation
We define in this section a set of raising and lowering operators that will serve as recursion relations (on operator form) between the different amplitudes. This language is very natural, and allows for computing amplitudes without invoking the whole mathematical machinery of the underlying theory. This is highly beneficial, since the theory is rather involved, even with the approximations adopted in this work [cf. the sums and integrals in Eqs. (5)- (8)]. In want of a better name we call the identity operator D, 4

and define
Dγ = γ identity element, The connection to the Roulette formalism is made by defining two basis-operators D + + and D + − from which all other operators may be constructed. We define basis-operators Here C + + and C + − are numerical factors given by (14). Note that C + + and C + − will depend on the particular value of m and s. The recursion relations are now given by applying the basis-operators to γ m s , requiring that they act such that D + + : γ m s → γ m+1 s+1 and D + − : γ m s → γ m+1 s−1 for all valid (m, s)-positions. In particular,

Recursion relations
The + and − signs are hence there to indicate whether we add or subtract to the number m (upper index) and s (lower index). Figure 1 illustrates how the operators propagate an amplitude to one of its neighbours. The above recursion relations are the same as the relations given in the previous section, Eqs. (10)- (13). Note that the coefficients are always evaluated at the end-point position (see the example below). A proof of the recursion relations can be found in "Appendix A". The two recursion relations are independent of each other and all non-zero modes α m s and β m s can therefore be obtained by these relations by starting from the lowest mode 5 γ 0 1 = α 0 1 +iβ 0 1 , which is known from Eq. (9).
We shall note in passing that any higher order operator D a+ b+ such that γ m+a s+b = D a+ b+ γ m s may be constructed from the basis-operators. We refer the interested reader to "Appendix C" for an explanation and a more detailed scrutiny of these operators. Here we continue instead with an example, to show that these operators work in a very intuitive manner. To illustrate, let us calculate γ 3 2 from γ 0 1 .
In the last equality we inserted for γ 0 1 from (20), just to show that we obtain the correct expression, as compared with γ 3 2 in (21). And again: Note that the coefficients are all calculated at the end-point of each step.

Inverse operators
We define inverses with natural notation in the following way. Recalling that D is the identity element we have and similarely This is a natural notation, as one thinks of the '+'s and '−'s as canceling against each other when they are on the same script level 6 . Indeed the notation is right in suggesting that for instance Also this structure is further explored in "Appendix C". By application to γ , we find explicit expressions for the inverses as follows. By (24) and (27) we find Endoving the inverse C + + with the same natural notation, C − − ≡ (C + + ) −1 (and similarly for C + − ) we hence find the corresponding operator identities.

Inverse basis-operators
Here we have invoked (19). The index form of C − − and C − + must be handled with care, however. The key lies in the fact that the inverse is taken at a different (m, s)position, as a consequence of the fact that the coefficients are always evaluated at the end-position of each step. This is easy to see when we write it out, as follows.
Comparing with (14), and recognizing that the inverse in the last step is just the functional inverse, we necessarily find The inverses of the recursion relations (25) are now given as Figure 1 illustrates the action of the basis-operators and their inverses.

Connection to existing literature
It is instructive at this point to compare with previous work on gravitational lensing.
Kaiser & Squires described a method for inverting the measured shear caused by the lens to obtain the surface mass density of a lens [1]. To do so they started with the lensing equation. Their results are very well described by Schneider, who extended their work in two papers [16,17]. Schneider finds a way to use only local data to reconstruct the lens-mass from the shear-field, and also points out the need for a generalized 'Kaiser-Squires' inversion procedrue, to account for stronger lensing effects. The standard procedure when starting from the lensing equation is to define the lensing potentialψ θ is a vector of the two angles parametrizing the lensing 2-sphere, andκ is the dimensionless surface mass density. We use a tilde (˜) to distinguish the notation from our own, since we have taken a different approach by starting from the GDE. Hence our lensing potential is defined differently (4). From the linearized lens mapping β = θ − ∇ψ one Here ,i refers to derivative w.r.t. θ i . From this, Kaiser and Squires obtained an inversion formula, yielding an expression for the surface mass density. Definingγ =γ 1 + iγ 2 one finds where D = (θ 2 1 −θ 2 2 +2iθ 1 θ 2 )/|θ | 4 . As mentioned, this formula gives the convergence (surface mass density) in terms of the shear-field generated by the lens.

Connection to roulette amplitudes
This section is just to show that we obtain similar expressions. Indeed, calculating the first few modes from (16) we find These equations have precisely the same structure as (38)-(40). 7 This means that the convergence-to-shear map that Kaiser and Squires have found, corresponds to the first few terms in the series expansion of roulette amplitudes. The higher order terms may in the same way be obtained from (16). For an explicit example where the mass (convergence) and the derivatives of the mass distribution is calculated, please refer to [12,Sec. IV], where a circularly symmetric lens in the weak-field and thin-lens approximations is considered. In that work, the expressions (42)-(44) were inverted in a Kaiser-Squires like fashion to obtain the mass of the lens (Eqs. 149-154 therein) and its derivatives in terms of the roulette modes. In this work, we do the same (but more generally), except we restrict attention to expressing the derivatives of the potential (and not the mass and its derivatives) in terms of the roulette modes. This is done in the next section. Before that however, a comment as to the relation between the variables of Eqs. (42)-(44) and Eqs. (38)-(40) is in place. Following [12,Eq. 104], one defines in the weakfield approximation with the Roulette expansion the amplification matrix A where, recall, the potential is given by Eq. (4). From the above expression we find that A is the same 8 as A in Eq. (37), and hence α 1 0 = − χ 2κ , α 1 2 = χ 2γ 1 and β 1 2 = χ 2γ 2 . 7 Note that the lensing potential ψ is in our case defined in Eq. (4). 8 on the level where one views A and A as operators acting on defined potentials ψ andψ, respectively.

Derivatives of the potential in terms of roulette modes
By Eqs. (5) and (6) we can, as mentioned, calculate all the higher order modes. To illustrate, we give in the following all the derivatives of the potential up to 4th order in terms of the roulette modes.
For brevity we used, in the above, notation such that ∂ n x ∂ m y ψ ≡ ψ nxmy , where n, m are positive integers.

Conclusion
In this work we have generalized the Kaiser-Squires relations to strong lenses. This was done deriving a set of recursion relations for gravitational lensing (in the flatsky, weak-field and thin-lens approximations) using the Roulette formalism. We have presented a simplified version of this formalism by formulating the theory in complex variables, and by introducing D-operators. We have found a very simple expression for the roulette amplitudes, invoking no sums or integrals (20). This is beneficial for computational purposes.
We expect that recursion relations similar to those found in this work, will also hold true in more general implementations of the Roulette expansions. One should therefore investigate if the same structures are seen in the absence of the flat-sky, weak-field and thin lens approximations used here. This, however, is left for future work.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Proof of recursion relations
In this appendix we derive the recursion relations (25), Take first the definition of γ m+1 s+1 , Eq. (16). According to it one finds that the term with ∂ m−k+1 x ∂ k y ψ is given by Similarly, starting from the recursion relation γ m+1 s+1 = (C + + ) m+1 s+1 ∂ c instead, one finds the corresponding terms Equating Eq. (48) with Eq. (49) and inserting for (C + + ) m+1 s+1 gives, after some straight forward algebraic manipulation, the equation Following the same procedure for γ m+1 s−1 similarely gives the equation The two relations (50) and (51) must hold for the recursion relations to be correct. In "Appendix B" they are shown to hold.

B.2 Additional relations
One may also note relations like which follow from manipulation with trigonometric identities. In this paper we omit the (straight forward) proofs, as we make no use of these relations. Take now the definition and recall that (a + b) n = n k=0 n k a n−k (ib) k . We therefore have the generic relation

C D-operators as a group
The D-operators have a very intuitive notation, and hence they are very easy to use. Now that we have the basis-operators D + − and D + + in store, we can construct any other D-operator by extending our formalism in the most natural of manners. With this definition the notation is endowed with an intuitive 'cancellation of signs' property. For instance D +++−++ +++−+ = D ++++ +++ = D 4+ 3+ . It is also straigth forward to show that the D-operators commute as operators applied to γ . For the basis-operators we find as expected. Note now some further properties. First of all, since D * γ = D * γ = 0 (where * ∈ {+, −}), and since any operator D a+ b+ can be decomposed into the two basis-operators D + − and D + + (and their inverses), alongside D * , we can conclude that any operator with an odd total number of signs will vanish.

Definition 2 (Set of proper D operators)
Consider the setD of all operators with an even total number of signs, as defined over allowed (m, s)-values. From the definition of C + + and C + − (with corresponding inverses) it is clear that all of these are non-zero. One must have that Straight forward algebraic manipulation now shows that n = (a + b)/2 and m = (a − b)/2.
Note that (M − S)/2 is always a natural number. 5. From here one may write down a coordinate expression. There will, however, be many paths connecting two γ s, and the expression will be path dependent. One must make sure to choose a path that contains only (m, s)-positions where γ is defined. Also remember to use the inverses as basis-operators if the exponents are negative 9 . With these precautions, the rest is straight forward. Here we suffice it to give an expression for the case where both exponents are non-negative. In order to ensure that we only go through valid (m, s)-positions we move along D + + first.
where A ± = (M ±S)/2 for brevity. Note the product sums, which follow as a consequence of iterating over successive (m, s) values as we move along. Very similar expressions will be found also for the cases where the exponents are negative. These expressions relate the γ s at any two valid positions in the (m, s)-configurationspace to each other.

E Derivation of the simpler expression for m s
From the last appendix, and recalling Eq. (9), it is evident that in order to express the modes in terms of the potential, all we need to do is to insert m = 0 and s = 1 in the last equation. Starting at the lowest rung, γ 1 0 , we will not need the inverses. It is cleanest now to use only the tuple (m, s) as variables. Letting therefore M → m and S → s − 1, Eq. (67) becomes where we have also inserted γ 0 1 = −χ∂ c ψ, and defined a ± = m + 1 ± s 2 , which implies m + 1 = a + + a − and s = a + − a − . (69) One may show that Hence we can always isolate the real and imaginary parts. Also note that