Diffraction of electromagnetic waves by an extended gravitational lens
aa r X i v : . [ g r- q c ] F e b Diffraction of electromagnetic waves by an extended gravitational lens
Slava G. Turyshev , Viktor T. Toth Jet Propulsion Laboratory, California Institute of Technology,4800 Oak Grove Drive, Pasadena, CA 91109-0899, USA and Ottawa, Ontario K1N 9H5, Canada (Dated: February 9, 2021)We continue our study of the optical properties of the solar gravitational lens (SGL). Taking thenext step beyond representing it as an idealized monopole, we now characterize the gravitationalfield of the Sun using an infinite series of multipole moments. We consider the propagation ofelectromagnetic (EM) waves in this gravitational field within the first post-Newtonian approximationof the general theory of relativity. The problem is formulated within the Mie diffraction theory. Wesolve Maxwell’s equations for the EM wave propagating in the background of a static gravitationalfield of an extended gravitating body, while accounting for multipole contributions. Using a wave-theoretical approach and the eikonal approximation, we find an exact closed form solution for theDebye potentials and determine the EM field at an image plane in the strong interference regionof the lens. The resulting EM field is characterized by a new diffraction integral. We study thissolution and show how the presence of multipoles affects the optical properties of the lens, resultingin distinct diffraction patterns. We identify the gravitational deflection angle with the individualcontributions due to each of the multipoles. Treating the Sun as an extended, axisymmetric, rotatingbody, we show that each zonal harmonics causes light to diffract into an area whose boundary isa caustic of a particular shape. The appearance of the caustics modifies the point-spread function(PSF) of the lens, thus affecting its optical properties. The new wave-theoretical solution allows thestudy gravitational lensing by a realistic lens that possesses an arbitrary number of gravitationalmultipoles. This angular eikonal method represents an improved treatment of realistic gravitationallensing. It may be used for a wave-optical description of many astrophysical lenses.
I. INTRODUCTION
Studied for over a century [1, 2], gravitational lensing today is well understood [3–5]. It occurs when light travelsin the vicinity of a gravitating body. In the post-Newtonian limit of the general theory of relativity, the gravitationalfield serves as a refracting medium [6, 7] that deflects light rays towards the body.Following a methodical approach, we began our investigation by treating the solar gravitational field as a sphericallysymmetric field of a gravitational monopole, or point mass [8, 9]. After passing by such a monopole, light rays arefocused in what we call the region of strong interference (Fig. 1), with impressive optical properties including significantlight amplification. However, even gravitational monopole lenses are subject to optical aberrations. As the deflectionangle is inversely proportional to the impact parameter, light rays with larger impact parameters with respect tothe lens are focused at larger distances from it. This causes spherical aberration, leading to blurred images and therequirement to employ appropriately designed deconvolution algorithms [10].With this model, we were able to establish the basic properties of the solar gravitational lens (SGL) and understandimage formation and image recovery. We considered gravitational lensing by the Sun as the means to obtain high-resolution images of faint objects, such a exoplanets. To enable practical applications of the SGL, we developed awave-optical treatment of the diffraction of light in the presence of the solar gravity field. We studied the impact ofthe solar corona on light propagation in the vicinity of the Sun. We showed that diffraction in the solar atmospheredefocuses EM waves for wavelengths greater than 1 mm, but its impact is negligible at optical and IR wavelengths[11, 12]. We extended our formulation to the case of extended sources at large but finite distances [13]. We studiedimage formation with the SGL [14, 15] and addressed the realistic sensitivity of prospective imaging observations [10].In addition, we studied the image recovery process and have shown that the SGL may be used for multipixel imagingof exoplanets [16] that may be conducted in the context of a realistic space mission [17].The next step is dictated by the realization that nothing is perfect in life, not even the Sun. Its rotation, theresulting oblateness, and its internal mass distribution result in a gravitational field that deviates from the idealizedmonopole. These deviations are small (in fact, the Sun is almost perfect), but given the distance and length scalesinvolved, their impact cannot be neglected. The dimensionless magnitude of corrections due to deviations from themonopole is of O (10 − ). This is not much until we consider that deflection of light by the SGL amounts to displacinga ray of light by at least as much as the solar radius by the time that ray approaches the focal region. An O (10 − )correction on this scale amounts to an additional deflection by tens if not hundreds of meters, which is quite significantwhen compared to the scale (typically, 1–10 km across) of the projected image size of a desired target.Therefore, it is necessary to develop a formalism to modify the point-spread function (PSF) of the SGL, taking !" !" FIG. 1: The different optical regions of the SGL (from [13]) with the strong interference region formed beyond 547.8 AU. into account that, on the one hand, we deal with very large distances (measured in light years for distant imagingtargets or in many hundreds of astronomical units (AU)) when it comes to the distance of the focal region from theSun) and, on the other hand, distances measured on the scale of meters or less (such as the telescope aperture or thecentimeter-scale Airy pattern that appears in the image plane, itself a result of observing a signal with a wavelength of O (1 µ m)). Consequently, even higher-order multipole moments (octupole, dodecapole, hexadecapole moments) of theSun may have to be considered for accurate image modeling and reconstruction of some exoplanetary targets. Thesemoments break the azimuthal symmetry of the PSF, introducing caustics instead of the regular Bessel J pattern [9].This is why we are turning the page on the chapter dealing with monopole gravitational lenses. With the presentpaper, we open a new, exciting area of investigation, aimed at developing a comprehensive description of realisticgravitational lenses possessing an arbitrary number of gravitational multipole moments.This paper is organized as follows: In Section II we discuss the solution of Maxwell’s equations in the curvedspacetime of the solar gravitational field, described at the first post-Newtonian approximation of the general theory ofrelativity. We develop a solution for the Debye potential using the eikonal approximation. In Section III, we formulatea generic solution for EM waves in the field of a static, extended gravitational lens. In Section IV we develop a generalsolution for the EM field, characterizing the scattering of EM waves on an extended lens. In Section V, we studythe EM field in the interference region. We develop a new integral formulation that describes light diffraction inthe strong interference region. In Section VI we discuss the results obtained and the next steps in our investigation.To aid with the flow of material in this paper, we placed some important derivations in appendices. Appendix Adiscusses an approach to Maxwell’s equations for EM waves propagating on the background of the static gravitationalfield of an extended lens. In Appendix B, we discuss the eikonal phase for i) a generic axisymmetric gravitating bodywhose gravitational potential given by a set of gravitational multipoles, and ii) generic spatial-trace free (STF) tensorsrepresenting bodies with arbitrary gravitational fields. Finally, in Appendix C, we explore the connection betweenour results and the path integral formalism. II. ELECTROMAGNETIC WAVES IN A STATIC GRAVITATIONAL FIELD
To describe the optical properties of the solar gravitational lens (SGL), we use a static harmonic metric in the firstpost-Newtonian approximation of the general theory of relativity. The line element for this metric may be given, inspherical coordinates ( r, θ, φ ), as [6, 18]: ds = u − c dt − u (cid:0) dr + r ( dθ + sin θdφ ) (cid:1) , (1) The notational conventions used in this paper are the same as in [7, 18]: Latin indices ( i, j, k, ... ) are spacetime indices that run from 0to 3. Greek indices α, β, ... are spatial indices that run from 1 to 3. In case of repeated indices in products, the Einstein summation ruleapplies: e.g., a m b m = P m =0 a m b m . Bold letters denote spatial (three-dimensional) vectors: e.g., a = ( a , a , a ) , b = ( b , b , b ). Thedot ( · ) and cross ( × ) are used to indicate the Euclidean inner product and cross product of spatial vectors; following the conventionof [6], these are enclosed in round and square brackets, respectively. Latin indices are raised and lowered using the metric g mn . TheMinkowski (flat) spacetime metric is given by γ mn = diag(1 , − , − , − γ µν a µ b ν = − ( a · b ). We use powers of the inverse ofthe speed of light, c − , and the gravitational constant, G as bookkeeping devices for order terms: in the low-velocity ( v ≪ c ), weak-field( r g /r = 2 GM/rc ≪
1) approximation, a quantity of O ( c − ) ≃ O ( G ), for instance, has a magnitude comparable to v /c or GM/c r .The notation O ( a k , b ℓ ) is used to indicate that the preceding expression is free of terms containing powers of a greater than or equal to k , and powers of b greater than or equal to ℓ . Other notations are explained in the paper. where, to the accuracy sufficient to describe light propagation in the solar system, the quantity u can be given interms of the Newtonian potential U as u = 1 + c − U + O ( c − ) , where U ( x ) = G Z ρ ( x ′ ) d x ′ | x − x ′ | . (2)The metric (1)–(2) allows us to consider effects on the propagation of light by the gravitational field of the Sun, dueto an arbitrary static gravitational field. We may also consider solar rotation, but its effect, although measurable, ismuch less than those of the solar monopole and quadrupole [19]. Furthermore, it was shown in [20] that to first orderin the gravitational constant G , a rotating and a nonrotating lens cannot be distinguished. Thus, solar rotation isnot included in our formalism.The gravitational field of the Sun is weak: its potential is GM/c r . × − everywhere in the solar system. Thisallows us to carry out calculations to the first post-Newtonian order, while dropping higher-order terms.We use the the generally covariant form of Maxwell’s equations for the electromagnetic (EM) field [7, 9] and considerthe propagation of an EM wave in the vacuum in the absence of charges and currents, i.e., j k = ( ρ, j ) = 0. As weshowed in [9, 11], for the metric (1) we obtain the following form for Maxwell’s equations:curl D = − u ∂ B c∂t + O ( G ) , div (cid:0) u D (cid:1) = O ( G ) , (3)curl B = u ∂ D c∂t + O ( G ) , div (cid:0) u B (cid:1) = O ( G ) , (4)where the differential operators curl F and div F are with respect to the usual 3-space Euclidean flat metric (see [9]for technical details). A. Representation of the EM field in terms of Debye potentials
To describe the problem of an EM wave propagating in the gravitational field of an extended lens that induces thestatic gravitational field with metric (1), we follow the Mie diffraction theory [21, 22] that allows one to determinethe three-dimensional structure of the EM field diffracted on a spherical obstruction. This technique is done basedrepresenting the Maxwell equations (3)–(4) in terms of the Debye potentials (see [9, 11] and references therein).Relying on the approach that we previously developed (see [9, 11]), in Appendix A we obtain the complete solutionof these equations in terms of the electric and magnetic Debye potentials [22], e Π and m Π. We follow closely thederivation in [9] (see Appendix E therein) and also in [11] (see Appendix A therein).We treat the lens as a compact gravitating body whose gravitational potential admits a representation in the formof an infinite series of zonal and tesseral harmonics (e.g., as given by (B1)). The result is a system of equations forthe components of a monochromatic EM field, characterized by the wavenumber k = ω/c :ˆ D r = 1 u n ∂ ∂r h r e Π u i + (cid:16) k u − u (cid:0) u (cid:1) ′′ (cid:17)h r e Π u io , (5)ˆ D θ = 1 u r ∂ (cid:0) r e Π (cid:1) ∂r∂θ + ikr sin θ ∂ (cid:0) r m Π (cid:1) ∂φ , (6)ˆ D φ = 1 u r sin θ ∂ (cid:0) r e Π (cid:1) ∂r∂φ − ikr ∂ (cid:0) r m Π (cid:1) ∂θ , (7)ˆ B r = 1 u n ∂ ∂r h r m Π u i + (cid:16) k u − u (cid:0) u (cid:1) ′′ (cid:17)h r m Π u io , (8)ˆ B θ = − ikr sin θ ∂ (cid:0) r e Π (cid:1) ∂φ + 1 u r ∂ (cid:0) r m Π (cid:1) ∂r∂θ , (9)ˆ B φ = ikr ∂ (cid:0) r e Π (cid:1) ∂θ + 1 u r sin θ ∂ (cid:0) r m Π (cid:1) ∂r∂φ , (10)where the electric and magnetic Debye potentials Π( r ) = ( e Π; m Π) satisfy the following wave equation: (cid:0) ∆ + k u (cid:1)h Π u i = O (cid:16) r g , J r Π u (cid:17) , (11)with the quantity u given by (2). The Newtonian potential U in (2) at this point is unconstrained and can describe anarbitrary (weak, static) gravitational field. Here, r g = 2 GM/c is the Schwarzschild radius of the lens; J characterizesthe quadrupole component of the gravitational potential, U , of an extended gravitational lens.Essentially the solution (5)–(10) together with (11) was obtained under the thin lens or eikonal approximationwhere the primary emphasis was on the effect of the gravitational field on the phase of the EM wave rather than itsamplitude. This approximation is well-justified as the source and the image plane are at very large distances fromthe lens. Our analysis showed that the effects of the higher order gravitational multipoles, starting from J , dependon the distance to the lens and, thus, may be neglected.As a result, the entire solution to Maxwell’s equations describing light propagation in the weak gravitational fieldwith the post-Newtonian metric tensor (1) depends on the solution of the wave equation for the Debye potential (11).Using the expression for u from (2), this equation is given as (see (A40)) (cid:16) ∆ + k (cid:0) Uc (cid:1)(cid:17)h Π u i = O (cid:16) r g , J r Π u (cid:17) . (12)Expressions (5)–(10) together (12) represent the solution of the Mie problem in terms of Debye potentials [21, 22], inthe presence of the gravitational field of an extended gravitating body, taken at the first post-Newtonian approximationof the general theory of relativity [8, 9] under the eikonal (or, essentially, the thin lens) approximation.The set of equations (5)–(10) with (12) determines the Debye potential for the entire problem. We see that thesolution of (12) now depends on the entire Newtonian potential, U ( r ) that may have arbitrary complexity. No exactsolution of this time-independent Schr¨odinger equation with exist. Thus, we need to develop an approximate solutionthat is suitable for our situation. We found an approach to develop such a solution using the eikonal approximation. B. Separating variables in the equation for the Debye potential
To consider the eikonal approximation, we present the Newtonian potential, U , as U ( r ) = GMr + δU ( r ) , (13)where the first term is the spherically symmetric monopole contribution and the second term, δU ( r ), represents thecombined contribution of all the other terms in a suitable expansion of U ( r ).If δU ( r ) is absent, (A40) reduces to the case of diffraction of the EM waves by a gravitational monopole (i.e.,Schr¨odinger’s equation with a Coulomb potential—see details in [9]): (cid:16) ∆ + k (cid:0) r g r (cid:1)(cid:17)h Π u i = 0 . (14)This equation describes light scattering that is dominated by a spherical relativistic potential due to a gravitationalmonopole (which is equivalent to an attractive Coulomb potential, discussed in quantum mechanics [23–25]). In ourcase, this equation describes the incident wave that travels towards the lens from the source.The solution to (14) is well known (see [9] for details). In this case, Eq. (12) is typically solved by separatingvariables [22], which, in spherical polar coordinates, takes the form [9, 12]:Π u = 1 r R ( r )Θ( θ )Φ( φ ) , (15)with integration constants and coefficients that are determined by boundary conditions. Direct substitution into (14)reveals that the functions R, Θ and Φ must satisfy the following ordinary differential equations: d Rdr + (cid:16) k (1 + 2 r g r ) − αr (cid:17) R = O ( r g ) , (16)1sin θ ddθ (cid:16) sin θ d Θ dθ (cid:17) + (cid:0) α − β sin θ (cid:1) Θ = O ( r g ) , (17) d Φ dφ + β Φ = O ( r g ) . (18)As we discussed in [9], the solution to (18) is given as usual [22, 24]:Φ m ( φ ) = e ± imφ → Φ m ( φ ) = a m cos( mφ ) + b m sin( mφ ) , (19)where β = m , m is an integer and a m and b m are integration constants.Equation (17) is well known for spherical harmonics. Single-valued solutions to this equation exist when α = l ( l + 1)with ( l > | m | , integer). With this condition, the solution to (17) becomesΘ lm ( θ ) = P ( m ) l (cos θ ) . (20)We now focus on the equation for the radial function (16), where, because of (17), we have α = ℓ ( ℓ + 1). As aresult, (16) takes the form d Rdr + (cid:16) k (1 + 2 r g r ) − ℓ ( ℓ + 1) r (cid:17) R = O ( r g ) . (21)The solution to this equation is given in the form of a Coulomb function F ℓ ( kr g , kr ) [9].Collecting results for Φ m ( φ ), Θ lm ( θ ) and R ℓ = F ℓ ( kr g , kr ), we can assemble the ultimate solution to (14), as wasdone in [9, 11]. This solution is used to describe the electric and magnetic potentials of the incident wave, e Π and m Π , which may be given in terms of a single potential Π ( r, θ ) (see [9] for details): (cid:18) e Π m Π (cid:19) = (cid:18) cos φ sin φ (cid:19) Π ( r, θ ) , where Π ( r, θ ) = E k ur ∞ X ℓ =1 i ℓ − ℓ + 1 ℓ ( ℓ + 1) e iσ ℓ F ℓ ( kr g , kr ) P (1) ℓ (cos θ ) + O ( r g ) . (22)Therefore, in the case when deviations from the monopole gravitational field represented by the term δU ( r ) in (13)are absent, we can find a solution for the Debye potential (14) by separating variables with the ansatz (15) that isused to deal with the Coulomb potential [8, 9]. The structure of the resulting solution (22) reflects the sphericalsymmetry that is preserved in this case. The presence of the monopole is manifested by the potential term 2 r g /r inthe equation for the radial function (16). Note that the other equations (17) and (18) are not affected by gravity.The situation changes drastically when the term δU ( r ) is present in (13). In this case, (12) becomes highly nonlinearand separation of variables (15) does not work. No exact solution of this equation is known. However, in some casesthis equation may be solved using well-justified approximation methods. One such method, the eikonal approximation,is particularly useful for high-energy atomic scattering [26–28] and it is applicable in our case, which corresponds tothe high-energy approximation in optical scattering [29–31]. C. Finding the Debye potential with the eikonal approximation
We may now use the result (22) as the basis to find solutions when U is not restricted to a monopole gravitationalfield. We extend the discussion in the preceding subsection by considering the complete post-Newtonian potential U of an extended body as the sum of two terms that includes the monopole field, GM/r , which is long-range [9], anddeviations from the monopole, which constitute a short-range potential, V sr ( r ) = δU ( r ) /c . This yields the followingform for the potential term in (A40): 4 Uc = 2 r g r + 4 V sr . (23)This decomposition allows us to proceed with solving (A40) that now takes takes the form (cid:16) ∆ + k (cid:0) r g r + 4 V sr ( r ) (cid:1)(cid:17) Π( r ) = O (cid:16) r g , J r Π (cid:17) , (24)where V sr is from (23). In explicit form this short-range potential is given by either by (B23) that is valid for anygeneric gravitational field, or expressed in terms of zonal harmonic coefficients J n using (B10), which is more suitableto describe the gravitational field of a rotating, axisymmetric mass, such as the Sun.To solve (24), we will treat V sr ( r ) as a perturbation to the monopole term and will use the eikonal approximation[7, 22, 30, 32, 33]. To implement this approach, we consider a trial solution in the formΠ( r ) = Π ( r ) φ ( r ) , (25)where Π ( r ) is the “free” Debye potential for the monopole gravitation given by (22) [9, 11]. In other words, in theeikonal approximation the Debye potential Π ( r ), becomes “distorted” in the presence of the potential V sr given inEq. (B23), by φ , a slowly varying function of r , such that (cid:12)(cid:12) ∇ φ (cid:12)(cid:12) ≪ k |∇ φ | . (26)When substituted into (24), the trial solution (25) yields n ∆Π ( r ) + k (cid:0) r g r (cid:1) Π ( r ) o φ ( r ) + Π ( r )∆ φ ( r ) ++ 2 (cid:0) ∇ Π ( r ) · ∇ φ ( r ) (cid:1) + 4 k V sr ( r )Π ( r ) φ ( r ) = O (cid:16) r g , J r Π (cid:17) . (27)As Π ( r ) is the solution of the homogeneous equation for the monopole gravitational field (22) [9, 11], the first termin (27) is zero. Then, we neglect the second term, Π ( r )∆ φ ( r ), because of (26). As a result, from the last two termswe have (cid:0) ∇ ln Π ( r ) · ∇ ln φ ( r ) (cid:1) = − k V sr ( r ) + O ( r g ) . (28)As we discussed above, we assume that contributions from deviations from the monopole are small and it issufficient to keep only terms to O (cid:0) r g , ( J /r )Π (cid:1) . Thus, to formally solve (28) we may present the solution for Π ( r )at a large distance from the monopole, which yields the well-known solution for the incident wave in the presence ofa gravitational monopole (see Eq. (23) in [9]):Π ( r ) = e ± ik (cid:0) z − r g ln k ( r − z ) (cid:1) + O ( r g ) . (29)To compute the gradient of Π ( r ), following [9], we represent the unperturbed trajectory of a ray of light as r ( t ) = r + k c ( t − t ) + O ( r g ) , (30)where k is the unit vector in the incident direction of the light ray’s propagation path and r represents the startingpoint. Following [9, 34, 35], we define b = [[ k × r ] × k ] to be the impact parameter of the unperturbed trajectoryof the light ray. The vector b is directed from the origin of the coordinate system toward the point of the closestapproach of the unperturbed path of light ray to that origin.With (30), we introduce the parameter τ = τ ( t ) along the path of the light ray (see details in Appendix B in [9]): τ = ( k · r ) = ( k · r ) + c ( t − t ) , (31)which may be positive or negative. Note that τ = z cos α where α is the angle between e z and k . Furthermore, τ = z when the z -axis of the chosen Cartesian coordinate system is oriented along the incident direction of the light ray.We can see that the quantity τ evolves from a negative value (representing a source at a large distance from the lens, α ≃ π ), through τ = 0 (the shortest distance from the lens where α = π/ α ≃ τ allows us to rewrite (30) as r ( τ ) = b + k τ + O ( r g ) , with || r ( τ ) || ≡ r ( τ ) = p b + τ + O ( r g ) . (32)Using (32), the gradient of Π (1) ( r ) from (29) may be computed as ∇ ln Π ( r ) = ± ik (cid:16) k (1 + r g r ) − r g b b (cid:0) τr (cid:1)(cid:17) + O ( r g ) . (33)As a result, (28) takes the form ± ik (cid:16)(cid:0) k (1 + r g r ) − r g b b (cid:0) τr (cid:1)(cid:1) · ∇ ln φ ( r ) (cid:17) = − V sr ( r ) + O ( r g ) . (34)As we want to identify the largest contribution from corrections to the monopole to light propagation, we keep onlylinear terms with respect to gravity. As a result, neglecting the r g -dependent terms in (34), we may present (28) as ± ik ( k · ∇ ) ln φ = − k V sr + O ( r g ) . (35)We may now compute the eikonal phase due to the short-range potential V sr . Using the representation of the lightray’s path as r = ( b , τ ) given by (32), we observe that (as was also shown in [9]) the gradient ∇ may be expressed interms of the variables along the path as ∇ = ∇ b + k d/dτ + O ( r g ), where ∇ b is the gradient along the direction ofthe impact parameter b and τ being the parameter taken along the path. Thus, the differential operator on the leftside of (35) is the derivative along the light ray’s path, namely ( k · ∇ ) = d/dτ .As a result, for (35) we have d ln φ ± dτ = ± ik V sr + O ( r g ) , (36)the solutions of which are φ ± ( b , τ ) = exp (cid:16) ± ik Z τ V sr ( b , τ ′ ) dτ ′ (cid:17) . (37)We therefore have the following two particular eikonal solutions of (24) for Π( r ):Π( r ) = Π ( r ) exp (cid:16) ± iξ b ( b , τ ) (cid:17) + O (cid:16) r g , J r Π (cid:17) , (38)where we introduced the eikonal phase ξ b ( b , τ ) = k Z τ V sr ( b , τ ′ ) dτ ′ + O ( r g ) . (39)The solution given by (38)–(39) was obtained under the eikonal condition (26) that allows us to consider the effectof the short-range potential due to gravitational multipoles (as shown in (23)–(24)) on the phase of the EM wave only,and not on its amplitude. This is similar to the thin lens approximation that is extensively used in the description ofmany problems on modern optics [22] and gravitational lensing [5]. That fact is captured by (39) where we assumethat light moves in a straight line before it reaches the lens and then it changes direction at τ = ( k · r ) = 0 and movesagain on a straight line towards the observer. Thus, the phase shift (39) occurs only on the second part of the path.Considering the structure of solutions (22) and (38), we note that the eikonal phase, ξ b ( b , s ) from (39), depends onthe vector of the impact parameter b and its orientation with respect to the solar rotational axis. Thus, the presenceof ξ b ( b , s ) in (38) is understood in the context of solution (22), where the sum over ℓ = kb also acts on the b -dependenteikonal phase. In general, this approach is similar to that of the Born approximation [25–27, 31] or path integrals inquantum mechanics [36–39]. This point will become more evident in Section V.In Appendix B 1, we compute the eikonal phases for two possible forms of the gravitational potential, valid in thegeneric case. In the case of the spatial trace-free (STF) multipole moments from (B5), the eikonal phase is given by(B27). However, in the case of the SGL, the gravitational potential of the Sun is that of an axisymmetric body bestcharacterized using zonal harmonics (B2) and may be expressed [40, 41] in terms of the usual dimensionless multipolemoments J ℓ : U = GMr n − ∞ X ℓ =2 J ℓ (cid:16) R ℓ r (cid:17) ℓ P ℓ (cid:16) s · x r (cid:17)o + O ( c − ) , (40)where s denotes the unit vector along the x -axis, P ℓ are the Legendre polynomials and the quantities M, J ..., J ℓ correspond to the multipole moments. Note that, in the case of an axisymmetric and rotating body with “north-southsymmetry” (i.e., a body that is symmetric under a reflection with respect to the plane of rotation), the expression(B2) contains only the ℓ = 2 , , , ... even moments.To determine the eikonal phase (39), we use a heliocentric coordinate system with its z -axis aligned with thewavevector k , so that k = (0 , , b = b n ξ ,coordinates on the image plane, x that is located at the distance z from the Sun, and the unit vector in the directionof the solar rotation axis, s : b = b (cos φ ξ , sin φ ξ , , (41) x = ρ (cos φ, sin φ, , (42) s = (sin β s cos φ s , sin β s sin φ s , cos β s ) . (43)In this coordinate system, the eikonal phase shift (39) accumulated by an EM wave propagating in the gravitationalfield of an axisymmetric body takes an elegant form (see discussion in Appendix B and the result (B21)): ξ b ( b , s ) = − kr g ∞ X n =2 J n n (cid:16) R ⊙ b (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] + O ( r g ) . (44)This result provides the context for our investigation below as it shows the explicit dependencies of the eikonal phaseon the orientation of the vector of the impact parameter with respect to the lens’ rotational axis. III. ELECTROMAGNETIC WAVE IN THE FIELD OF A STATIC EXTENDED GRAVITY LENS
Our next goal is to find a solution to the EM field in that region. We accomplish this objective using the approachdeveloped for classical diffraction theory, by finding the set of equations that determine the EM field via Debyepotentials and then matching these equations with the incident wave.
A. Solution for the radial function R ℓ ( r, b ) At this point, we already have all the key components needed to develop the solution for the Debye potentials in thecase of the long-range, spherically symmetric gravitational field produced by the solar monopole, and the static long-range gravitational field produced by deviations from the monopole, characterized using zonal harmonics. Following[9, 11], a particular solution for the Debye potential, Π, is obtained by combining results for Φ( φ ) from (19) and Θ( θ )from (20). The solution for the Debye potential takes the formΠ u = 1 r ∞ X ℓ =0 ℓ X m = − ℓ µ ℓ R ℓ ( r, b ) (cid:2) P ( m ) l (cos θ ) (cid:3)(cid:2) a m cos( mφ ) + b m sin( mφ ) (cid:3) + O (cid:16) r g , J r Π (cid:17) , (45)where the yet to be constructed R ℓ ( r, b ) is the radial function and µ ℓ , a m , b m are arbitrary and as yet unknownconstants. Note that the structure of the solution (45) preserves the angular symmetries of the monopole case givenby Φ( φ ) and Θ( θ ). The presence of the gravitational zonal harmonics is accounted for by the generalized radialfunction R ℓ ( r, b ) that now depends on b via the eikonal phase, as shown in (38).In the vacuum, the solutions for the electric and magnetic potentials of the incident wave, e Π and m Π , were foundto be given in terms of a single potential Π ( r, θ ) that is given by (22). In other words, the incident EM wave is notaffected by the gravitational field from the zonal harmonics of the extended Sun. Its form is identical to that of thefree EM wave propagating in monopole gravity, discussed in [9].Considering deviations from spherical symmery, we notice that, for large r , the potential V sr ( r ) in (24) can beneglected in comparison to the Coulomb potential U c ( r ) = 2 k r g /r and this equation reduces to the Coulomb equationdiscussed in [9] with the solution given by (22). The solution of (24) that is regular at the origin can thus bewritten asymptotically as a linear combination of the regular and irregular Coulomb wave functions F e ( kr g , kr )and G ℓ ( kr g , kr ), respectively [11, 26–28, 42], which are solutions of (24) in the absence of the potential V sr ( r ).Asymptotically, at large values of the argument ( kr ), these functions behave as [9, 11]: F ℓ ( kr g , kr ) ∼ sin (cid:16) k ( r + r g ln 2 kr ) + ℓ ( ℓ + 1)2 kr − πℓ σ ℓ (cid:17) , (46) G ℓ ( kr g , kr ) ∼ cos (cid:16) k ( r + r g ln 2 kr ) + ℓ ( ℓ + 1)2 kr − πℓ σ ℓ (cid:17) . (47)In the case of centrally symmetric potentials, since the Coulomb potential falls off slower than the centrifugalpotential (i.e., the ℓ ( ℓ + 1) /r term in (46) and (47)) at large distances, it dominates the asymptotic behavior of theeffective potential in every partial wave. Hence, we can generally look for a solution satisfying the following boundaryconditions [28]: R ℓ ( r ) ∼ r → nr ℓ +1 , (48) R ℓ ( r ) ∼ r →∞ F ℓ ( kr g , kr ) + tan δ ℓ G ℓ ( kr g , kr ) ∝ kr →∞ sin (cid:16) k ( r + r g ln 2 kr ) + ℓ ( ℓ + 1)2 kr − πℓ σ ℓ + δ ℓ (cid:17) , (49)where n is a normalization factor and F ℓ ( kr g , kr ) and G ℓ ( kr g , kr ) are solutions of (24) in the absence of the potential V sr ( r ), which, as we discussed above, are respectively regular and irregular at the origin. The real quantities δ ℓ ( k )introduced by these equations are the phase shifts due to the short-range potential V sr ( r ) (B23) in the presence ofthe Coulomb potential U c ( r ) = 2 k r g /r in (24). We note that δ ℓ ( k ) fully describes the non-Coulombic part of thescattering and vanishes when this short-range potential is absent.In the case of generic gravitational fields, we can satisfy the conditions (48)–(49) by choosing the function R ℓ ( r )as a linear combination of the two solutions (38), where δ ℓ ( k ) is replaced by the eikonal phase, ξ b ( b ). One way to dothat is by relying on the two solutions to (38), taken in the form of the incident and scattered waves [43], which arecorrespondingly given by the functions H − ℓ ( kr g , kr ) and H + ℓ ( kr g , kr ), and to show explicit dependence on the eikonalphase shift, ξ b ( b ), which can be captured in the following form: R ℓ ( r, b ) = 12 i (cid:16) H + ℓ ( kr g , kr ) e iξ b ( b ) − H − ℓ ( kr g , kr ) e − iξ b ( b ) (cid:17) , (50)where the Coulomb–Hankel functions H ( ± ) ℓ are related to the Coulomb functions by H ± ℓ ( kr g , kr ) = G ℓ ( kr g , kr ) ± iF ℓ ( kr g , kr ) (for discussion, see Appendix A of [9]) and their asymptotic behavior is given by (see Appendix F of [9]): H ± ℓ ( kr g , kr ) ∼ kr →∞ exp n ± i (cid:16) k ( r + r g ln 2 kr ) + ℓ ( ℓ + 1)2 kr − πℓ σ ℓ (cid:17)o , (51)where ξ b ( b ) in (50) is the eikonal phase shift that is accumulated by the EM wave along its entire path. The expressionfor this quantity is given by (39), which, for axisymmetric body, is computed by (44).The form of the radial function R ℓ from (50) captures our expectation that, in the presence of a potential V sr from(B23), the Coulomb–Hankel functions (which represent the radial free-particle wavefunction solutions of the homo-geneous equation (24)), become “distorted” by this short-range potential due to the gravitational mass multipoles.We can verify that R ℓ in the form of (50) also satisfies the asymptotic boundary conditions (48)–(49). Indeed, as thegravitational potential for the inner region of the Sun vanishes, the eikonal phase ξ b is zero for r < R ⊙ . Therefore, as r →
0, the radial function (50) becomes R ℓ ( r, b ) → F ℓ ( kr g , kr ), where the function F ℓ ( kr g , kr ) obeys the condition(48). Next, we consider another limit, when r → ∞ . Using the asymptotic behavior of H ± ℓ from (51), we see that, as r → ∞ , the radial function obeys the asymptotic condition (49) taking the form where the phase shift δ ℓ is given bythe eikonal phase ξ b introduced by (39).We may put the result (50) it in the following equivalent form: R ℓ ( r, b ) = cos ξ b ( b ) F ℓ ( kr g , kr ) + sin ξ b ( b ) G ℓ ( kr g , kr ) , (52)which explicitly shows the phase shift, ξ b ( b ), induced by the short-range extended gravity potential, clearly satisfyingthe boundary condition (49) with the quantity ξ ( b ) from (39) being the anticipated phase shift.To match the potentials (45) of the incident and scattered waves outside, the latter must be expressed in a similarform but with arbitrary coefficients. Only the function F ℓ ( kr g , kr ) may be used in the expression for the potentialinside the sphere since G ℓ ( kr g , kr ) becomes infinite at the origin. On the other hand, the scattered wave must vanish atinfinity. The Coulomb–Hankel functions H + ℓ ( kr g , kr ) are characterized by precisely this property, which makes themsuitable as representations of scattered waves. For large values of the argument kr , the result behaves as e ik ( r + r g ln 2 kr ) and the Debye potential Π ∝ e ik ( r + r g ln 2 kr ) /r for large r . Thus, at large distances from the sphere the scattered waveis spherical (with the ln term in the phase due to the modification by the Coulomb potential), with its center at theorigin r = 0. Accordingly, we use it in the expression for the scattered wave.Collecting results for the functions Φ( φ ) and Θ( θ ), respectively given by (19) and (20), and R ℓ ( r, b ) = H + ℓ ( kr g , kr ) e iξ b ( b ) from (38), to O ( r g , ( J /r )Π), we obtain the Debye potential for the scattered wave:Π s = ur ∞ X ℓ =0 ℓ X m = − ℓ a ℓ H + ℓ ( kr g , kr ) e iξ b ( b ) (cid:2) P ( m ) ℓ (cos θ ) (cid:3)(cid:2) a ′ m cos( mφ ) + b ′ m sin( mφ ) (cid:3) , (53)where a ℓ , a ′ m , b ′ m are arbitrary and as yet unknown constants.Representing the potential via F ℓ ( kr g , kr ) is appropriate. The trial solution to (24) for the electric and magneticDebye potentials relies on the radial function R ℓ ( r, b ) given by (52) and has the formΠ in = ur ∞ X ℓ =0 ℓ X m = − ℓ b ℓ n cos ξ b ( b ) F ℓ ( kr g , kr ) + sin ξ b ( b ) G ℓ ( kr g , kr ) o(cid:2) P ( m ) ℓ (cos θ ) (cid:3)(cid:2) a m cos( mφ ) + b m sin( mφ ) (cid:3) , (54)where b ℓ , a m , b m are arbitrary and yet unknown constants.The boundary (continuity) conditions (see discussion in [9, 22]), imposed on the quantities (A41) at some radius r = R ⋆ ⊙ = R ⊙ + r g , are written in full as ∂∂r h r e Π u + r e Π s √ ǫu i(cid:12)(cid:12)(cid:12) r = R ⋆ ⊙ = ∂∂r h r e Π in √ ǫu i(cid:12)(cid:12)(cid:12) r = R ⋆ ⊙ , (55) ∂∂r h r m Π u + r m Π s √ µu i(cid:12)(cid:12)(cid:12) r = R ⋆ ⊙ = ∂∂r h r m Π in √ µu i(cid:12)(cid:12)(cid:12) r = R ⋆ ⊙ , (56) h r e Π u + r e Π s √ ǫu i(cid:12)(cid:12)(cid:12) r = R ⋆ ⊙ = h r e Π in √ ǫu i(cid:12)(cid:12)(cid:12) r = R ⋆ ⊙ , (57) h r m Π u + m Π s √ µu i(cid:12)(cid:12)(cid:12) r = R ⋆ ⊙ = h r m Π in √ µu i(cid:12)(cid:12)(cid:12) r = R ⋆ ⊙ . (58)0We now make use of the symmetry of the geometry of the problem [22] and by applying the boundary conditions(55)–(58). We recall that we can use a single Debye potential Π in (53) and (54) to represent electric and magneticfields. We find that the constants a m and b m for the electric Debye potentials are a = 1, b = 0 and a m = b m = 0for m ≥
2. For the magnetic Debye potentials, we obtain a = 0, b = 1 and a m = b m = 0 for m ≥
2. The values areidentical for a ′ m and b ′ m .As a result, the solutions for the electric and magnetic potentials of the scattered wave, e Π s and m Π s , may be givenin terms of a single potential Π s ( r, θ ) (see [9] for details), which, to O ( r g , ( J /r )Π), is given by (cid:18) e Π s m Π s (cid:19) = (cid:18) cos φ sin φ (cid:19) Π s ( r, θ ) , where Π s ( r, θ ) = ur ∞ X ℓ =1 a ℓ H + ℓ ( kr g , kr ) e iξ b ( b ) P (1) ℓ (cos θ ) . (59)In a relevant scattering scenario, the EM wave and the Sun are well separated initially, so the Debye potentialfor the incident wave can be expected to have the same form as for the pure monopole case that includes only theCoulomb potential that is given by (22). Therefore, the Debye potential for the inner region has the form: (cid:18) e Π in m Π in (cid:19) = (cid:18) cos φ sin φ (cid:19) Π in ( r, θ ) , (60)with the potential Π in given, to O ( r g , ( J /r )Π), asΠ in ( r, θ ) = ur ∞ X ℓ =1 b ℓ n cos ξ b ( b ) F ℓ ( kr g , kr ) + sin ξ b ( b ) G ℓ ( kr g , kr ) o P (1) ℓ (cos θ ) . (61)We thus expressed all the potentials in the series (45) and any unknown constants can now be determined easily.If we now substitute the expressions (22), (59) and (60)–(61) into the boundary conditions (55)–(58), we obtain thefollowing linear relationships between the coefficients a ℓ and b ℓ : h E k i ℓ − ℓ + 1 ℓ ( ℓ + 1) e iσ ℓ F ′ ℓ ( kr g , kr ) + a ℓ (cid:16) H + ℓ ( kr g , kr ) e iξ b ( b ) (cid:17) ′ i(cid:12)(cid:12)(cid:12) r = R ⋆ ⊙ = b ℓ R ′ ℓ ( r, b ) (cid:12)(cid:12)(cid:12) r = R ⋆ ⊙ , (62) h E k i ℓ − ℓ + 1 ℓ ( ℓ + 1) e iσ ℓ F ℓ ( kr g , kr ) + a ℓ H + ℓ ( kr g , kr ) e iξ b ( b ) i(cid:12)(cid:12)(cid:12) r = R ⋆ ⊙ = b ℓ R ℓ ( r, b ) (cid:12)(cid:12)(cid:12) r = R ⋆ ⊙ , (63)where R ℓ ( r ) is from (52) and ′ = d/dr . We now define, for convenience, α ℓ and β ℓ as a ℓ = E k i ℓ − ℓ + 1 ℓ ( ℓ + 1) e iσ ℓ α ℓ and b ℓ = E k i ℓ − ℓ + 1 ℓ ( ℓ + 1) e iσ ℓ β ℓ . (64)From (62)–(63), we have: F ′ ℓ ( R ⋆ ⊙ ) + α ℓ H + ℓ ′ ( R ⋆ ⊙ ) e iξ b ( b ) = β ℓ R ′ ℓ ( R ⋆ ⊙ , b ) , (65) F ℓ ( R ⋆ ⊙ ) + α ℓ H + ℓ ( R ⋆ ⊙ ) e iξ b ( b ) = β ℓ R ℓ ( R ⋆ ⊙ , b ) , (66)where F ℓ ( R ⋆ ⊙ ) = F ℓ ( kr g , kR ⋆ ⊙ ) and H + ℓ ( R ⋆ ⊙ ) = H + ℓ ( kr g , kR ⋆ ⊙ ) with similar definitions for the derivatives of thesefunctions. Equations (65)–(66) may now be solved to determine the two sets of coefficients α ℓ and β ℓ : α ℓ = e − iξ b ( b ) F ℓ ( R ⋆ ⊙ ) R ′ ℓ ( R ⋆ ⊙ , b ) − F ′ ℓ ( R ⋆ ⊙ ) R ℓ ( R ⋆ ⊙ , b ) R ℓ ( R ⋆ ⊙ , b ) H + ℓ ′ ( R ⋆ ⊙ ) − R ′ ℓ ( R ⋆ ⊙ , b ) H + ℓ ( R ⋆ ⊙ ) , (67) β ℓ = F ℓ ( R ⋆ ⊙ ) H + ′ ℓ ( R ⋆ ⊙ ) − F ′ ℓ ( R ⋆ ⊙ ) H + ℓ ( R ⋆ ⊙ ) R ℓ ( R ⋆ ⊙ , b ) H + ℓ ′ ( R ⋆ ⊙ ) − R ′ ℓ ( R ⋆ ⊙ ) H + ℓ ( R ⋆ ⊙ , b ) . (68)Taking into account the asymptotic behavior of all the functions involved: namely (51) for H + ℓ and (46)–(47) for F ℓ and G ℓ , we have the following solution for the coefficients α ℓ and β ℓ : α ℓ = sin ξ b ( b ) , β ℓ = e iξ b ( b ) , (69)with ξ b ( b ) is the phase shift induced by the gravitational multipoles to the phase of the EM wave propagating throughthe solar system.1Therefore, using the value for a ℓ from (64), together with α ℓ from (69), we determine that the solution for thescattered potential (59) takes the formΠ s ( r, θ ) = E k ur ∞ X ℓ =1 i ℓ − ℓ + 1 ℓ ( ℓ + 1) e iσ ℓ sin ξ b ( b ) H + ℓ ( kr g , kr ) e iξ b ( b ) P (1) ℓ (cos θ ) , (70)which we can present asΠ s ( r, θ ) = E ik ur ∞ X ℓ =1 i ℓ − ℓ + 1 ℓ ( ℓ + 1) e iσ ℓ H + ℓ ( kr g , kr ) (cid:16) e iξ b ( b ) − (cid:17) P (1) ℓ (cos θ ) . (71)In the region outside the Sun, r > R ⋆ ⊙ , we may take the asymptotic form for the Coulomb–Hankel function andpresent (71) as Π s ( r, θ ) = − E k ur e ik ( r + r g ln 2 kr ) ∞ X ℓ =1 ℓ + 1 ℓ ( ℓ + 1) e i (2 σ ℓ + ℓ ( ℓ +1)2 kr ) (cid:16) e iξ b ( b ) − (cid:17) P (1) ℓ (cos θ ) . (72)As a result, using (22) and (71), we present the Debye potential in the region outside the Sun, r > R ⊙ , in thefollowing form:Π out ( r, θ ) = Π ( r, θ ) + Π s ( r, θ ) == E k ur ∞ X ℓ =1 i ℓ − ℓ + 1 ℓ ( ℓ + 1) e iσ ℓ n F ℓ ( kr g , kr ) + 12 i (cid:16) e iξ b ( b ) − (cid:17) H + ℓ ( kr g , kr ) o P (1) ℓ (cos θ ) . (73)Similarly, substituting the value for b ℓ from (64), together with β ℓ from (69), we determine the solution for theinner Debye potential (69) in the formΠ in ( r, θ ) = E k ur ∞ X ℓ =1 i ℓ − ℓ + 1 ℓ ( ℓ + 1) e i ( σ ℓ + ξ b ( b )) n cos ξ b ( b ) F ℓ ( kr g , kr ) + sin ξ b ( b ) G ℓ ( kr g , kr ) o P (1) ℓ (cos θ ) . (74)As solar gravity is rather weak, we may use the asymptotic expressions for F ℓ , G ℓ and H ± ℓ for r ≥ R ⊙ . Therefore,the radial function R ℓ ( r, b ) from (50) (or, equivalently, from (52)) may be given as R ℓ ( r, b ) = 12 i (cid:16) H + ℓ ( kr g , kr ) e iξ b ( b ) − H − ℓ ( kr g , kr ) e − iξ b ( b ) (cid:17) = e − iξ b ( b ) n F ℓ ( kr g , kr ) + 12 i (cid:0) e iξ b ( b ) − (cid:1) H + ℓ ( kr g , kr ) o , (75)where ξ b = ξ b ( b ) is the eikonal phase.As a result, outside the Sun, we may present (74) in the following equivalent form:Π in ( r, θ ) = E k ur ∞ X ℓ =1 i ℓ − ℓ + 1 ℓ ( ℓ + 1) e iσ ℓ n F ℓ ( kr g , kr ) + 12 i (cid:16) e iξ b ( b ) − (cid:17) H + ℓ ( kr g , kr ) o P (1) ℓ (cos θ ) . (76)The solution for the Debye potential, Π( r, θ ) from (76), describing the propagation of the EM wave on the back-ground of the static gravitational monopole and the short-range multipole gravitational field takes the formΠ in ( r, θ ) = E k ur ∞ X ℓ =1 i ℓ − ℓ + 1 ℓ ( ℓ + 1) e iσ ℓ F ℓ ( kr g , kr ) P (1) ℓ (cos θ ) ++ E ik ur ∞ X ℓ =1 i ℓ − ℓ + 1 ℓ ( ℓ + 1) e iσ ℓ (cid:16) e iξ b ( b ) − (cid:17) H (+) ℓ ( kr g , kr ) P (1) ℓ (cos θ ) + O (cid:16) r g , J r Π (cid:17) . (77)The first term in (76) is the Debye potential of an EM wave propagating in a vacuum but modified by the gravityof extended Sun. The second term represents the effect of the solar gravitational multipoles on the propagation ofthe EM waves. Notice that, as the distance increases, this term approaches the form of the Debye potential Π s forthe scattered EM field given by (72).Thus, we have identified all the Debye potentials involved in the Mie problem [21], namely the potential Π givenby (22) representing the incident EM field, the potential Π s from (72) describing the scattered EM field, and thepotential Π in from (76) total field.2 IV. GENERAL SOLUTION FOR THE EM FIELD
To describe the scattering of light by the extended Sun, we use solutions for the Debye potential representing thescattered EM wave (72), and the EM wave (77). The presence of the Sun itself is not yet captured. For this, we needto set additional boundary conditions that describe the interaction of the Sun with the incident radiation. Similarlyto [9, 12], we apply the fully absorbing boundary conditions that represent the physical size and the surface propertiesof the Sun [11, 44].We begin with the area that lies outside the Sun where three regions are present, namely (i) the shadow region, (ii)the geometric optics region, and (iii) the interference region. Clearly, as far as imaging with the SGL is concerned, theinterference region is of the greatest importance. This is where the SGL focuses light coming from a distant object,forming an image.
A. Fully absorbing boundary conditions
Boundary conditions representing the opaque Sun were introduced in [45] and were used in [9, 12]. Here we usethese conditions again. Specifically, to set the boundary conditions, we rely on the semiclassical analogy between thepartial momentum, ℓ , and the impact parameter, b , that is given as ℓ = kb [24, 25].To set the boundary conditions, we require that rays with impact parameters b ≤ R ⋆ ⊙ = R ⊙ + r g are completelyabsorbed by the Sun [9]. Thus, the fully absorbing boundary condition signifies that all the radiation intercepted bythe body of the Sun is fully absorbed by it and no reflection or coherent reemission occurs. All intercepted radiation istransformed into some other forms of energy, notably heat. Thus, we require that no scattered waves exist with impactparameter b ≪ R ⋆ ⊙ or, equivalently, for ℓ ≤ kR ⋆ ⊙ . Such formulation relies on the concept of the semiclassical impactparameter b and its relationship with the partial momentum, ℓ , as ℓ = kb . (A relevant discussion on this relationbetween ℓ and b is on p. 29 of [46] with reference to [47].) In terms of the boundary conditions, this means that weneed to subtract the scattered waves from the incident wave for ℓ ≤ kR ⋆ ⊙ , as was discussed in [9]. Furthermore, as itwas shown in [44], the fully absorbing boundary conditions introduce a fictitious EM field that precisely compensatesthe incident field in the area behind the Sun. This area has the shape of a rotational hyperboloid that starts directlyat the solar surface behind the Sun and extends to the vertex of the hyperboloid at z = R ⊙ / r g ≃ B. The Debye potential for the region outside the Sun
To implement the boundary conditions for the EM wave outside the Sun, we realize that the total EM field in thisregion is given as the sum of the incident and scattered waves, Π = Π + Π s , with these two potentials given by (22)and (72), correspondingly. Accordingly, we use (73), which represents the Debye potential in the region of interestandis given asΠ( r, θ ) = Π ( r, θ ) + Π s ( r, θ ) = E k ur ∞ X ℓ =1 i ℓ − ℓ + 1 ℓ ( ℓ + 1) e iσ ℓ n F ℓ ( kr g , kr ) + 12 i (cid:0) e iξ b ( b ) − (cid:1) H + ℓ ( kr g , kr ) o P (1) ℓ (cos θ ) . (78)Next, relying on the representation of the regular Coulomb function F ℓ via incoming, H + ℓ , and outgoing, H − ℓ , wavesas F ℓ = ( H + ℓ − H − ℓ ) / i (discussed in [9] and also by the expression given after (50)), we may express the Debyepotential (78) asΠ( r, θ ) = E ik ur ∞ X ℓ =1 i ℓ − ℓ + 1 ℓ ( ℓ + 1) e iσ ℓ n e iξ b ( b ) H + ℓ ( kr g , kr ) − H − ℓ ( kr g , kr ) o P (1) ℓ (cos θ ) . (79)This form of the combined Debye potential is convenient for implementing the fully absorbing boundary conditionsdiscussed in Sec. IV A. Specifically, subtracting from (79) the outgoing wave (i.e., ∝ H (+) ℓ ) for the impact parameters b ≤ R ⋆ ⊙ or equivalently for ℓ ∈ [1 , kR ⋆ ⊙ ], we haveΠ( r, θ ) = E ik ur ∞ X ℓ =1 i ℓ − ℓ + 1 ℓ ( ℓ + 1) e iσ ℓ n e iξ b ( b ) H + ℓ ( kr g , kr ) − H − ℓ ( kr g , kr ) o P (1) ℓ (cos θ ) −− E ik ur kR ⋆ ⊙ X ℓ =1 i ℓ − ℓ + 1 ℓ ( ℓ + 1) e iσ ℓ e iξ b ( b ) H + ℓ ( kr g , kr ) P (1) ℓ (cos θ ) , (80)3or, equivalently, coming back to the form (78),Π( r, θ ) = Π ( r, θ ) + E ik ur ∞ X ℓ =1 i ℓ − ℓ + 1 ℓ ( ℓ + 1) e iσ ℓ (cid:0) e iξ b ( b ) − (cid:1) H + ℓ ( kr g , kr ) P (1) ℓ (cos θ ) −− E ik ur kR ⋆ ⊙ X ℓ =1 i ℓ − ℓ + 1 ℓ ( ℓ + 1) e iσ ℓ e iξ b ( b ) H + ℓ ( kr g , kr ) P (1) ℓ (cos θ ) . (81)This is a rather complex expression. It requires the tools of numerical analysis to fully explore its behavior andthe resulting EM field [46–48]. However, in most practically important applications, we need to know the field inthe forward direction. Furthermore, our main interest is to study the largest impact of the extended gravity on lightpropagation, which corresponds to the smallest values of the impact parameter. In this situation, we may simplifythe result (81) by taking into account the asymptotic behavior of the function H + ℓ ( kr g , kr ), considering the field atlarge heliocentric distances, such that kr ≫ ℓ , where ℓ is the order of the Coulomb function (see p. 631 of [49]). For kr → ∞ and also for r ≫ r t = p ℓ ( ℓ + 1) /k (see [9, 12]), such an expression is given in the form [11]):lim kr →∞ H ± ℓ ( kr g , kr ) ∼ exp h ± ik (cid:0) r + r g ln 2 kr (cid:1) + ℓ ( ℓ + 1)2 kr + σ ℓ − πℓ (cid:17)i + O (cid:0) ( kr ) − , r g (cid:1) , (82)which includes the contribution from the centrifugal potential in the radial equation (21) (see e.g., Appendix C of[11], Appendix A in [50] or [48]). In fact, expression (82) extends the argument of (51) to shorter distances, closerto the turning point of the potential (see the relevant discussion in Appendix F of [9]). By including the extendedcentrifugal term in (82) (i.e., shown by the terms with various powers of ℓ ( ℓ + 1) / kr ), we can now better describethe bending of the trajectory of a light ray under the combined influence of extended gravity.We use the approximate behavior of H + ℓ given by (82) and use it in (81) to present the solution for the Debyepotential in the following form:Π( r, θ ) = Π ( r, θ ) + ue ik ( r + r g ln 2 kr ) r n E k kR ⋆ ⊙ X ℓ =1 ℓ + 1 ℓ ( ℓ + 1) e i (cid:0) σ ℓ + ℓ ( ℓ +1)2 kr (cid:1) P (1) ℓ (cos θ ) −− E k ∞ X ℓ = kR ⋆ ⊙ ℓ + 1 ℓ ( ℓ + 1) e i (cid:0) σ ℓ + ℓ ( ℓ +1)2 kr (cid:1)(cid:0) e i ξ b ( b ) − (cid:1) P (1) ℓ (cos θ ) o + O (cid:16) r g , J r Π (cid:17) == Π ( r, θ ) + Π bc ( r, θ ) + Π G ( r, θ ) . (83)The first term in (83), Π ( r, θ ), is the Debye potential that represents the incident EM wave propagating in thevacuum on the background of a post-Newtonian gravity field produced by a gravitational monopole. The solution forΠ ( r, θ ) is known and is given by (22) in the form of infinite series with respect to partial momenta, ℓ (see [9, 11]).The second term in (83), Π bc ( r, θ ), is due to the physical obscuration introduced by the Sun and was derived byapplying the fully absorbing boundary conditions. This term is responsible for the geometric shadow behind the Sun.The third term in (83), Π G ( r, θ ), quantifies the contribution of the extended gravitational field to the scattering ofthe EM wave.With the solution for the Debye potential given by (83), and with the help of (5)–(10) (also see [9]), we may nowcompute the EM field in the various regions involved. Given the smallness of the ratio r g J R ⊙ /r , we may neglectthe distance-dependent effects of the solar extended gravity on the amplitude of the EM wave. Thus, the extendedgravity contributes to the delay of the EM wave and is fully accounted for by the solution for the Debye potentials.Therefore, we can use the following expressions to construct the EM field in the static, gravity field produced by anextended gravity (see details in [9, 11]): ˆ D r ˆ B r ! = (cid:18) cos φ sin φ (cid:19) e − iωt α ( r, θ, φ ) , ˆ D θ ˆ B θ ! = (cid:18) cos φ sin φ (cid:19) e − iωt β ( r, θ, φ ) , ˆ D φ ˆ B φ ! = (cid:18) − sin φ cos φ (cid:19) e − iωt γ ( r, θ, φ ) , (84)with the quantities α, β and γ computed from the known Debye potential, Π, as α ( r, θ, φ ) = 1 u n ∂ ∂r h r Π u i + k u h r Π u io + O (cid:16)(cid:0) u (cid:1) ′′ (cid:17) , (85) β ( r, θ, φ ) = 1 u r ∂ (cid:0) r Π (cid:1) ∂r∂θ + ik (cid:0) r Π (cid:1) r sin θ , (86)4 γ ( r, θ, φ ) = 1 u r sin θ ∂ (cid:0) r Π (cid:1) ∂r + ikr ∂ (cid:0) r Π (cid:1) ∂θ . (87)This completes the solution for the Debye potentials on the background of a spherically symmetric, static gravita-tional field of the Sun. We will use (84)–(87) to compute the relevant EM fields. C. EM field in the shadow region
In the shadow behind the Sun (i.e., for impact parameters b ≤ R ⋆ ⊙ ) the EM field is represented by the Debyepotential of the shadow, Π sh , which is given asΠ sh ( r, θ ) = Π ( r, θ ) + ue ik ( r + r g ln 2 kr ) r E k kR ⋆ ⊙ X ℓ =1 ℓ + 1 ℓ ( ℓ + 1) e i (cid:0) σ ℓ + ℓ ( ℓ +1)2 kr (cid:1) P (1) ℓ (cos θ ) + O (cid:16) r g , J r Π (cid:17) , (88)where Π ( r, θ ) is well represented by (22). As discussed in [9, 44], the potential (88) produces no EM field. In otherwords, there is no light in the shadow. Furthermore, as the solar boundary is rather diffuse, there is no expectationfor a Poisson–Arago bright spot to form in this region. D. EM field outside the shadow
In the region behind the Sun but outside the solar shadow (i.e., for light rays with impact parameters b > R ⊙ )which includes both the geometric optics and interference regions (in the immediate vicinity of the focal line), theEM field is derived from the Debye potential given by the remaining terms in (83) to O ( r g , J /r ) asΠ( r, θ ) = Π ( r, θ ) − ue ik ( r + r g ln 2 kr ) r E k ∞ X ℓ = kR ⋆ ⊙ ℓ + 1 ℓ ( ℓ + 1) e i (cid:0) σ ℓ + ℓ ( ℓ +1)2 kr (cid:1)(cid:16) e i ξ b ( b ) − (cid:17) P (1) ℓ (cos θ ) , (89)where for the geometric optics region the potential Π ( r, θ ) is well represented by (22).Expression (89) is our main result for the regions outside the shadow region, i.e., b ≥ R ⋆ ⊙ . It contains all theinformation needed to describe the total EM field originating from an incident Coulomb-modified plane wave thatpassed through the region of the extended solar gravity field, characterized by the distance dependence that diminishesas r − or faster.To evaluate the total solution for the Debye potential (89), we present it in the following compact form:Π( r, θ ) = Π ( r, θ ) + E f G ( r, θ, φ ) ue ik ( r + r g ln 2 kr ) r , (90)where the extended gravity scattering amplitude f G ( r, θ, φ ) is given by f G ( r, θ, φ ) = − k ∞ X ℓ = kR ⋆ ⊙ ℓ + 1 ℓ ( ℓ + 1) e i (cid:0) σ ℓ + ℓ ( ℓ +1)2 kr (cid:1)(cid:16) e i ξ b ( b ) − (cid:17) P (1) ℓ (cos θ ) + O (cid:16) r g , J r Π (cid:17) . (91)We note that because of the contribution from the centrifugal potential in (82), the scattering amplitude f p ( r, θ ) isnow also a function of the heliocentric distance [12]. This is not the case in typical problems describing nuclear andatomic scattering [24, 25, 51, 52]. However, as we observed in [9, 12, 50], when we are interested in the trajectories oflight rays, the presence of such dependence and especially the ∝ /r term in the phase of the scattering amplitude (91)allows us to properly describe the bending of the light rays in the presence of gravity together with the contributionfrom deviations from spherical symmetry.As a result, the Debye potential takes the formΠ G ( r, θ ) = E f G ( r, θ, φ ) ue ik ( r + r g ln 2 kr ) r , (92)with the extended gravity scattering amplitude f G ( r, θ ) given by (91). We use these expressions to derive the compo-nents of the EM field produced by this wave. For this, we substitute (92)–(91) in the expressions (85)–(87) to derive5the factors α ( r, θ, φ ) , β ( r, θ, φ ) and γ ( r, θ, φ ), which to O ( r g , J /r ) are computed to be: α ( r, θ, φ ) = − E ue ik ( r + r g ln 2 kr ) k r ∞ X ℓ = kR ⋆ ⊙ ( ℓ + ) e i (cid:0) σ ℓ + ℓ ( ℓ +1)2 kr (cid:1)(cid:16) e i ξ b ( b ) − (cid:17) P (1) ℓ (cos θ ) n − ℓ ( ℓ + 1)4 k r o , (93) β ( r, θ, φ ) = E ue ik ( r + r g ln 2 kr ) ikr ×× ∞ X ℓ = kR ⋆ ⊙ ℓ + ℓ ( ℓ + 1) e i (cid:0) σ ℓ + ℓ ( ℓ +1)2 kr (cid:1)(cid:16) e i ξ b ( b ) − (cid:17)n ∂P (1) ℓ (cos θ ) ∂θ (cid:16) − ℓ ( ℓ + 1)2 u k r (cid:17) + P (1) ℓ (cos θ )sin θ o , (94) γ ( r, θ, φ ) = E ue ik ( r + r g ln 2 kr ) ikr ×× ∞ X ℓ = kR ⋆ ⊙ ℓ + ℓ ( ℓ + 1) e i (cid:0) σ ℓ + ℓ ( ℓ +1)2 kr (cid:1)(cid:16) e i ξ b ( b ) − (cid:17)n ∂P (1) ℓ (cos θ ) ∂θ + P (1) ℓ (cos θ )sin θ (cid:16) − ℓ ( ℓ + 1)2 u k r (cid:17)o , (95)where we neglected small terms that behave as ∝ i/ ( u kr ); terms ∝ ikr g /ℓ were also omitted because of the largepartial momenta involved, ℓ ≥ kR ⋆ ⊙ . Terms in both of these groups are negligably small when compared to the leadingterms in each of these expressions above (a similar conclusion was reached in [11, 12].)This is an important result as it allows us to describe the EM field in all the regions of interest for the SGL, namelythe strong and weak interference regions and the region of geometric optics. V. EM FIELD IN THE INTERFERENCE REGION
Results from the previous section allow us to study optical properties of the SGL in the case of extended gravitationallens. Our primary concern is the strong interference region: the area behind the Sun, reachable by light rays withimpact parameters b > R ⋆ ⊙ . The focal region of the SGL begins where r > b / r g and 0 ≤ θ ≃ p r g /r . The EMfield is derived from the Debye potential (90)–(91), given by the factors α , β and γ from (93)–(95). In the stronginterference region, these expressions take the following form [9, 11, 13]: α ( r, θ, φ ) = − E ue ik ( r + r g ln 2 kr ) k r ∞ X ℓ = kR ⋆ ⊙ ( ℓ + ) e i (cid:0) σ ℓ + ℓ ( ℓ +1)2 kr (cid:1)(cid:16) e i ξ b ( b ) − (cid:17) P (1) ℓ (cos θ ) (cid:16) O (cid:0) r g r , r g , J (cid:1)(cid:17) , (96) γ ( r, θ, φ ) = β ( r, θ, φ ) = E ue ik ( r + r g ln 2 kr ) ikr ×× ∞ X ℓ = kR ⋆ ⊙ ℓ + ℓ ( ℓ + 1) e i (cid:0) σ ℓ + ℓ ( ℓ +1)2 kr (cid:1)(cid:16) e i ξ b ( b ) − (cid:17)n ∂P (1) ℓ (cos θ ) ∂θ + P (1) ℓ (cos θ )sin θ o(cid:16) O (cid:0) r g r , r g , J (cid:1)(cid:17) . (97)We recognize that (96)–(97) represent the scattered EM field in the interference region. As evident in the structureof the expression (90), the term ( e i ξ b ( b ) −
1) present in these expressions leads to formation of two waves – that onethat is ∝ e i ξ b ( b ) is the EM wave due to diffraction of light by the gravitational multipoles, while the ∝ − ∝ − e i ξ b ( b ) − ℓ ≫ O ( θ )) to be evaluated with the method of stationary phase (with β ( r, θ, φ ) = γ ( r, θ, φ )): α ( r, θ, φ ) = − E ue ik ( r + r g ln 2 kr ) k r Z ∞ ℓ = kR ⋆ ⊙ ℓdℓe i (cid:0) σ ℓ + ℓ kr +2 ξ b ( b ) (cid:1) P (1) ℓ (cos θ ) (cid:16) O (cid:0) θ, r g r , r g , J (cid:1)(cid:17) , (98) γ ( r, θ, φ ) = E ue ik ( r + r g ln 2 kr ) ikr Z ∞ ℓ = kR ⋆ ⊙ dℓℓ e i (cid:0) σ ℓ + ℓ kr +2 ξ b ( b ) (cid:1)n ∂P (1) ℓ (cos θ ) ∂θ + P (1) ℓ (cos θ )sin θ o(cid:16) O (cid:0) θ, r g r , r g , J (cid:1)(cid:17) . (99)To evaluate these expressions in the interference region and for 0 ≤ θ ≃ p r g /r , we use the asymptotic represen-tation for P ℓ (cos θ ) and ℓ ≫ P ℓ (cos θ ) = r θ sin θ J (cid:0) ℓθ (cid:1) + O ( θ ) . (100)6For improved explicit two-term uniformly valid asymptotic form of this expression, check [56]. We use the expression P (1) ℓ (cos θ ) = − ∂P ℓ (cos θ ) ∂θ = ℓJ ( ℓθ ) + θJ ( ℓθ ) + O ( θ ) , (101)to derive P (1) ℓ (cos θ )sin θ = ℓ (cid:16) J ( ℓθ ) + J ( ℓθ ) (cid:17) , dP (1) ℓ (cos θ ) dθ = ℓ (cid:16) J ( ℓθ ) − J ( ℓθ ) (cid:17) . (102)Using expressions (101) and (102) in (98)–(99), we have α ( r, θ, φ ) = − E ue ik ( r + r g ln 2 kr ) k r Z ∞ ℓ = kR ⋆ ⊙ ℓ dℓ J ( ℓθ ) e i (cid:0) σ ℓ + ℓ kr +2 ξ b ( b t ) (cid:1)(cid:16) O (cid:0) θ, r g r , r g , J (cid:1)(cid:17) , (103) γ ( r, θ, φ ) = E ue ik ( r + r g ln 2 kr ) ikr Z ∞ ℓ = kR ⋆ ⊙ ℓdℓ J ( ℓθ ) e i (cid:0) σ ℓ + ℓ kr +2 ξ b ( b ) (cid:1)(cid:16) O (cid:0) θ, r g r , r g , J (cid:1)(cid:17) . (104)In the case of a pure gravitational monopole, the eikonal phase ξ b ( b ) in (103)–(104) is absent. In that case, theseintegrals can be evaluated using the method of stationary phase, leading to the well-known result [9, 11, 57] with thePSF ∝ J of a monopole (i.e., a point mass or spherically-symmetric) lens.However, in the presence of ξ b ( b ), the method of stationary phase is not applicable as the expressions (103)–(104)now have angular dependence that is not captured by the integrals. Therefore, we need a method that can addressthis to by transforming these integrals into an appropriate form that captures such a dependence.To evaluate these integrals, we developed what we call the angular eikonal method . This approach entails replacingthe Bessel functions J and J with their integral representations : J ( ℓθ ) = 12 π Z π dφ ξ e − iℓθ cos( φ ξ − φ ) , J ( ℓθ ) = i π Z π dφ ξ cos( φ ξ − φ ) e − iℓθ cos( φ ξ − φ ) . (105)These expressions recognize the fact that, to describe the geometry of the problem, we selected a heliocentric coordinatesystem whose z axis is co-linear with the wavevector, k , of the incident EM wave. The expressions in (105) are integralsover the azimuthal angle φ ξ , representing the orientation of the unit vector of the impact parameter b , as given by(43). The expressions (105) allow us to rewrite (103)–(104), to O (cid:0) θ, r g /r, r g (cid:1) , in the following 2-dimensional form: α ( r, θ, φ ) = E ue ik ( r + r g ln 2 kr ) ik r π Z π dφ ξ cos( φ ξ − φ ) Z ∞ ℓ = kR ⋆ ⊙ ℓ dℓe i (cid:0) σ ℓ + ℓ kr +2 ξ b ( b ) − ℓθ cos( φ ξ − φ ) (cid:1) , (106) γ ( r, θ, φ ) = E ue ik ( r + r g ln 2 kr ) ikr π Z π dφ ξ Z ∞ ℓ = kR ⋆ ⊙ ℓdℓe i (cid:0) σ ℓ + ℓ kr +2 ξ b ( b ) − ℓθ cos( φ ξ − φ ) (cid:1) . (107)The step presented above correctly captures the functional form of the integrand, which is azimuthally perturbedby the eikonal phase shift, 2 ξ b ( b ), whose presence breaks the spherical symmetry present in the case of a gravitationalmonopole. Technically, this step could have been done much earlier, in the Debye potential of the incident wave(22) that still possesses the symmetries representative of Coulomb-scattering. However, doing it that early wouldobscure the presentation of the overall solution. As it is known, solving the time-independent Schr¨odinger equationin the presence of the Coulomb potential (representing a point source) is a well-posed problem. As demonstrated by(15)–(22), this problem reduces to solving the relevant wave equation by implementing separation of variables thatresults in a well-known solution [9]. In the case when the scattering potential is not spherically symmetric, separationof variables, in general, is not possible. Thus, other methods are needed.For gravitational lensing in a weak gravitational field, characterized by a scattering potential with only smalldeviations from spherical symmetry, we may use the eikonal approximation to identify the eikonal phase shift thatcorresponds to that particular scattering potential (see details in Section II C). This eikonal phase shift is effectively arepresentation of the well-known thin lens approximation [5]. However, our variables in (103)–(104) were still given interms of the monopole case. This is where we recognized that the integral expressions (105) may be used to solve the Note that we can use the same representations of these functions with positive sign in the phase. The result is identical as it onlyreplaces the integrand with its complex conjugate, but it leaves the real-valued result unaffected. dφ ξ in (106)–(107) act not only the monopole part of the phase, − ℓθ cos( φ ξ − φ ) as in (103)–(104),but on the entire phase, which now includes the eikonal phase shift term 2 ξ b ( b ) due to the gravitational multipoles.This outlines the logic behind the angular eikonal method .Lastly, we mentioned that the resulting quantities (106)–(107) determine the EM field in the strong interferenceregion of the SGL. Below, we find that these expressions can be evaluated using the method of stationary phase.Furthermore, as we know [9, 11], in the interference region the factor α determining the radial components of the EMfield is very small, behaving as α ( r, θ, φ ) ≃ p r g /r γ ( r, θ, φ ). Thus, this factor will yield a negligible contribution tothe Poynting vector and it may be omitted. Therefore, in the discussion below we focus on the factor β ( r, θ, φ ) only. A. Integral over the vector impact parameter
As we mentioned above, the integral over dφ ξ in (107) properly acts not only on the monopole term of the phase, − ℓθ cos( φ ξ − φ ), but on the entire phase 2 σ ℓ + ℓ k ˜ r + 2 ξ b − ℓθ cos( φ ξ − φ ), that now includes the contributions fromthe parts that perturb the spherically symmetric gravitational potential via the eikonal phase, ξ b . In the case of anaxisymmetric gravitational field, this perturbation is given by (B21): ξ b ( b ) = − kr g ∞ X n =2 J n n (cid:16) R ⊙ b (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] . (108)For convenience, we introduce ξ b ( b ) = − kr g ψ ( b ) , ψ ( b ) = ∞ X n =2 J n n (cid:16) R ⊙ b (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] . (109)Furthermore, for ℓ ≫ kr g , we evaluate σ ℓ as [44]: σ ℓ = − kr g ln ℓ. (110)This form agrees with the other known forms of σ ℓ [58, 59] that are approximated for large ℓ .We rely on the semiclassical approximation that connects the partial momentum, ℓ , to the impact parameter, b forsmall angles θ (or large distances from the Sun, R ⊙ /r < b/r ≪ ℓ ≃ kb, (111)we obtain ϕ ( b ) = 2 σ ℓ + ℓ kr + 2 ξ b ( b ) − ℓθ cos( φ ξ − φ ) (cid:1) == k n b r − bθ cos( φ ξ − φ ) − r g (cid:16) ln kb + ∞ X n =2 J n n (cid:16) R ⊙ b (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] (cid:17)o . (112)or, compactly, using (109): ϕ ( b ) = k n b r − bθ cos( φ ξ − φ ) − r g (cid:0) ln kb + ψ ( b ) (cid:1)o . (113)We recognize that the vector of the impact parameter has the form b = b (cos φ ξ , sin φ ξ , θ to a point on the image plane with coordinates ( r, θ, φ ) that has the form θ = θ (cos φ, sin φ, bθ cos( φ ξ − φ ) = ( b · θ ), therefore ϕ ( b ) = k n r (cid:0) b − brθ cos( φ ξ − φ ) (cid:1) − r g (cid:0) ln kb + ψ ( b ) (cid:1)o = ≃ k n r (cid:0) b − r θ (cid:1) − r g (cid:0) ln kb + ψ ( b ) (cid:1)o . (114)Thus, the phase ϕ ( b ) represents the Fermat potential that governs the gravitational lensing phenomena [3–5, 60].8 b e b ke φ ξ FIG. 2: The basis vectors used to characterize the vector deflection angle. e b and k are in the plane spanned by the incidentlight ray and the center of the lens; e φ ξ is normal to this plane. The multipole moments change the amount by which theincident ray is deflected compared to the effect of the monopole (dashed arrow), but also lift the ray out of the plane spannedby the incident ray and the center of the lens. As a result, we can present (107) as γ ( r, θ, φ ) = E ue ik ( z + r g ln 2 kr ) kir π Z d b exp h ik (cid:16) r ( b − r θ ) − r g (cid:0) ln kb + ψ ( b ) (cid:1)(cid:17)i . (115)The integral in (115) is known rather well. It was obtained using different methods and tools by several authors.For instance, a similar integral formula for the lensed wave amplitude was obtained using the scalar theory of light in[57, 61–63]; by using the Fresnel–Kirchhoff diffraction formula [22]; and it was also obtained using the path integralformalism [36, 37] in [38, 39]. However, all previous efforts discussed primarily a monopole case. Our expression (115)generalizes these previously obtained results via the presence of the eikonal phase shift term, − r g ψ ( b ), to the caseof a lens with arbitrary multipole structure, which is explicitly captured by (109).We also note that all the previous results were obtained using the scalar theory, considering only the amplitudeof the EM wave. A unique feature of our approach is that we are able to reconstruct the entire vector structureof the EM field (e.g., using (84) together with (85)–(87)). This is an important capability when we consider thethree-dimensional behavior of a vector theory, for instance, polarization of the EM wave as it propagates through arefractive medium without spherical symmetry. Thus, our approach supersedes and generalizes all previous resultsobtained for gravitation lensing in post-Newtonian gravity.To put the entire problem in the proper context related to the geometry of light propagation in the vicinity of agravitating body, we consider the total gravitational delay, d ( b ), acquired by the EM wave as it travels in the gravityfield of the extended lens. This delay contributes to the phase shift ξ b ( b ) = kd ( b ) and using (112) is identified as d ( b ) = − r g (cid:16) ln kb + ∞ X n =2 J n n (cid:16) R ⊙ b (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] (cid:17) . (116)This is the generalization of the classic Shapiro time delay to the case of an extended axisymmetric gravitational lens.This delay corresponds to the total gravitational deflection angle acquired by a light ray or, equivalently, rotation ofthe wavefront of the EM wave. Using the expression (32) for the radius vector of the EM wave, together with b givenby (41), we compute this angle as θ g = − ∇ d ( b ) = − n e b ∂d ( b ) ∂b + e φ ξ ∂d ( b ) b ∂φ ξ + k ∂d ( b ) ∂τ o , (117)where the basis vector e b is the unit vector in the direction of the vector of the impact parameter b and e φ ξ is theunit vector in the azimuthal direction and is orthogonal to b and k (Fig. 2).Note that expression (108) for the eikonal phase obtained in Appendix B 2 a was derived under the thin lensapproximation where the distances travelled by the EM wave from the source to the lens, τ = ( k · r ), and from thelens to the observer, τ = ( k · r ), are much larger than the impact parameter, namely b/ | τ | ≪ b/ | τ | ≪
1. Thisis the reason why (116) does not depend on τ , thus yielding a vanishing derivative with respect to τ in (117).With these considerations in mind, we compute the vector of the total angle of the gravitational deflection of lightas the light ray passes in the vicinity of an extended axisymmetric lens: θ g = 2 r g b n e b − ∞ X n =2 J n (cid:16) R ⊙ b (cid:17) n sin n β s (cid:16) e b cos[ n ( φ ξ − φ s )] + e φ ξ sin[ n ( φ ξ − φ s )] (cid:17)o . (118)9The first term in (118) is the Einstein deflection angle in the gravity field of a spherically symmetric matterdistribution (i.e., in the presence of a monopole or point mass). The second term with J n describes the effect of themultipole moments as a sum of a) an additional deflection toward or away from the optical axis (the line parallel tothe incoming ray of light that intersects the lens at the center), and b) a deflection away from the plane defined bythe incoming ray and the center of the lens.To further appreciate the geometry that is represented by this e φ ξ term, consider the J case, when n = 2. Thereare four principal directions (which depend on the orientation of the impact parameter given by the angle φ ξ ) forwhich this second term is zero. These directions correspond to the cusps well-known astroid caustic of the quadrupolelens (see Sec. V E below). For all other angles, light is deflected away from the optical axis, lifted out of the planespanned by the direction of the incident light ray and the center of the lens (Fig. 2). These rays of light never intersectthe optical axis; therefore, an observer at the optical axis sees only light from the principal directions. Thus we caninstantly see how Eq. (118) gives rise to the famous Einstein-cross that appears in images formed by gravitationallenses that do not possess spherical symmetry. This result demonstrates the utility and power of the angular eikonalmethod.This result (118) is new. It generalizes previous expressions that characterized the deflection of light by thegravitational field of a compact object. These earlier results mostly dealt with the monopole [45, 64, 65], (see [19] andreferences therein), and with the quadrupole contributions [66–68]. Our expression describes the deflection of light inthe presence of an axisymmetric expended body with an arbitrary set of the gravitational multipoles. B. Reducing the double integral using the method of stationary phase
We continue to work with the integral (115). We evaluate one of the two integrals using the method of stationaryphase. We deal with the integral over the magnitude of the impact parameter, b . Introducing θ = ρ/r , we see that(115) has the following explicit form that is useful for practical consideration: γ ( r, θ, φ ) = E ue ik ( r + r g ln 2 kr ) kir π Z π dφ ξ Z ∞ b = R ⋆ ⊙ bdb e ik (cid:0) r (cid:0) b − bρ cos( φ ξ − φ ) (cid:1) − r g (cid:0) ln kb + ψ ( b ) (cid:1)(cid:1) . (119)We recognize that (119) is a double integral with respect to the impact parameter, b : namely, d b = dφ ξ bdb .One may be tempted to try to evaluate this integral using the 2-dimensional method of stationary phase. However,the presence of higher multipoles leads to appearance of caustics so that such a stationary phase solution will notbe accurate, especially at the cusps. Thus, only one of the two integrals in (119) should be approximate using themethod of stationary phase. We choose to approximate the integral over db , leaving the integral over dφ ξ unchanged.Considering (119), we see that the points of stationary phase, where dϕ ( b ) /db = 0, are given by the equation dϕdb = k n r (cid:0) b − ρ cos( φ ξ − φ ) (cid:1) − r g b (cid:16) − ∞ X n =2 J n (cid:16) R ⊙ b (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] (cid:1)(cid:17)o = 0 , (120)which may be transformed as b − bρ cos( φ ξ − φ ) − r g r (cid:16) − ∞ X n =2 J n (cid:16) R ⊙ b (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] (cid:1)(cid:17) = 0 . (121)We solve this equation iteratively, by presenting the impact parameter as b = b [0] + b [1] , where b [0] is the solutioninvolving only the monopole term, while b [1] is due to the eikonal phase. Substituting this trial solution in (121) andequating same orders we get: b [0]2 − b [0] ρ cos( φ ξ − φ ) − r g r = 0 , (122) b [1] (cid:16) b [0] − ρ cos( φ ξ − φ ) (cid:17) + 2 r g r ∞ X n =2 J n (cid:16) R ⊙ b [0] (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] = 0 . (123)From (122), we have b [0] = ± q r g r + (cid:0) ρ cos( φ ξ − φ ) (cid:1) + ρ cos( φ ξ − φ ) . (124)0In the region of strong interference, the relations 0 ≤ θ ≃ p r g /r are satisfied, so that this solution may be givenonly to O ( ρ ). Also, as the magnitude of the impact parameter may only be positive, we choose the positive sign in(124), which yields b [0] = p r g r + ρ cos( φ ξ − φ ) + O ( ρ ) . (125)Substituting this solution into (123), we get b [1] p r g r + r g r ∞ X n =2 J n (cid:16) R ⊙ p r g r (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] = 0 . (126)Thus, we have b [1] = − p r g r ∞ X n =2 J n (cid:16) R ⊙ p r g r (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] + O (cid:0) ρ p r g r J n (cid:1) . (127)Ultimately, we have the following solution for the impact parameter, b = b [0] + b [1] + O ( r g , ρ , ρJ n / p r g r ): b = p r g r + ρ cos( φ ξ − φ ) − p r g r ∞ X n =2 J n (cid:16) R ⊙ p r g r (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] . (128)We compute the stationary phase, ϕ ( b ) from (113) for the values of b given by (128): ϕ ( b ) = k n r g − r g ln k p r g r − r r g r (cid:16) ρ cos( φ ξ − φ ) + p r g r ∞ X n =2 J n n (cid:16) R ⊙ p r g r (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] (cid:17)o . (129)Computing the second derivative of the phase ϕ ( b ) from (120) with respect to b , we have (as we do not account forthe impact of the multipoles on the amplitude of the EM field, we need to know this quantity only to O ( J n )): d ϕdb = k (cid:16) r + 2 r g b + O ( J n ) (cid:17) . (130)Now, using b from (128), to O ( r g , ρ , J n ), we have ϕ ′′ ( b ) = 2 kr (cid:16) − ρ cos( φ ξ − φ ) p r g r (cid:17) ⇒ s π | ϕ ′′ ( b ) | = r πrk (cid:16) ρ cos( φ ξ − φ ) p r g r (cid:17) . (131)We now may compute the amplitude of the integrand in (119), which for b from (128) may be given as A ( b ) s π | ϕ ′′ ( b ) | = b s π | ϕ ′′ ( ℓ ) | = p r g r (cid:16) ρ cos( φ ξ − φ ) p r g r (cid:17)r πrk (cid:16) ρ cos( φ ξ − φ ) p r g r (cid:17) == p r g r r πrk (cid:16) O (cid:16) r g , ρ p r g r , J n (cid:17)(cid:17) . (132)As a result, the expression for the factor γ ( r, θ, φ ) given by (119) takes the form γ ( r, θ, φ ) = E p πkr g e iσ e ikz ×× π Z π dφ ξ exp h − ik r r g r (cid:16) ρ cos( φ ξ − φ ) + p r g r ∞ X n =2 J n n (cid:16) R ⊙ p r g r (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] (cid:17)i , (133)valid to the order of O (cid:0) r g , ρ/ p r g r, ( J /r )Π (cid:1) and the constant σ = − kr g ln( kr g /e ) − π (see [9, 11] for details.)We can present expression (133) in the following compact form: γ ( r, θ, φ ) = E p πkr g e iσ e ikz B ( ρ, φ ) + O (cid:16) r g , ρ p r g r , J r Π (cid:17) , (134)1where B ( ρ, φ ) = B ( x ), with x = ρ (cos φ, sin φ,
0) being the coordinates on the image plane, is the complex amplitudeof the EM field that has the form B ( x ) = 12 π Z π dφ ξ exp h − ik (cid:16)r r g r ρ cos( φ ξ − φ ) + 2 r g ∞ X n =2 J n n (cid:16) R ⊙ p r g r (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] (cid:17)i . (135)The quantity B ( x ) is the complex amplitude of the EM field after it scatters on the gravitational field of an extendedaxisymmetric lens that is represented by a set of gravitational multipoles. If the presence of the gravitational zonalharmonics be neglected, the result (135) reduces to the familiar form for the monopole (see [9] and references therein), J ( k p r g /rρ ). Eq. (135) is a new diffraction integral formula that extends the previous wave-theoretical descriptionof gravitational lensing phenomena to the case of a lens with an arbitrary axisymmetric gravitational potential. Thisresult offers a new, powerful tool to study gravitational lensing in the limit of weak gravitational fields, at the firstpost-Newtonian approximation of the general theory of relativity. C. The EM field in the interference region
Now we are ready to present the components of the EM field in the interference region. The total field in accordwith (84) to O ( r g , r g /r ) has the form ˆ D θ ˆ B θ ! = ˆ B φ − ˆ D φ ! = E p πkr g e iσ B ( ρ, φ ) e i ( kz − ωt ) (cid:18) cos φ sin φ (cid:19) , (136)with radial components of the EM field behaving as ( ˆ D r , ˆ B r ) ≃ O ( r g , r g /r ). The radial components of the EM fieldare negligibly small compared to the other two components, which is consistent with the fact that while passingthrough the gravity field of higher multipoles the EM wave preserves its transverse structure.Expression (136) describes the EM field in the interference region of the SGL in the spherical coordinate system.To study this field in the image plane, we need to transform this result to a cylindrical coordinate system. To dothat, we follow the approach demonstrated in [9], where instead of spherical coordinates ( r, θ, φ ), we introduced acylindrical coordinate system ( ρ, φ, z ), more convenient for these purposes. In the region r ≫ r g , this was done bydefining R = ur = r + r g / O ( r g ) and introducing the coordinate transformations ρ = R sin θ, z = R cos θ , which,from (1), result in the following line element: ds = u − c dt − u (cid:0) dr + r ( dθ + sin θdφ ) (cid:1) = u − c dt − (cid:0) dρ + ρ dφ + nu dz (cid:1) + O ( r g ) . (137)As a result, using (136), for a high-frequency EM wave (i.e., neglecting terms ∝ ( kr ) − ) and for r ≫ r g , we derivethe field near the optical axis, which up to terms of O ( ρ /z ), takes the form (cid:18) E ρ H ρ (cid:19) = (cid:18) H φ − E φ (cid:19) = E p πkr g e iσ B ( ρ, φ ) e i ( kz − ωt ) (cid:18) cos φ sin φ (cid:19) , (138)with ( E z , H z ) = O ( ρ/z ) and where r = p z + ρ = z (1 + ρ / z ) = z + O ( ρ /z )) and θ = ρ/z + O ( ρ /z ).Note that these expressions were obtained using the approximations (102) and are valid for forward scattering when θ ≤ p r g /r , or when ρ ≤ r g . For completeness, one may obtain a more general expression that will be valid for muchlarger deviations from the optical axis, say ρ ∼ R ⊙ . This work is ongoing and results will be reported.Considering the image plane, we see that the quantity B ( ρ, φ ) in (138) given by (135) is a function of the coordinateson the image plane, x = ρ (cos φ, sin φ, B ( ρ, φ ) ≡ B ( x ), is given by a single integral (135).This is our main result. It determines the amplitude of the EM field in the image plane in the strong interferenceregion of the SGL. This function determines the structure of the point-spread function (PSF) of the SGL, whichgoverns the optical properties of the SGL as far as imaging is concerned. This expression describes light from adistant point source, projected onto the image plane by the SGL. Furthermore, it is presented in a form using unitsand parameters that relate directly to physically relevant quantities, making the result readily applicable to studygravitational lensing by real astrophysical objects, such as the Sun.2 a) b)FIG. 3: Even and odd caustics representing individual contributions of the multipoles of a gravitational field to the PSF ofthe extended axisymmetric gravitational lens, obtained through numerical integration of PSF = | B ( x ) | with B ( x ) from (141).Images of a point source formed in the image plane of the lens. From top left, clockwise: a) J , J , J and J ; b) monopole, J , J and J . For the odd-numbered caustics, a change in sign flips the image in the north-south direction.FIG. 4: Interaction between caustics and the effects of sign, calculated by numerically integrating | B ( x ) | using (141). Top rowdepicts the effect of J , distorting the J astroid starting with a negative value similar in magnitude to J , going through 0 andreaching a positive value. Bottom row depicts the effect of J on J in a similar fashion. These images demonstrate that thesign of the J caustic reverses its vertical, “north–south” orientation, whereas the sign of the J caustic determines if it is theastroid’s vertical or horizontal pair of cusps that are “split” as the astroid is stretched in the horizontal vs. vertical direction. D. Multipole contributions
Using the result (138), we may now compute the energy flux in the image region of the lens. With overline andbrackets denoting time-averaging and ensemble averaging, the relevant components of the time-averaged Poyntingvector for the EM field in the image volume may be given in the following form (see [9–11] for details): S z ( x ) = c π (cid:10) [Re E × Re H ] z (cid:11) = c π E πkr g (cid:10)(cid:0) Re (cid:2) B ( x ) e i ( kz − ωt ) (cid:3)(cid:1) (cid:11) , (139)3 FIG. 5: The PSF of the SGL with multipole gravitational moments, obtained by integrating | B ( x ) | using (141), using realisticsolar parameters. These examples show light from a λ = 1 µ m point source, projected by the SGL to 650 AU from the Sun.Left: sin β s = 0 . β s ∼ . ◦ from the solar axis of rotation) in an 8 × β s = 0 .
387 ( β s ∼ . ◦ ) in a 120 ×
120 marea, at 6 cm resolution. with ¯ S ρ = ¯ S φ = 0 for all practical purposes. Defining the light amplification as usual [9–11], µ z ( x ) = S z ( x ) / | S ( x ) | ,where S ( x ) being the Poynting vector carried by a plane wave in a vacuum in a flat space-time, we have the lightamplification factor of the lens that, for short wavelengths (i.e., kr g ≫
1) is given as µ z ( x ) = 2 πkr g | B ( x ) | , (140)with | B ( x ) | = B ( x ) B ∗ ( x ), where B ∗ ( x ) is the complex conjugate of B ( x ), given by (135) that we repeat here forconvenience: B ( x ) = 12 π Z π dφ ξ exp h − ik (cid:16)r r g r ρ cos( φ ξ − φ ) + 2 r g ∞ X n =2 J n n (cid:16) R ⊙ p r g r (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] (cid:17)i . (141)As we see, the amplification of the lens is driven by the factor 2 πkr g in (140). However, in the case of the monopole,the complex amplitude of the EM field (141) was ∝ J ( p r g /rρ ); thus it was reaching its maximum of 1 on the opticalaxis, ρ = 0. In the case of an extended gravitating body, the complex amplitude is given by (141), where B ( x ), ingeneral, is a complex quantity whose magnitude is | B ( x ) | <
1. As we see in Fig. 3, it reaches its maximum value noton the optical axis, but on the caustic that is formed in the image plane. For lenses dominated by the contribution ofa single multipole moment, these caustics acquire the shapes of hypocycloids (e.g., the astroid, characterizing the J quadrupole). However, when several multipole moments are present, their interaction results in more complex shapes;see Fig. 4 for some examples. In general, all the light from an extended source is still there on the image plane, butnow it is scrambled, that for imaging will require deconvolution tools.The quantity | B ( x ) | is the point-spread function (PSF) that characterizes the optical properties of the gravitationallens and can be used to assess its imaging capabilities. The PSF of the lens is extended from the J ( k p r g /rρ ) formfor the monopole lens, discussed in [9] and takes the form of | B ( ρ, φ ) | that now provides a complete descriptionof the intensity distribution in the image plane and accounts for gravitational lensing by an arbitrary axisymmetricgravitational potential.4 FIG. 6: The PSF of the SGL in color, evaluated in multiple wavelengths between 400 and 675 nm in 25 nm increments, eachassigned an approximate RGB color. The parameters used are sin β s = 0 .
05 ( β s ∼ . ◦ ) at 650 AU from the Sun, φ s = 30 ◦ ; a2 × E. Extended Sun contribution to image formation
When applying these results to the SGL, we need to recognize the fact that the Sun axisymmetric rotating bodythat also has north-south symmetry. As such, it will have only even multipole moments J n . The solar multipolemoments are determined using available tracking data from interplanetary spacecraft: J = (2 . ± . × − [69],and J = − . × − , J = − . × − , J = 1 . × − [40]. The deflection of light by these multipolesmay lead to light rays missing the optical axis by many meters, resulting in large caustics on the image plane inthe strong interference region of the SGL. With the contribution from J being the dominant one to consider (seeFig. 5), depending on the target’s position with respect to the solar rotational axis (captured by the angle β s ), someof these multipoles may be needed for developing a comprehensive physical model needed for image deconvolution.The multipole moments of the Sun may also be varying temporally [70], which requires further analysis. On theother hand, the magnitudes of light deflection due to J and higher multipoles are very small at IR, optical or longerwavelengths. These fall within the diffraction pattern of the solar monopole, and thus may be omitted.With these considerations in mind, the most comprehensive form of the complex amplitude of the EM field in thestrong interference region of the SGL is given by (141) were multipole summation is from n = 2 to n = 8, which iscorrect to the order of O ( J ).Clearly, there is no closed-form analytical solution for this integral. It can, however, be readily evaluated usingnumerical methods. It is also clear that, as J , J , J are small, the resulting diffraction pattern will be dominated bythe quadrupole with other multipoles contributing only small corrections (see Fig. 5).The result (141) depends on the wavelength of incident light. However, the geometric shape of the resulting caustic iswavelength-independent. This becomes evident when we evaluate (141) in multiple wavelengths (Fig. 6). Wavelength-dependent features (represented by approximate RGB color in this figure) are clearly evident both inside and outsidethe caustic boundary. However, the caustic boundary’s location does not change, and the cusps, in particular, areachromatic white.Finally, we mention that expression (141) may be easily generalized to the case of extended sources at large butfinite distances from the Sun. Examining this integral, we see that it contains the expression ρ cos( φ ξ − φ ), whichmay be transformed as ρ cos( φ ξ − φ ) = ( n ξ · x ) , where x = ρ (cos φ, sin φ, . (142)The form (141) allows us to extend the new formulation to the case of sources at large but finite distances, z . A5formal way to extend the result (141) to the case of extended source is to rotate the coordinate system by a smallangle (¯ z/z ) x ′ , as discussed in [13], where ¯ z is the heliocentric distance to the image in the strong interference regionof the SGL, z is the heliocentric distance to the source plane and x ′ is a particular point on that source plane. As aresult, to deal with extended sources we start with (141) and extend the argument as follows: x ⇒ x + ¯ zz x ′ , where x ′ = ρ ′ (cos φ ′ , sin φ ′ , , (143)which, equivalently, may be expressed as ( n ξ · x ) → ( n ξ · x ) + (¯ z/z )( n ξ · x ′ ), where n ξ = (cos φ ξ , sin ξ ξ , B ( x , x ′ ) = 12 π Z π dφ ξ exp h − ik nr r g r (cid:16) n ξ · (cid:0) x + ¯ zz x ′ (cid:1)(cid:17) ++ 2 r g ∞ X n =2 J n n (cid:16) R ⊙ p r g r (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] oi . (144)This expression allows us to consider imaging of extended bodies that are positioned at large, but finite distancesfrom the SGL, with the SGL now treated as that produced by a gravitating body that is axisymmetric and rotatingthus admitting characterization of its external gravitational field by zonal harmonics. VI. DISCUSSION AND CONCLUSIONS
This paper represents a continuation of our efforts to provide a reliable, accurate, complete theoretical descriptionof the image formation capabilities of gravitational lenses within the post-Newtonian approximation of the generaltheory of relativity. This work is especially relevant to our on-going work on the study of the optical properties of theSGL in the context of use for a resolved imaging of distant faint sources.In previous papers [8, 9], we offered a complete wave-theoretical description of the SGL under the simplifyingassumption that the Sun’s gravitational field is accurately represented as a gravitational monopole that was modeledas a point mass. Clearly, this is not exactly the case: the actual gravitational field of the Sun deviates from themonopole slightly. Though the effect is very small compared to the size of the solar system, it has considerable impacton the image formation capabilities of the SGL. Therefore, an accurate and complete description of the SGL mustproperly take into account these small deviations from spherical symmetry.It was long understood that the tools of geometric optics are limited when it comes to caustics and the full wave-optical treatment is required [71, 72]. It was in light of this limitation that we developed our new method to describegravitational lensing within the weak-field and slow motion (i.e., a post-Newtonian) approximation of the generaltheory of relativity. The new method addresses light propagation in a weak gravitational field of arbitrary shape, notrestricted by spherical symmetry. Our formalism allows us to describe the contribution of deviations from sphericalsymmetry on the optical properties of corresponding lens using the language of spherical harmonics. In particular,we can use zonal harmonics in particular in the case of an axisymmetric body, such as the Sun.Key to our approach is what we dubbed the angular eikonal method : a convolution of the eikonal phase (which isused to characterize deviations from spherical symmetry) and integral representations of Bessel-functions (that relyon the symmetries that exist in the case of a monopole lens). This allows us to correctly capture the functionaldependence of the integrand and, in effect, to solve the wave equations within a slightly modified symmetry that isextended from spherical to an azimuthally perturbed one (that is due to the presence of the multipole moments.)The method is consistent with the thin lens or eikonal approximations used with in the scalar theory of diffraction[7, 22, 29–31] that are frequently used to describe gravitational lensing.Our method preserves the structure and vector nature of the EM field and allows us to treat the diffracted EM fieldusing regular tools of modern optics [22]. The entire diffraction behavior is captured in the form of a single integral(141), which extends the set of analytical tools developed for gravitational lenses. The approach that we present andand the resulting expressions are applicable to a wide variety of astrophysical lenses. Applying the method to the caseof an axisymmetric body, represented using zonal harmonics, we arrive at our main result, Eq. (144), which reducesthe problem of the finding the EM field in the image plane placed in the string interference region of the SGL. Thissolution preserves the vector nature of the EM field, thus, going beyond the approaches relying on the scalar theories.An important outcomes of the new solution is that it allows one evaluate the behavior of the PSF of an extendedgravitating lens. Applying this method to the SGL, we treat the solar gravitational field as that produced by anaxisymmetric rotating body whose external gravity field is determined by the infinite set of zonal spherical harmonics.6The PSF of the SGL is now determined by a single, well-behaved integral that can be readily evaluated using numericalmethods, especially near the optical axis of the gravitational lens in what we call the region of strong interference.Concerning the imaging of extended sources with the SGL of the extended Sun, we note that the total energydeposited in the image plane is still almost the same as it was in the case of the monopole SGL. However, the PSFof the extended SGL scrambles light on the image plane more than it did in the case of treating the Sun as the pointmass. This will adversely affect the signal-to-noice ratio as far as as the realistic imaging capabilities of the extendedSGL are concerned. The impact on the observing scenario and the integration time are being investigated.Concluding, we note that the new method can be used to investigate image formation processes for extendedsources by the SGL, at a variety of wavelengths, using physically realistic observational scenarios. The approachthat we presented may also be used in reverse: observing astrophysical lensing of distant objects may allow one toreconstruct the multipoles of the gravitational field of the lens and infer its mass distribution, possibly offering a newpractical method in modern astrophysics. Our solution may also help in other areas, such as the modeling of particlecollisions in high energy particle physics experiments on potentials with complex structure. Results of our studies inthese and other directions, once available, will be published elsewhere.
Acknowledgments
This work in part was performed at the Jet Propulsion Laboratory, California Institute of Technology, undera contract with the National Aeronautics and Space Administration. VTT acknowledges the generous support ofPlamen Vasilev and other Patreon patrons. [1] A. Einstein, Annalen der Physik , 146 (1916).[2] A. Einstein, Science , 506 (1936).[3] S. Liebes, Phys. Rev. , B835 (1964).[4] S. Refsdal, MNRAS , 307 (1964).[5] P. S. Schneider, J. Ehlers, and E. Falco, Gravitational Lenses (Springer-Verlag Berlin Heidelberg, 1992).[6] V. A. Fock,
The Theory of Space, Time and Gravitation (Fizmatgiz, Moscow (in Russian), 1959), [English translation(1959), Pergamon, Oxford].[7] L. D. Landau and E. M. Lifshitz,
The Classical Theory of Fields. (7th edition. Nauka: Moscow (in Russian), 1988).[8] S. G. Turyshev, Phys. Rev. D , 084041 (2017), arXiv:1703.05783 [gr-qc].[9] S. G. Turyshev and V. T. Toth, Phys. Rev. D , 024008 (2017), arXiv:1704.06824 [gr-qc].[10] S. G. Turyshev and V. T. Toth, Phys. Rev. D , 024038 (2020), arXiv:2002.06492 [astro-ph.IM].[11] S. G. Turyshev and V. T. Toth, Phys. Rev. D , 024044 (2019), arXiv:1810.06627 [gr-qc].[12] S. G. Turyshev and V. T. Toth, J. Optics , 045601 (2019), arXiv:1805.00398 [physics.optics].[13] S. G. Turyshev and V. T. Toth, Phys. Rev. D , 084018 (2019), arXiv:1908.01948 [gr-qc].[14] S. G. Turyshev and V. T. Toth, Phys. Rev. D , 044025 (2020), arXiv:1909.03116 [gr-qc].[15] S. G. Turyshev and V. T. Toth, Phys. Rev. D , 044048 (2020), arXiv:1911.03260 [gr-qc].[16] V. T. Toth and S. G. Turyshev, Phys. Rev. D (2020), arXiv:2012.05477 [gr-qc].[17] S. G. Turyshev, M. Shao, V. T. Toth, and et al., Direct multipixel imaging and spectroscopy of an exoplanet with a solargravity lens mission (2020), arXiv:1908.01948 [gr-qc].[18] S. G. Turyshev and V. T. Toth, Int. J. Mod. Phys.
D24 , 1550039 (2015), arXiv:1304.8122 [gr-qc].[19] C. M. Will,
Theory and Experiment in Gravitational Physics (Cambridge University Press, Cambridge, UK, 1993).[20] H. Asada and M. Kasai, Prog. Theor. Phys. , 95 (2020).[21] G. Mie, Annalen der Physik , 377 (1908).[22] M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (Cambridge University Press; 7th edition, October 13, 1999).[23] L. I. Schiff,
Quantum mechanics. (McGraw-Hill, 1968).[24] L. D. Landau and E. M. Lifshitz,
Quantum mechanics. Non-Relativistic Theory. (4th edition. Nauka: Moscow (in Russian),1989).[25] A. Messiah,
Quantum Mechanics, Vol 1 (John Wiley & Sons, 1968).[26] H. Friedrich,
Theoretical Atomic Physics, 3-ed (Springer-Verlag, Berlin, Heidelberg, 2006).[27] H. Friedrich,
Scattering Theory (Springer-Verlag, Berlin, Heidelberg, 2013).[28] P. G. Burke,
R-Matrix Theory of Atomic Collisions. Application to Atomic, Molecular and Optical Processes (Springer-Verlag, Berlin, Heidelberg, 2011).[29] S. K. Sharma, T. K. Roy, and D. J. Sommerford, Journal of Physics D: Applied Physics , 1685 (1988).[30] S. K. Sharma and D. J. Sommerford, Il Nuovo Cimento D , 719 (1990). [31] S. K. Sharma and D. J. Sommerford, Light Scattering by Optically Soft Particles: Theory and Applications (Springer-Verlag,Berlin, Heidelberg, New York, 2006).[32] M. D. Semon and J. R. Taylor, Phys. Rev. A , 33 (1977).[33] W. T. Grandy Jr., Scattering of Waves from Large Spheres (Cambridge University Press, 2000).[34] S. M. Kopeikin, J. Math. Phys. , 2587 (1997).[35] S. M. Kopeikin, M. Efroimsky, and G. Kaplan, Relativistic Celestial Mechanics of the Solar System (Wiley-VCH, Berlin,2011).[36] R. P. Feynman, Rev. of Mod. Phys. , 367 (1948).[37] R. Feynman and A. Hibbs, Quantum Mechanics and Path Integrals (McGraw-Hill, New York, 1965).[38] T. T. Nakamura and S. Deguchi, Prog. Theor. Phys. Supp. , 137 (1999).[39] K. Yamamoto, Int. J. Astron. Astrophys. , 221 (2017).[40] I. W. Roxburgh, Astron. Astrophys. , 688 (2001).[41] C. Le Poncin-Lafitte and P. Teyssandier, Phys. Rev. D , 044029 (2008).[42] M. H. Hull Jr. and D. G. Breit, in Handbuch der Physik, vol. 41/1 , edited by S. Flugge (Springer, Berlin, Gottingen,Heidelberg, 1959), p. 408.[43] I. J. Thompson and F. M. Nunes,
Nuclear Reactions for Astrophysics: Principles, Calculation and Applications of Low-Energy Reactions (Cambridge University Press, 2009), 1st ed.[44] S. G. Turyshev and V. T. Toth, Phys. Rev. D , 104015 (2018), arXiv:1805.10581 [gr-qc].[45] E. Herlt and H. Stephani, Int. J. Theor. Phys. , 45 (1976).[46] W. T. Grandy Jr., Scattering of Waves from Large Spheres (Cambridge University Press, Cambridge, UK, 2005).[47] H. C. van de Hulst,
Light Scattering by Small Particles (Dover Publications, New York, 1981).[48] M. Kerker,
The scattering of light, and other electromagnetic radiation (Academic Press, New York, 1969).[49] P. M. Morse and H. Feshbach,
Methods of Theoretical Physics, Part I (McGraw-Hill Science, New York, 1953).[50] S. G. Turyshev and V. T. Toth, Phys. Rev. A , 033810 (2018), arXiv:1801.06253 [physics.optics].[51] P. G. Burke, Potential Scattering in Atomic Physics. (Springer, 2011).[52] R. G. Newton,
Scattering Theory of Waves and Particles (Dover Books on Physics. 2-nd Edition, 2013).[53] H. Bateman, A. Erd´elyi, and Bateman Manuscript Project,
Higher transcendental functions , vol. 1 of
Higher TranscendentalFunctions (McGraw-Hill, New York, 1953).[54] G. A. Korn and T. M. Korn,
Mathematical Handbook for Scientists and Engineers: Definitions, Theorems, and Formulasfor Reference and Review (McGraw-Hill Book Co., New York, 1968).[55] M. Abramowitz and I. A. Stegun,
Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables. (Dover Publications, New York; revised edition, 1965).[56] L. Bakaleinikov and A. Silbergleit, J. Math. Phys. , 083503 (2020).[57] Y. Nambu, in J. Phys. Conf. Ser. (2013), vol. 410 of
J. Phys. Conf. Ser , p. 012036.[58] W. J. Cody and K. E. Hillstrom, Mathematics of Computation , 671 (1970).[59] J. C. A. Barata, L. F. Canto, and M. S. Hussein, Braz. J. Phys. , 50 (2011).[60] P. Schneider, C. Kochanek, and J. Wambsganss, Gravitational Lensing: Strong, Weak and Micro: Saas-Fee AdvancedCourse 33 (Springer, Berlin, 2006).[61] S. Deguchi and W. D. Watson, Astrophys. J. , 440 (1987).[62] Y. Nambu, Int. J. Astron. Astrophys. , 1 (2013).[63] N. Matsunaga and K. Yamamoto, JCAP , 023 (2006).[64] S. Deguchi and W. D. Watson, Astrophys. J. , 30 (1986).[65] Y. Nambu, J. Phys. Conf. Ser. , 012036 (2013).[66] S. A. Klioner, Sov. Astron. , 523 (1991).[67] S. A. Klioner and S. M. Kopeikin, Astron. J. , 897 (1992).[68] S. Zschocke and S. A. Klioner, CQG , 015009 (2010).[69] R. S. Park, W. M. Folkner, A. S. Konopliv, J. G. Williams, D. E. Smith, and M. T. Zuber, Astron. J. , 121 (2017).[70] J. Rozelot and S. Eren, Adv. Space Res. , 2821 (2020).[71] M. V. Berry and C. Upstill, Optics Laser Technology , 257 (1982).[72] M. V. Berry, in Huygens’ Principle 1690-1990: Theory and Applications , edited by H. Block, H.A. Ferwerda, and H.K.Kuiken (Elsevier Science Publishers B.V., 1992), pp. 97–111.[73] L. Blanchet and T. Damour, Philos. Trans. R. Soc. London Ser. A , 379 (1986).[74] L. Blanchet and T. Damour, Ann. Inst. Henri Poincar´e , 377 (1989).[75] B. Linet and P. Teyssandier, Phys. Rev. D (2002).[76] S. Mathis and C. Le Poncin-Lafitte, Astron. & Astrophys. , 889 (2007).[77] K. S. Thorne, Rev. Mod. Phys. , 299 (1980).[78] S. G. Turyshev, V. T. Toth, and M. V. Sazhin, Phys. Rev. D , 024020 (2013), arXiv:1212.0232 [gr-qc].[79] S. G. Turyshev, M. V. Sazhin, and V. T. Toth, Phys. Rev. D , 105029 (2014), arXiv:1402.7111 [gr-qc].[80] S. G. Turyshev, W. Farr, W. M. Folkner, A. R. Girerd, H. Hemmati, T. W. J. Murphy, J. G. Williams, and J. J. Degnan,Exper. Astron. , 209 (2010), arXiv:1003.4961 [gr-qc].[81] N. F. Mott, Proc. Royal Soc. London Ser. A , 542 (1928).[82] W. Gordon, Zeitschrift f¨ur Physik , 180 (1928).[83] L. D. Landau and E. M. Lifshitz, Mechanics. (4th edition. Nauka: Moscow (in Russian), 1988).[84] M. Chaichian and A. Demichev,
Path Integrals in Physics. Volume I: Stochastic Processes and Quantum Mechanics (IOP Publishing, Bristol and Philadelphia, 2001).
Appendix A: Representation of the field in terms of Debye potentials
To represent the EM field equations in terms of Debye potentials, we start with (3)–(4), where we treat gravityto be static, thus ˙ u = 0. Assuming, as usual (we follow closely the discussion presented in [9, 11, 22], adapted forthe gravitational lens), the time dependence of the field in the form exp( − iωt ) where k = ω/c , the time-independentparts of the electric and magnetic vectors satisfy Maxwell’s equations:curl D = ik u B + O ( G ) , (A1)curl B = − ik u D + O ( G ) . (A2)As shown in [9], in spherical coordinates, the field equations (A1)–(A2) to O ( G ) become − ik u ˆ D r = 1 r sin θ (cid:16) ∂∂θ ( r sin θ ˆ B φ ) − ∂∂φ ( r ˆ B θ ) (cid:17) , (A3) − ik u ˆ D θ = 1 r sin θ (cid:16) ∂ ˆ B r ∂φ − ∂∂r ( r sin θ ˆ B φ ) (cid:17) , (A4) − ik u ˆ D φ = 1 r (cid:16) ∂∂r ( r ˆ B θ ) − ∂ ˆ B r ∂θ (cid:17) , (A5) ik u ˆ B r = 1 r sin θ (cid:16) ∂∂θ ( r sin θ ˆ D φ ) − ∂∂φ ( r ˆ D θ ) (cid:17) , (A6) ik u ˆ B θ = 1 r sin θ (cid:16) ∂ ˆ D r ∂φ − ∂∂r ( r sin θ ˆ D φ ) (cid:17) , (A7) ik u ˆ B φ = 1 r (cid:16) ∂∂r ( r ˆ D θ ) − ∂ ˆ D r ∂θ (cid:17) , (A8)while the remaining two equations from Eq. (3)–(4) take the form ∂∂r (cid:16) u r sin θB r (cid:17) + ∂∂θ (cid:16) u r sin θB θ (cid:17) + ∂∂φ (cid:16) u rB φ (cid:17) = 0 , (A9) ∂∂r (cid:16) u r sin θD r (cid:17) + ∂∂θ (cid:16) u r sin θD θ (cid:17) + ∂∂φ (cid:16) u rD φ (cid:17) = 0 . (A10)Our goal is to find a general solution to these equations in the form of a superposition of two linearly independentsolutions (cid:0) e D , e B (cid:1) and (cid:0) m D , m B (cid:1) that satisfy the following relationships: e ˆ D r = ˆ D r , e ˆ B r = 0 , (A11) m ˆ D r = 0 , m ˆ B r = ˆ B r . (A12)With ˆ B r = e ˆ B r = 0, (A4) and (A5) become ik u e ˆ D θ = 1 r ∂∂r (cid:0) r e ˆ B φ (cid:1) , (A13) ik u e ˆ D φ = − r ∂∂r (cid:0) r e ˆ B θ (cid:1) . (A14)Substituting these relationships into (A7) and (A8) we obtain ∂∂r h u ∂∂r (cid:0) r e ˆ B θ (cid:1)i + k u ( r e ˆ B θ ) = − ik sin θ ∂ e ˆ D r ∂φ , (A15) ∂∂r h u ∂∂r (cid:0) r e ˆ B φ (cid:1)i + k u ( r e ˆ B φ ) = ik ∂ e ˆ D r ∂θ . (A16)9From div( u e B ) = 0 given by Eq. (4) (which in the expanded form is given by Eq. (A9)) and using our assumptionthat e ˆ B r = 0, we have ∂∂θ (cid:0) u sin θ e ˆ B θ (cid:1) + ∂∂φ (cid:0) u e ˆ B φ (cid:1) = 0 , (A17)which ensures that (A6) is also satisfied at the needed level of accuracy. As we know, this equation is valid for aspherically symmetric gravitational field. Terms that characterize deviations from the monopole in the generic formof the Newtonian potential, U , lack spherical symmetry. For these terms, the condition (A17) may be satisfied onlyapproximately. Indeed, after substitution from (A13), (A14), in (A6), we have1 r sin θ (cid:16) ∂∂θ (cid:0) r sin θ e ˆ D φ (cid:1) − ∂∂φ (cid:0) r e ˆ D θ (cid:1)(cid:17) == − ik u r sin θ n ∂∂r h ru (cid:16) ∂∂θ (cid:0) u sin θ e ˆ B θ (cid:1) + ∂∂φ (cid:0) u e ˆ B φ (cid:1)(cid:17)i ++ 2 ∂∂r h r (cid:16) sin θ e ˆ B θ ∂ ln u ∂θ + e ˆ B φ ∂ ln u ∂φ (cid:17)i − r (cid:16) sin θ e ˆ B θ ∂ ln u ∂r∂θ + e ˆ B φ ∂ ln u ∂r∂φ (cid:17)o == 1 ik u r sin θ n ∂∂r h r (cid:16) sin θ e ˆ B θ ∂ ln u ∂θ + e ˆ B φ ∂ ln u ∂φ (cid:17)i − r (cid:16) sin θ e ˆ B θ ∂ ln u ∂r∂θ + e ˆ B φ ∂ ln u ∂r∂φ (cid:17)o . (A18)The first term in this expression vanishes because of (A17) (as it was in the case of a monopole, see [9]). Consideringthe remaining terms, and taking into account the form of u from (2) with the Newtonian potential, U , given either(B5) or (B4), we see that the following relation is true: ∂∂θ (cid:0) r sin θ e ˆ D φ (cid:1) − ∂∂φ (cid:0) r e ˆ D θ (cid:1) == 1 ik u n ∂∂r h r (cid:16) sin θ e ˆ B θ ∂ ln u ∂θ + e ˆ B φ ∂ ln u ∂φ (cid:17)i − r (cid:16) sin θ e ˆ B θ ∂ ln u ∂r∂θ + e ˆ B φ ∂ ln u ∂r∂φ (cid:17)o ≃ ik u O (cid:16) J R ⊙ r (cid:17) , (A19)where J is the gravitational quadrupole moment of the mass distribution inside the lens (which typically is the largestterm after the monopole). Clearly, as r gets larger, this expression vanishes, justifying the validity of (A6), namelylim r →∞ r sin θ (cid:16) ∂∂θ (cid:0) r sin θ e ˆ D φ (cid:1) − ∂∂φ (cid:0) r e ˆ D θ (cid:1)(cid:17) ≃ ik u r sin θ O (cid:16) J R ⊙ r (cid:17) → . (A20)In all practical scenarios, the limit (A20) is satisfied for r & R ⊙ . Thus, for scenarios relevant for the SGL, (A18) isequal to 0. As a result, when describing the SGL and considering light propagation in a weak gravitational field, wemay neglect the effect of the gravitational field on the amplitude of the EM wave. In this case, our primary interest isthe phase of the wave, thus this expression constitutes the condition consistent with the eikonal approximation. Thecomplementary case with m ˆ D r = 0 is treated identically, in accordance with (A12).When the radial magnetic field vanishes, the solution is called the electric wave (or transverse magnetic wave);correspondingly, when the radial electric field vanishes, the solution is called the magnetic wave (or transverse electricwave). These can both be derived from the corresponding Debye scalar potentials e Π and m Π.Given e ˆ B r = 0, e ˆ D φ and e ˆ D θ in (A6) can be represented as a scalar field’s gradient: e ˆ D φ = 1 r sin θ ∂U∂φ + 1 ik u r sin θ O (cid:16) J R ⊙ r (cid:17) , e ˆ D θ = 1 r ∂U∂θ + 1 ik u r O (cid:16) J R ⊙ r (cid:17) , (A21)where U is some function. Introducing the electric Debye potential e Π that relates to U as U = 1 u ∂∂r (cid:0) r e Π (cid:1) . (A22)We use this expression in (A21) and obtain e ˆ D θ = 1 u r n ∂ (cid:0) r e Π (cid:1) ∂r∂θ − ∂ ln u ∂θ ∂∂r (cid:0) r e Π (cid:1)o , e ˆ D φ = 1 u r sin θ n ∂ (cid:0) r e Π (cid:1) ∂r∂φ − ∂ ln u ∂φ ∂∂r (cid:0) r e Π (cid:1)o , (A23)which satisfy Eq. (A6) with e ˆ B r = 0 and thus also (A19).0It can be seen that (A13) and (A14) are satisfied by e ˆ B φ = ikr ∂ (cid:0) r e Π (cid:1) ∂θ + O (cid:16) J r e Π (cid:17) , e ˆ B θ = − ikr sin θ ∂ (cid:0) r e Π (cid:1) ∂φ + O (cid:16) J r e Π (cid:17) . (A24)If we substitute both of (A24) into (A3) we obtain e ˆ D r = − u r sin θ h ∂∂θ (cid:16) sin θ ∂ ( r e Π) ∂θ (cid:17) + 1sin θ ∂ ( r e Π) ∂φ i + O (cid:16) J r e Π (cid:17) . (A25)Substituting expressions (A24) into (A15)–(A16) yields − ik sin θ ∂∂φ n ∂∂r h u ∂∂r (cid:0) r e Π (cid:1)i + k u ( r e Π) − e ˆ D r o = O (cid:16) J r e Π (cid:17) , (A26) ik ∂∂θ n ∂∂r h u ∂∂r (cid:0) r e Π (cid:1)i + k u ( r e Π) − e ˆ D r o = O (cid:16) J r e Π (cid:17) , (A27)i.e., the derivative of the same expression with respect to both φ and θ vanishes. This is clearly satisfied if we set theexpression itself to O (cid:0) ( J /r ) e Π (cid:1) . Dividing by u and using (A25) leads to1 u ∂∂r h u ∂ ( r e Π) ∂r i + 1 u r sin θ n ∂∂θ (cid:16) sin θ ∂ ( r e Π) ∂θ (cid:17) + 1sin θ ∂ ( r e Π) ∂φ o + k ( r e Π) = O (cid:16) J r e Π (cid:17) . (A28)Defining u ′ = ∂u/∂r and u ′′ = ∂ u/∂r , this equation may be rewritten as1 r ∂∂r (cid:16) r ∂∂r h e Π u i(cid:17) + 1 r sin θ n ∂∂θ (cid:16) sin θ ∂∂θ h e Π u i(cid:17) + 1sin θ ∂ ∂φ h e Π u io + (cid:16) k u − u (cid:0) u (cid:1) ′′ (cid:17)h e Π u i = O (cid:16) J r e Π u (cid:17) , (A29)which is the wave equation for the quantity e Π /u : (cid:16) ∆ + k u − u (cid:0) u (cid:1) ′′ (cid:17)h e Π u i = O (cid:16) r g , J r e Π u (cid:17) . (A30)We are concerned only with the field produced by the extended gravitational field. Thus, the quantity u has thefrom u ( r ) = 1 + U/c + O ( c − ) , as given by (2). With this, we can rewrite (A30) as (cid:16) ∆ + k (cid:0) Uc (cid:1) − u (cid:0) u (cid:1) ′′ (cid:17)h e Π u i = O (cid:16) r g , J r e Π u (cid:17) . (A31)Equation (A31) is similar to the Schr¨odinger equation of quantum mechanics, used to describe scattering on theCoulomb potential. However, this equation has an extra potential of − u ( u − ) ′′ ≃ r g /r . It is known [25] that thepresence of potentials of ∝ /r in (A31) does not alter the asymptotic behavior of the solutions. Neglecting the u ( u − ) ′′ ≃ r − term in (A31) reduces this equation to the form of a time-independent Schr¨odinger equation thatdescribes scattering in a Newtonian potential: (cid:16) ∆ + k (cid:0) Uc (cid:1)(cid:17)h e Π u i = O (cid:16) r g , J r e Π u (cid:17) . (A32)In the case of the SGL, we will always be at distances that are much larger than the Sun’s Schwarzschild radius.Therefore, we may neglect the term u ( u − ) ′′ ≃ r g /r in (A31). We use (A32) for the purposes of establishing theproperties of the EM wave diffracted by the solar gravitational lens. An identical equation may be obtained for m Π.This solution is consistent with the eikonal approximation, the use of which to describe the scattering of high-energyparticles or processes related to the diffraction of light is well-justified.By means of (A28), Eq. (A25) may be written as e ˆ D r = ∂∂r h u ∂ ( r e Π) ∂r i + k u ( r e Π) + O (cid:16) J r e Π u (cid:17) . (A33)It can be verified by substituting (A23)–(A28) and (A33) into (A3)–(A8) that we have obtained a solution of our setof equations. In a similar way, we may consider the magnetic wave. We find that this wave can be derived from apotential m Π which satisfies the same differential equation (A28) as e Π.1The complete solution of the EM field equations is obtained by adding the two fields (as discussed in [21, 22, 48]),namely D = e D + m D ; and B = e B + m B . This yields, to O (cid:0) r g , ( J /r )( e Π /u ) (cid:1) ,ˆ D r = 1 u n ∂ ∂r h r e Π u i + (cid:16) k u − u (cid:0) u (cid:1) ′′ (cid:17)h r e Π u io , (A34)ˆ D θ = 1 u r ∂ (cid:0) r e Π (cid:1) ∂r∂θ + ikr sin θ ∂ (cid:0) r m Π (cid:1) ∂φ , (A35)ˆ D φ = 1 u r sin θ ∂ (cid:0) r e Π (cid:1) ∂r∂φ − ikr ∂ (cid:0) r m Π (cid:1) ∂θ , (A36)ˆ B r = 1 u n ∂ ∂r h r m Π u i + (cid:16) k u − u (cid:0) u (cid:1) ′′ (cid:17)h r m Π u io , (A37)ˆ B θ = − ikr sin θ ∂ (cid:0) r e Π (cid:1) ∂φ + 1 u r ∂ (cid:0) r m Π (cid:1) ∂r∂θ , (A38)ˆ B φ = ikr ∂ (cid:0) r e Π (cid:1) ∂θ + 1 u r sin θ ∂ (cid:0) r m Π (cid:1) ∂r∂φ . (A39)Both potentials Π = (cid:0) e Π , m Π (cid:1) are solutions of the differential equation (A30), which, in the case of the weak gravitycharacteristic for the SGL, is given by (A32) and in terms of potential Π takes the form: (cid:16) ∆ + k (cid:0) Uc (cid:1)(cid:17)h Π u i = O (cid:16) r g , J r Π u (cid:17) . (A40)This completes decomposition of the Maxwell equations (3)–(4) on the curved background in the weak gravitationalfield of the solar system. Eqs. (A34)–(A39) together with (A40) is our primary result that will use throughout thispaper. Again, note that in (A40), we discarded the term u ( u − ) ′′ ∼ /r , representing the tail of the gravitationalpotential, as insignificant (see discussions in Appendix F of Ref. [9] and in Appendix C of [11]).Finally, for the components ˆ D θ , ˆ D φ and ˆ B θ , ˆ B φ to be continuous over a spherical surface at some large distancefrom the origin, r = R ⋆ , it is evidently sufficient that the four quantities [9] ǫ ( r e Π) , µ ( r m Π) , ∂ ( r e Π) ∂r , ∂ ( r m Π) ∂r , (A41)shall also be continuous over this surface. Thus, our boundary conditions also split into independent conditions on e Πand m Π. Our problem is thus reduced to the problem of finding two mutually independent solutions of the equations(A28) with prescribed boundary conditions.
Appendix B: Computing the eikonal phase1. Different forms of the gravitational potential
Before we proceed with solving (A40), we recognize that the gravitational potential U from (2) in spherical coor-dinates ( r ≡ | x | , φ, θ ) may be given in the most general case in the form of spherical harmonics: U = G Z σ ( t, x ′ ) d x ′ | x − x ′ | + O ( c − ) = GMr (cid:16) ∞ X ℓ =2 + ℓ X k =0 (cid:16) Rr (cid:17) ℓ P ℓk (cos θ )( C ℓk cos kφ + S ℓk sin kφ ) (cid:17) + O ( c − ) , (B1)where σ ( t, x ′ ) is the relativistic mass density inside the body, M is its mass, R is its radius, P ℓk are the Legendrepolynomials, while C ℓk and S ℓk are relativistic normalized spherical harmonic coefficients that characterize the body.In the case of an axisymmetric body (i.e., the Sun), its external gravitational potential is reduced to the k = 0zonal harmonics and may be expressed [40, 41] in terms of the usual dimensionless multipole moments J ℓ : U = GMr n − ∞ X ℓ =2 J ℓ (cid:16) R ℓ r (cid:17) ℓ P ℓ (cid:16) k · x r (cid:17)o + O ( c − ) , (B2)2where k denotes the unit vector along the x -axis, P ℓ are the Legendre polynomials and the quantities M, J ..., J ℓ correspond to the generalized Blanchet–Damour mass multipole moments in general relativity [73, 74]. Furthermore,in the case of an axisymmetric and rotating body with “north-south symmetry”, such as the Sun, the expression (B2)contains only the ℓ = 2 , , , ... even moments [40].Following [75], we take into account the identity ∂ ℓ ∂z ℓ (cid:16) r (cid:17) = ( − ℓ ℓ ! r ℓ P ℓ (cid:16) k · x r (cid:17) , z = x , (B3)and present U as the following expansion in a series of derivatives of 1 /rU = GM n r − ∞ X ℓ =2 ( − ℓ ℓ ! J ℓ R ℓ ∂ ℓ ∂z ℓ (cid:16) r (cid:17)o + O ( c − ) . (B4)As we shall see below, this form is much more convenient for the computation of integrals involving U .Considering the generic case, it was shown [76] that the scalar gravitational potential (B1) may equivalently begiven in the following form: U = GM n r + ∞ X ℓ =2 ( − ℓ ℓ ! T ∂ ℓ ∂x (cid:16) r (cid:17)o + O ( c − ) , (B5)where r = | x | , M is the post-Newtonian mass of the body, and T are the symmetric trace-free (STF) massmultipole moments of the body [34, 73, 74, 77] defined as M = Z V d x σ ( x ) , T = 1 M Z V d x σ ( x ) x , (B6)where x = x , the angle brackets < ... > denote the STF operator, and V means the total volumeof the isolated gravitating body under consideration. The dipole moment T a is absent in the expansion (B5) sincewe took the origin of the coordinates to be at the center of mass of the body.
2. Computing the eikonal phase
Based on the form of the post-Newtonian potentials (B1), (B5), and (B4), it is convenient to separate the monopoleterm from the rest of the multipoles. As we know [9], the action of the monopole term is similar to that of theCoulomb potential which is a long-range potential that is felt as far as the source. The remaining multipole termsform the short-range potential, V sr c , yielding the decomposition U/c = r g / r + V sr for the Newtonian potential,which allows to present the potential term in (A40) in the following form:4 Uc = 2 r g r + 4 V sr . (B7)The short-range potential forms the eikonal phase given by (39) that has the form ξ b ( τ ) = k Z τ V sr ( b , τ ′ ) dτ ′ + O ( r g ) . (B8) a. Computing the eikonal phase for an axisymmetric body Here we develop an expression for the eikonal phase in the case of an axisymmetric body, with its potential givenby (B4). In this case, the decomposition of the post-Newtonian potential takes the from4 Uc = 2 r g r − r g ∞ X ℓ =2 ( − ℓ ℓ ! J ℓ R ℓ ∂ ℓ ∂s ℓ (cid:16) r (cid:17) . (B9)In the case the short-range potential, V sr from (B5) is given as V sr = − r g ∞ X ℓ =2 ( − ℓ ℓ ! J ℓ R ℓ ∂ ℓ ∂s ℓ (cid:16) r (cid:17) . (B10)3We now compute the leading term of this expansion. For that, we define the vector s to be a unit vector in thedirection of the axis of rotation. Remembering that r = √ b + τ + O ( r g ) from (32), we evaluate directional derivatives ∂/∂s along s ≡ k , which have the form ∂∂s = ( s · ∇ ) = (cid:16) s · ∂∂ r (cid:17) . (B11)This relation allows us to compute the relevant partial derivatives for the leading tersm in (B4): ∂∂s r = − ( s · r ) r , ∂ ∂s r = 3( s · r ) r − r , ∂ ∂s r = − (cid:16) s · r ) r − s · r ) r (cid:17) , (B12) ∂ ∂s r = 3 (cid:16) s · r ) r − s · r ) r + 3 r (cid:17) , ∂ ∂s r = − (cid:16) s · r ) r − s · r ) r + 15( s · r ) r (cid:17) , (B13) ∂ ∂s r = 45 (cid:16) s · r ) r − s · r ) r + 105( s · r ) r − r (cid:17) , (B14) ∂ ∂s r = − (cid:16) s · r ) r − s · r ) r + 315( s · r ) r − s · r ) r (cid:17) , (B15) ∂ ∂s r = 315 (cid:16) s · r ) r − s · r ) r + 6930( s · r ) r − s · r ) r + 35 r (cid:17) . (B16)Using these expressions in (B10), we have2 V sr ( b , τ ) = − r g n J R (cid:16) s · r ) r − r (cid:17) + J R (cid:16) s · r ) r − s · r ) r (cid:17) ++ J R (cid:16) s · r ) r − s · r ) r + 3 r (cid:17) + J R (cid:16) s · r ) r − s · r ) r + 15( s · r ) r (cid:17) ++ J R (cid:16) s · r ) r − s · r ) r + 105( s · r ) r − r (cid:17) ++ J R (cid:16) s · r ) r − s · r ) r + 315( s · r ) r − s · r ) r (cid:17) ++ J R (cid:16) s · r ) r − s · r ) r + 6930( s · r ) r − s · r ) r + 35 r (cid:17) ++ ∞ X ℓ =9 ( − ℓ ℓ ! J ℓ R ℓ ∂ ℓ ∂s ℓ (cid:16) r (cid:17)o . (B17)We can now substitute result (B17) into expression (B8) and integrate it. We observe that, technically, it is morestraightforward to compute the eikonal phase shift integration along the entire path from τ to τ . Note that this wayone computes the double shift, 2 ξ b ( τ ). This integration results in many terms that depend on the distance to thesource, r = p b + τ , and that to the image plane, r = √ b + τ . The resulting expression is greatly simplified inthe case when both the source and the observer on the image plane are located at very large distances from the lensand the following inequalities are satisfied: b/ √ b + τ ≃ b/τ ≪ b/ p b + τ ≃ b/τ ≪
1. This step, essentially,constitutes the thin lens approximation. This allows us to greatly simplify the result of the integration, yielding ξ b ( τ ) = − kr g n J (cid:16) Rb (cid:17) h s · b ) b + ( s · k ) − i + J (cid:16) Rb (cid:17) h ( s · b ) b (cid:16) s · b ) b + 3( s · k ) − (cid:17)i ++ J (cid:16) Rb (cid:17) h(cid:0) ( s · b ) b + ( s · k ) − (cid:1) ( s · b ) b + (cid:0) ( s · k ) − (cid:1) i ++ J (cid:16) Rb (cid:17) h ( s · b ) b (cid:16) ( s · b ) b + (cid:0) ( s · k ) − (cid:1) ( s · b ) b + 5 (cid:0) ( s · k ) − (cid:1) (cid:17)i ++ J (cid:16) Rb (cid:17) h ( s · b ) b + (cid:0) ( s · k ) − (cid:1) ( s · b ) b + (cid:0) ( s · k ) − (cid:1) ( s · b ) b + (cid:0) ( s · k ) − (cid:1) i + If needed, one can use all those terms to evaluate the eikonal phase, ξ b ( τ ), for shorter distances, when τ ∼ τ ≃ b . For problems relatedto gravitational lensing this is unnecessary, but may be needed for some solar system spacecraft tracking applications [78–80]. J (cid:16) Rb (cid:17) h ( s · b ) b (cid:16) ( s · b ) b + (cid:0) ( s · k ) − (cid:1) ( s · b ) b + (cid:0) ( s · k ) − (cid:1) ( s · b ) b + 7 (cid:0) ( s · k ) − (cid:1) (cid:17)i ++ J (cid:16) Rb (cid:17) h ( s · b ) b + (cid:0) ( s · k ) − (cid:1) ( s · b ) b + (cid:0) ( s · k ) − (cid:1) ( s · b ) b ++ (cid:0) ( s · k ) − (cid:1) ( s · b ) b + (cid:0) ( s · k ) − (cid:1) i ++ ∞ X ℓ =9 ( − ℓ ℓ ! J ℓ R ℓ Z ττ ∂ ℓ ∂s ℓ (cid:16) r (cid:17) dτ ′ o + O ( r g ) . (B18)Note that a similar result for the quadrupole J term was obtained in [66–68]. Expression (B18) extends all theprevious computations to the higher order terms including J .We use the heliocentric spherical coordinate system and define the vectors of impact parameter, b , the wavevector, k , and the unit vector long the solar rotational axis, s , as follows: b = b (cos φ ξ , sin φ ξ , , k = (0 , , , s = (sin β s cos φ s , sin β s sin φ s , cos β s ) . (B19)With these definitions, the eikonal phase (B24) for the case of axisymmetric body whose gravitational potential isgiven by (B4) as ξ b ( r ) = − kr g n J (cid:16) R ⊙ b (cid:17) sin β s cos[2( φ ξ − φ s )] + J (cid:16) R ⊙ b (cid:17) sin β s cos[3( φ ξ − φ s )] ++ J (cid:16) R ⊙ b (cid:17) sin β s cos[4( φ ξ − φ s )] + J (cid:16) R ⊙ b (cid:17) sin β s cos[5( φ ξ − φ s )] ++ J (cid:16) R ⊙ b (cid:17) sin β s cos[6( φ ξ − φ s )] + J (cid:16) R ⊙ b (cid:17) sin β s cos[7( φ ξ − φ s )] ++ J (cid:16) R ⊙ b (cid:17) sin β s cos[8( φ ξ − φ s )] + ∞ X n =9 n J n (cid:16) R ⊙ b (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] o + O ( r g ) . (B20)Assuming that the pattern evident in these expressions continues for higher multipoles, we obtain the followingcompact expression for the eikonal phase: ξ b ( b , s ) = − kr g ∞ X n =2 J n n (cid:16) R ⊙ b (cid:17) n sin n β s cos[ n ( φ ξ − φ s )] + O ( r g ) . (B21)Note that the sum in (B21) contains contributions from all multipole moments, n = 2 , , , ... and is valid for anyaxisymmetric body with respect to the z = x axis represented by s . If in addition to being axisymmetric, thatbody also has “north-south” symmetry (symmetry under a reflection with respect to the plane of rotation), that sumcontains only even terms, n = 2 , , , ..., [40]. b. Computing the eikonal phase using STF tensors Using the representations (B1), (B5) or (B4), it is convenient to separate the monopole term in the potential U . Infact, to determine the solution to (A40), similarly to [11, 12], we first separate the monopole contribution and presentthe U -dependent term in (A40) as4 Uc = 2 r g r + 2 r g ∞ X ℓ =2 ( − ℓ ℓ ! T ∂ ℓ ∂x (cid:16) r (cid:17) , (B22)where r g = 2 GM/c is the Schwarzschild radius of the body and the short-range potential V sr from (B5) is given by V sr = r g ∞ X ℓ =2 ( − ℓ ℓ ! T ∂ ℓ ∂x (cid:16) r (cid:17) . (B23) To derive the results in a compact form we used multiple angle formulae: https://mathworld.wolfram.com/Multiple-AngleFormulas.html and also V sr ( r ) from (B23), we reduced the problem to evaluating a single integral to determine the Debye potentialΠ( r ) from (25), which is a great simplification. Given the fact that b is constant and by taking the short-rangepotential V sr ( r ) from (B23), we evaluate (39) as ξ b ( r ) = kr g ∞ X ℓ =2 ( − ℓ ℓ ! T Z ττ ∂ ℓ ∂x (cid:16) r (cid:17) dτ ′ . (B24)In fact, we may generalize expression ∇ = ∇ b + k d/dτ + O ( r g ) and write ∂ ℓ ∂x ≡ ∇ = ℓ X p =0 ℓ ! p !( ℓ − p )! k ∂ p ∂τ p + O ( r g ) , (B25)where a new shorthand notation ∂ a ≡ ∂/∂ b a has been used and τ is defined by (31).Using representation (B25), we can compute the relevant integral (with r = √ b + τ and r = p b + τ , asdiscussed in Sec. II C) Z ττ ∂ ℓ ∂x (cid:16) r (cid:17) dτ ′ = Z ττ ℓ X p =0 ℓ ! p !( ℓ − p )! k ∂ p ∂τ ′ p (cid:16) √ b + τ ′ (cid:17) dτ ′ == k ln r + τr + τ + ℓ X p =1 ℓ ! p !( ℓ − p )! k h ∂ p − ∂τ p − (cid:16) r (cid:17) − ∂ p − ∂τ p − (cid:16) r (cid:17)i . (B26)As a result, the eikonal phase (B24) takes the form: ξ b ( r ) = kr g ∞ X ℓ =2 ( − ℓ ℓ ! T n k ln r + τr + τ ++ ℓ X p =1 ℓ ! p !( ℓ − p )! k h ∂ p − ∂τ p − (cid:16) r (cid:17) − ∂ p − ∂τ p − (cid:16) r (cid:17)io + O ( r g ) . (B27)Note that expression (B27) is derived for an arbitrary matter distribution within a gravitating body. One may usethis result to derive the eikonal phase for any of the terms in the Newtonian potential (B1). Such a generic potentialmay be suitable for analysis of the gravitational lensing by a lens with an arbitrary intrinsic matter distribution. Appendix C: Using the path integral formalism
As we mentioned before, Eq. (12) is nearly identical to the time-independent Schr¨odinger equation that in nucleardescribes the scattering problem on a potential U [25, 81, 82]. Here we further explore this analogy. For that, wenote that in the absence of the scattering potential, U , solution (12) may be given in the from of a plane EM wavegiven as ψ ( r ) = E e i k · r . Next, we introduce a cylindrical coordinate system ( ρ, φ, z ), whose z -axes is directed alongthe wavevector k . Then, by defining the amplification factor due to scattering on the gravitational potential U of thelens as µ ( r ) = [Π( r ) /u ] /ψ ( r ), we rewrite (12) as n ∆ ψ ( r ) + k ψ ( r ) o µ ( r ) + ψ ( r )∆ µ ( r ) + 2 (cid:0) ∇ ψ ( r ) · ∇ µ ( r ) (cid:1) + k U ( r ) c µ ( r ) ψ ( r ) = 0 . (C1)As ψ ( r ) is the solution of the homogeneous wave equation in flat, vacuum spacetime, the first term in (C1) is zero.Then, we can divide the remaining terms of (C1) by ψ ( r ), which yields:∆ µ ( r ) + 2 (cid:16) ∇ ln ψ ( r ) · ∇ µ ( r ) (cid:17) + k U ( r ) c µ ( r ) = 0 . (C2)Clearly, ∇ ln ψ ( r ) = ik k , where k = ω/c is the wavenumber and k is the unit vector in the direction of the wavevector.From discussion in Section II C, we know that ( k · ∇ ) = d/dτ . We remember the form of the Laplacian in the cylindricalcoordinate system ( ρ , z ) that in our case is given as∆ µ ( r ) = ∆ ρ µ ( r ) + ∂ µ ( r ) ∂z , where ∆ ρ µ ( r ) = 1 ρ ∂∂ρ (cid:16) ρ ∂µ ( r ) ∂ρ (cid:17) + 1 ρ ∂ µ ( r ) ∂φ . (C3)6Substituting these results into (C2), we have ∂ µ ( r ) ∂z + ∆ ρ µ ( r ) + 2 ik dµ ( r ) dτ + k U ( r ) c µ ( r ) = 0 . (C4)We assume that k/ | ∂ ln µ/∂z | ∼ (scale at which µ varies)/(wavelength) ≫
1, we neglect the first term compared withthe second term, which constitutes the eikonal approximation. Then Eq. (C4) takes a familiar form: i dµ ( τ, ρ ) dτ = (cid:16) − k ∆ ρ − kU ( τ, ρ ) c (cid:17) µ ( τ, ρ ) , (C5)which is the Schr¨odinger equation with the “time” coordinate τ , the “particle mass” k , and the “time-dependentpotential” − kc − U ( τ, ρ ). The corresponding Lagrangian that yields the classical equation of motion is given as L ( τ, ρ , ˙ ρ ) = k (cid:16) ˙ ρ + 2 c − U ( τ, ρ ) (cid:17) , (C6)where ˙ ρ = d ρ /dτ . This Lagrangian describes the motion of a pendulum in a gravity field with potential − kc − U ( τ, ρ ).This is a mechanical analogy for forming the caustics on the image plane of the SGL. This is similar to a motion ofa connected pendulum where each of the multipoles characterized by a unique natural spacial frequency, affects themotion of the entire pendulum in a carefully prescribed fashion [83].In the path integral formulation [37], the solution to (C5) may formally be written as µ ( r ) = Z D ρ ( τ ) exp h i Z ττ L ( τ, ρ , ˙ ρ ) dτ i . (C7)Following the established rules of evaluating path integrals [37, 84], we have Z ττ k ˙ ρ dτ = k τ (cid:16) ρ ( τ ) − ρ (0) (cid:17) ≃ k r ( b − r θ ) , (C8)with ρ ( τ ) = r θ , ρ (0) = b and where we realize that for very small angles θ = ρ/r , τ = ( k · x ) ≃ r is a validapproximation. Integrating (C8), we also use a thin lens approximation while assuming that the effect of the lenson light is instantaneous and affects light only after it has passed through the lens. Initially the light continues on astraight line, so that ρ (0) − ρ ( τ ) = 0, then, there is a sudden path change after which the light continues on a differentstraight line until it reaches the observer on the image plane. Integrating the potential terms, we use representation(23), that is, 4 U /c = 2 r g /r + 4 V sr and the approach presented in Appendix B 2 a, which yeilds Z ττ k c − U ( τ, ρ ) dτ = kr g ln 2 kr + kr g ln 2 kr − kr g (cid:0) ln kb + ψ ( b ) (cid:1) , (C9)where ψ ( b ) is given by (109). As a result, after applying the appropriate normalization factor A = p k/ πir to eachof the two dimensions involved [37, 84], the expression (C7) results in µ ( r ) = E e ikr g ln 2 kr kir π Z d b exp h ik (cid:16) r ( b − r θ ) − r g (cid:0) ln kb + ψ ( b ) (cid:1)(cid:17)i . (C10)Combining this expression with ψ ( r ) = E e i k · r , we get for the Debye potential Π( r ) an expression that is equivalentto (115) for the factor γ ( r, θ, φr, θ, φ