Diffraction of electromagnetic waves in the gravitational field of the Sun

TL;DR

This paper develops an exact wave-theoretical model of electromagnetic wave diffraction in the Sun's gravitational field, providing insights into the solar gravitational lens's optical properties for potential exoplanet imaging.

Contribution

It presents a closed-form solution for EM wave diffraction in the Sun's gravitational field using the first post-Newtonian approximation, advancing the understanding of the solar gravitational lens.

Findings

Derived a solution valid for all distances and angles using confluent hypergeometric functions.

Developed a wave-theoretical description of the solar gravitational lens.

Analyzed the optical properties of the SGL for high-resolution exoplanet imaging.

Abstract

We consider the propagation of electromagnetic (EM) waves in the gravitational field of the Sun within the first post-Newtonian approximation of the general theory of relativity. We solve Maxwell's equations for the EM field propagating on the background of a static mass monopole and find an exact closed form solution for the Debye potentials, which, in turn, yield a solution to the problem of diffraction of EM waves in the gravitational field of the Sun. The solution is given in terms of the confluent hypergeometric function and, as such, it is valid for all distances and angles. Using this solution, we develop a wave-theoretical description of the solar gravitational lens (SGL) and derive expressions for the EM field and energy flux in the immediate vicinity of the focal line of the SGL. Aiming at the potential practical applications of the SGL, we study its optical properties and discuss its suitability for direct high-resolution imaging of a distant exoplanet.

Equations

$$d{s}^2$$

$$u=1+c^{-2}U+\mathcal{O}(c^{-4}),\text{where}U(\mathrm{x})=G\int \frac{\rho (x^{\prime })d^3{x}^{\prime }}{|\mathrm{x}-{\mathrm{x}}^{\prime }|}.$$

$$\displaystyle\partial_{l}F_{ik}+\partial_{i}F_{kl}+\partial_{k}F_{li}=0,\qquad\frac{1}{\sqrt{-g}}\partial_{k}\Big{(}\sqrt{-g}F^{ik}\Big{)}=-\frac{4\pi}{c}j^{i},$$

$$\mathrm{D}=\frac{1}{\sqrt{g_{00}}}\mathrm{E}=u\mathrm{E},\mathrm{B}=\frac{1}{\sqrt{g_{00}}}\mathrm{H}=u\mathrm{H}.$$

$$\mathrm{cur}{\mathrm{l}}_{\kappa }\mathbf{E}$$

$$\mathrm{cur}{\mathrm{l}}_{\kappa }\mathbf{H}$$

$$\mathrm{curl}\mathrm{D}$$

$$\mathrm{curl}\mathrm{B}$$

$$\Delta \mathrm{D}$$

$$\Delta \mathrm{B}$$

$$\mathrm{D}=\psi \mathrm{d}{e}^{-i\omega t}\mathrm{and}\mathrm{B}=\psi \mathrm{b}{e}^{-i\omega t},$$

$$\displaystyle\Big{\{}\Delta\psi+k^{2}u^{4}\psi\Big{\}}{\boldsymbol{\mathrm{d}}}$$

$$\displaystyle\Big{|}\frac{|{\boldsymbol{\mathrm{\nabla}}}^{2}{\boldsymbol{\mathrm{d}}}|}{|{\boldsymbol{\mathrm{d}}}|}\Big{|}^{1/2},\quad\frac{|{\boldsymbol{\mathrm{\nabla}}}{\boldsymbol{\mathrm{d}}}|}{|{\boldsymbol{\mathrm{d}}}|},\quad{\rm and}\quad\frac{1}{r}\ll|{\boldsymbol{\mathrm{\nabla}}}\ln\psi|,$$

$$|\nabla \mathrm{D}|\sim |k\mathrm{D}|,|{\nabla }^2\mathrm{D}|\sim |k^2\mathrm{D}|.$$

$$\Delta \psi +k^2{u}^4\psi$$

$$(\nabla \psi \cdot \nabla )\mathrm{d}$$

$$u(\mathrm{r})$$

$$\Delta \psi +k^2(1+\frac{2r_g}{r})\psi$$

$$(\nabla \psi \cdot \nabla )\mathrm{d}$$

$$\displaystyle\psi(\boldsymbol{\mathrm{r}})=\psi_{0}e^{ikz}{}_{1}F_{1}\big{(}ikr_{g},1,ik(r-z)\big{)},$$

$$\psi (\mathrm{r})$$

$${\psi }_0=E_0{e}^{\frac{\pi }{2}k{r}_g}\Gamma (1-ik{r}_g).$$

$${\psi }_{\mathtt{inc}}(\mathrm{r})$$

$${\psi }_{\mathtt{s}}(\mathrm{r})$$

$$\nabla \psi =i\psi \nabla \phi +\mathcal{O}(r_g^2)=i\psi K^0\kappa +\mathcal{O}(r_g^2),$$

$$(\nabla \phi \nabla )\mathrm{d}=K^0(\kappa \cdot \nabla )\mathrm{d}=(\frac{d\mathrm{r}}{d\lambda }\cdot \nabla )\mathrm{d}=\frac{d\mathrm{d}}{d\lambda }=\frac{d{x}^0}{d\lambda }\frac{d\mathrm{d}}{d{x}^0}=K^0\frac{d\mathrm{d}}{d{x}^0}=K^0\frac{d\mathrm{d}}{d\ell }.$$

$$\frac{d\mathrm{d}}{d\ell }$$

$$\mathrm{r}=(\mathrm{r}\cdot \mathrm{k})\mathrm{k}+[\mathrm{k}\times [\mathrm{r}\times \mathrm{k}]]+{\mathrm{r}}_{\mathtt{G}}+\mathcal{O}(r_g^2)=\mathrm{k}\ell +{\mathrm{b}}_0+{\mathrm{r}}_{\mathtt{G}}+\mathcal{O}(r_g^2),$$

$$\mathrm{d}=(\mathrm{d}\cdot \mathrm{k})\mathrm{k}+[\mathrm{k}\times [\mathrm{d}\times \mathrm{k}]]+{\mathrm{d}}_{\mathtt{G}}+\mathcal{O}(r_g^2)=d_{||0}\mathrm{k}+{\mathrm{d}}_{\perp 0}+{\mathrm{d}}_{\mathtt{G}}+\mathcal{O}(r_g^2),$$

$$(\mathrm{d}\cdot \mathrm{r})$$

$$\frac{d{\mathrm{d}}_{\mathtt{G}}}{d\ell }$$

$$\mathrm{d}$$

Full text

I Introduction

According to Einstein’s general theory of relativity Einstein (1915, 1916), gravitation induces refractive properties on spacetime Fock (1959), with massive objects acting as lenses by bending photon trajectories Turyshev (2008) and amplifying brightness of faint sources. Experimental confirmation of the general relativistic gravitational bending of light nearly a century ago Eddington (1919); Dyson et al. (1920) unambiguously established that celestial bodies act as gravitational lenses, deflecting light from distant sources. The properties of gravitational lenses, including light amplification and the appearance of ringlike images (Einstein rings), are well established Chwolson (1924); Einstein (1936) and have a rich literature Liebes (1964); Refsdal (1964); Refsdal and Surdej (1994); Narayan and Bartelmann (1996); Schneider et al. (1992); Wambsganss (1998); Gould (2001); Eshleman (1979); Kraus (1986). Compact, opaque and spherical bodies acting as gravitational lenses could be used as diffractive telescopes to form images of distant objects at extreme resolution Labeyrie (1994).

Unlike an optical lens, a gravitational lens is astigmatic, with the bending angle inversely proportional to the impact parameter of a light ray with respect to the lens. Therefore, such a lens has no single focal point but a focal line. Although all the bodies in the solar system may act as gravitational lenses Turyshev (2008), only the Sun is massive and compact enough for the focus of its gravitational deflection to be within the range of a realistic deep space mission. Its focal line begins at $\sim$547.8 astronomical units (AU). A probe positioned beyond this distance from the Sun could use the solar gravitational lens (SGL) to magnify light from distant objects on the opposite side of the Sun Eshleman (1979); Maccone (2009).

In recent years, the unique properties of the SGL garnered increasing attention. On the one hand, the discovery of numerous exoplanets by the Kepler telescope, including those that may be Earth-like Torres et al. (2015), created interest in methods to image these distant worlds. On the other hand, the success of the Voyager-1 spacecraft, operating at a distance of nearly 140 AU from the Sun, demonstrates the feasibility of long-duration deep-space missions to the outer solar system, including regions where images are formed by the SGL. The idea of using the SGL for direct megapixel high-resolution imaging of an object of extreme interest, such as a habitable exoplanet, was only recently suggested Shao (2013). It was extensively discussed within the context of a recent study at the Keck Institute for Space Studies Stone et al. (2015). In the past, only the amplification properties of the SGL under a set of idealized physical conditions were explored, considering only the gain of a combined receiver consisting of a large parabolic radio antenna, at the focus of which there was a single pixel detector situated on the focal line of the SGL Eshleman (1979); Maccone (2009, 2011). The SGL’s imaging properties, where the image occupies many pixels in the immediate vicinity of the focal line, are still not fully explored (except perhaps for some introductory considerations on geometric raytracing Koechlin et al. (2005, 2006)), especially in a deep-space mission context. In addition, the SGL’s potential for high-resolution spectroscopy should also be considered.

The reason for the large amplification of the SGL is the fact that, as a typical gravitational lens, the SGL forms a folded caustic Gaudi and Petters (2002a, b) in its focal area. As the wavelength of light is much smaller than the Schwarzschild radius of the Sun, the wavefront in the focal region of the SGL is dominated by the caustic and singularities typical for geometric optics. In reality, the geometric singularities are softened and decorated on fine scales by wave effects Berry and Upstill (1982); Berry (1992). Despite leading to divergent results, geometric optics may be used to predict the focal line, and make qualitative arguments about the magnification and the size of the image. However, designing a telescope entails addressing practical questions concerning the magnification, resolution, field of view (FOV), and the plate scale of the imaging system. These parameters are usually estimated by a wave optics approach and are needed to assess the imaging potential of the SGL. Recently, we reported on a method Turyshev (2017) of providing a wave-theoretical description of the SGL, demonstrating that with its light amplification power of $\sim 10^{11}$ (for $\lambda=1~{}\mu$m) and angular resolution of $\lesssim 10^{-10}$ arcsec, the SGL may be used for direct megapixel imaging of an exoplanet. In this paper we provide details of this derivation.

This paper is structured as follows: In Sec. II, we consider propagation of electromagnetic waves in the Einstein’s general theory of relativity (GR). We establish a set of equations that guide the evolution of EM wave in the presence of a static gravitational monopole. We solve these equations in the post-Newtonian approximation of the GR. In Sec. III we find exact solutions for the Debye potentials for the EM waves traversing the field of a static gravitational monopole. We derive the components of the entire EM field and determine the components of the relevant Poynting vector. Our results yield a wave-optical description of a monopole gravitational lens and are valid for any distances and angles, including those in the immediate vicinity of the focal line. In Sec. IV we provide preliminary considerations for imaging with the Solar Gravitational Telescope (SGT) and its potential application for direct multipixel imaging and spectroscopy of an exoplanet. In Section V we discuss our results and the potential of using the SGL for remote investigations of faint distant objects. In an attempt to streamline the discussion, we placed some important but technically lengthy derivations into appendices. Appendix A contains a summary of results concerning the $(3+1)$ decomposition of a general Riemannian metric and relevant useful relations. Appendix B is devoted to a description of light propagation in a weak, static gravitational field. We solve the geodesic equation and model the phase evolution in the context of geometric optics. We also discuss spherical waves in the post-Newtonian gravity. In Appendix C we present useful properties of the confluent hypergeometric functions. Appendix D discusses Coulomb functions. Appendix E introduces Debye potentials as a means to represent the electromagnetic field. Finally, Appendix F discusses the Wentzel–Kramers–Brillouin (WKB) approximation.

II Electromagnetic waves in a static gravitational field

To describe the optical properties of the solar gravitational lens (SGL), we use a static harmonic metric111The notational conventions used in this paper are the same as in Landau and Lifshitz (1988); Turyshev and Toth (2015): Latin indices ($i,j,k,...$) are spacetime indices that run from 0 to 3. Greek indices $\alpha,\beta,...$ are spatial indices that run from 1 to 3. In case of repeated indices in products, the Einstein summation rule applies: e.g., $a_{m}b^{m}=\sum_{m=0}^{3}a_{m}b^{m}$. Bold letters denote spatial (three-dimensional) vectors: e.g., ${\boldsymbol{\mathrm{a}}}=(a_{1},a_{2},a_{3}),{\boldsymbol{\mathrm{b}}}=(b_{1},b_{2},b_{3})$. The dot ($\cdot$) and cross ($\times$) are used to indicate the Euclidean inner product and cross product of spatial vectors; following the convention of Fock (1959), these are enclosed in round and square brackets, respectively. Latin indices are raised and lowered using the metric $g_{mn}$. The Minkowski (flat) spacetime metric is given by $\gamma_{mn}={\rm diag}(1,-1,-1,-1)$, so that $\gamma_{\mu\nu}a^{\mu}b^{\nu}=-({\boldsymbol{\mathrm{a}}}\cdot{\boldsymbol{\mathrm{b}}})$. We use powers of the inverse of the speed of light, $c^{-1}$, and the gravitational constant, $G$ as bookkeeping devices for order terms: in the low-velocity ($v\ll c$), weak-field ($r_{g}/r=2GM/rc^{2}\ll 1$) approximation, a quantity of ${\cal O}(c^{-2})\simeq{\cal O}(G)$, for instance, has a magnitude comparable to $v^{2}/c^{2}$ or $GM/c^{2}r$. The notation ${\cal O}(a^{k},b^{\ell})$ 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 $\ell$. Other notations are explained in the paper. in the first post-Newtonian approximation of the general theory of relativity. The line element for this metric may be given, in spherical coordinates $(r,\theta,\phi)$, as Fock (1959); Turyshev and Toth (2015):

[TABLE]

where, to the accuracy sufficient to describe light propagation in the solar system, the quantity $u$ can be given in terms of the Newtonian potential $U$ as

[TABLE]

The metric (1)–(2) allows us to consider the largest effects of the gravitational field of the Sun on propagation of light, those due to the static distribution of matter inside the Sun. One may also want to consider including solar rotation, but its effect, although measurable, is much less than those of the solar monopole and quadrupole Will (1993). Thus, the solar spin is not present in the metric above. Nevertheless, if needed, one can always consider the effect of the solar rotation on the properties of the SGL using the same methods that are developed in this paper. Also, the gravitational field of the Sun is weak: its potential is $GM/c^{2}r\lesssim 2\times 10^{-6}$ everywhere in the solar system. This allows us to carry out calculations to the first post-Newtonian order, while dropping higher-order terms.

The generally covariant form of Maxwell’s equations for the electromagnetic (EM) field is well known:

[TABLE]

where $F_{ik}$ is the antisymmetric Maxwell tensor of the EM field Landau and Lifshitz (1988), $g_{mn}$ is a Riemann metric tensor with $g=\det g_{mn}$ its determinant, and $\partial_{k}$ are coordinate derivatives.

Note that in this paper we study the propagation of the EM waves on the background of the Sun without accounting for the corona. That is to say, we do not consider contributions of the solar plasma to light propagation. The refractive properties of the solar corona are such that for high-frequency EM waves such as visible light, one may neglect the refractive effects of the solar plasma Turyshev and Andersson (2003). This may not be the case for any noise contribution to an image due to the brightness of the corona. These issues will be addressed elsewhere. Here we consider only a purely gravitational case, accounting only for the shadow due to a spherical Sun, but ignoring the corona.

II.1 Maxwell’s equations in three-dimensional form

To study the problem of gravitational lensing, we need to present equations (3) in a three-dimensional form. To this effect, we consider a $(3+1)$ decomposition of a generic metric $g_{mn}$ (e.g., using methods discussed in §84 of Landau and Lifshitz (1988)). We introduce quantities describing physical vectors of the EM field, namely the 3-vectors ${\boldsymbol{\mathrm{E}}},{\boldsymbol{\mathrm{D}}}$ and antisymmetric 3-tensors $B_{\alpha\beta}$ and $H_{\alpha\beta}$: $E_{\alpha}=F_{0\alpha},\,D^{\alpha}=-\sqrt{g_{00}}F^{0\alpha},\,B_{\alpha\beta}=F_{\alpha\beta},\,H^{\alpha\beta}=\sqrt{g_{00}}F^{\alpha\beta}$ (see the Problem in §90 of Landau and Lifshitz (1988)). These quantities are not independent. In the case of a static metric, such as that given by (1), for which $g_{0\alpha}=0$ and $\partial_{0}g_{mn}=0$, they are related by the following identities:

[TABLE]

Given the definitions above, Eqs. (3) can be written in the following three-dimensional form:

[TABLE]

where the differential operators ${\rm curl}_{\kappa}{\boldsymbol{\mathrm{F}}}$ and ${\rm div}_{\kappa}{\boldsymbol{\mathrm{F}}}$, for the static metric (1), are taken with respect to the three-dimensional metric tensor $\kappa_{\alpha\beta}=-g_{\alpha\beta}$ (see (146)–(147) and (154)–(155) in Appendix A for details).

We consider the propagation of an EM wave in the vacuum where no sources or currents exist, i.e., $j^{k}=(\rho,{\boldsymbol{\mathrm{j}}})=0$. For the metric (1), using the definitions (4) together with (147) and (154)–(155), we obtain the following form for Maxwell’s equations (5)–(6):

[TABLE]

where the differential operators ${\rm curl}\,{\boldsymbol{\mathrm{F}}}$ and ${\rm div}\,{\boldsymbol{\mathrm{F}}}$ are now with respect to the usual 3-space Euclidean flat metric.

Using the standard identities of vector calculus involving the $\nabla$ operator Korn and Korn (1968); Callen (2006) and a bit of algebra, one can verify that ${\boldsymbol{\mathrm{D}}}$ and ${\boldsymbol{\mathrm{B}}}$ obey the following wave equations:

[TABLE]

All the properties of a propagating EM wave in the presence of a weak and static post-Newtonian gravitational field are encoded in (9)–(10). Note that the last two terms in (9) and (10) are important for establishing the directional and polarization properties of EM field represented by the vectors ${\boldsymbol{\mathrm{D}}}$ and ${\boldsymbol{\mathrm{B}}}$.222Equations (9) and (10) are rather well known. In fact, they are similar to (5)–(6) in Chapter 1.2 of Born and Wolf (October 13, 1999), written for an EM wave propagating in a refractive medium. A form of Maxwell’s equations, similar to (9)–(10), appears any time when one deals with EM waves propagating in a medium with a variable index of refraction, such as in the case of optical waveguides Bures (2009); Chew (2015). This form emphasizes the fact that a weak gravitational field also induces effective refractive properties on spacetime Fock (1959). These properties may be investigated using the tools of classical optics Landau and Lifshitz (1988); Born and Wolf (October 13, 1999). As we show in this paper, omitting these terms (e.g., as in Deguchi and Watson (1986)) may lead to the loss of important information about the propagation direction and the amplitude of the EM field. These equations can be used to study propagation of EM waves in the presence of a weak and static gravitational field. In particular, in the case of solving the problem of diffraction of the EM waves, they can be used to describe both incident and scattered waves. This is the knowledge that helps us study the properties of the EM field in the image plane when dealing with the imaging properties of the SGL.

II.2 Solving Maxwell’s equations

We look for a solution to the wave equations (9)–(10) for the fields ${\boldsymbol{\mathrm{D}}}$ and ${\boldsymbol{\mathrm{B}}}$ in the following generic form:

[TABLE]

where $\psi(\boldsymbol{\mathrm{r}})$ is a scalar function representing the intensity of a monochromatic EM wave along the path of its propagation, ${\boldsymbol{\mathrm{d}}}(\boldsymbol{\mathrm{r}})$ and ${\boldsymbol{\mathrm{b}}}(\boldsymbol{\mathrm{r}})$ are unit vectors specifying the direction of the wave’s propagation and its polarization, and $\omega$ is the frequency of the wave. Although (11) gives the two fields as complex quantities, the actual physical fields ${\boldsymbol{\mathrm{D}}}$ and ${\boldsymbol{\mathrm{B}}}$ are given by the real part of these expressions. Then, for example, the wave equation (9) can be presented in terms of equations for the new quantities $\psi$ and $\boldsymbol{\mathrm{d}}$ as

[TABLE]

where $k=\omega/c$ is the wave number, as usual. As we intend to work with optical frequencies, Eq. (12) may be simplified. For high-frequency propagation, the representation of ${\boldsymbol{\mathrm{D}}}$ given in (11) implies Mo and Papas (1971) that

[TABLE]

which means that $L_{d}\gg L_{\psi}$ and $r\gg L_{\psi}$, where $L_{d}$ and $L_{\psi}$ represent the typical length scales over which the changes in ${\boldsymbol{\mathrm{d}}}$ and $\psi$, respectively, are significant Mo and Papas (1971) (same applies to ${\boldsymbol{\mathrm{b}}}$). In other words, we can see that ${\boldsymbol{\mathrm{d}}}$ (and ${\boldsymbol{\mathrm{b}}}$) vary slowly, but $\psi$ varies rapidly when $k\rightarrow\infty$, resulting in the following relationships:

[TABLE]

Thus, in the case of high-frequency EM wave propagation, the following two equations hold simultaneously333Note that representations similar to (15)–(16) occur when raytracing methods are used to describe the propagation of high-frequency EM waves in optical waveguides Chew (2015). The numerical tools developed in that area may be quite useful to model imaging with the SGL.:

[TABLE]

Below, we focus our discussion on the largest contribution to the gravitational deflection of light, namely that produced by the field of a gravitational monopole. In this case, the Newtonian potential in (2) is given as $U({\boldsymbol{\mathrm{r}}})={r_{g}}/{2r}+{\cal O}(r^{-3},c^{-4}),$ where $r_{g}=2GM/c^{2}$ is the Schwarzschild radius of the source444If needed, our approach, in conjunction with the tools developed in Kopeikin et al. (2011); Turyshev et al. (2014), may be used to account for the contributions from higher order gravitational multipole moments. For details, see Appendices B.1 and B.2.. Therefore, the quantity $u$ in (1) and its logarithmic gradient ${\boldsymbol{\mathrm{\nabla}}}\ln u^{2}$ have the form

[TABLE]

As a result, the system of Eqs. (15)–(16) takes the form

[TABLE]

Experiments in the presence of weak gravitational fields, such as those present in our own solar system Turyshev (2008), are often described using geodesic equations. These equations determine the direction of light propagation and related relativistic frequency shifts Turyshev et al. (2013, 2014). However, geodesic equations provide no information about gravitationally induced changes in the intensity of light. In the solar system, such changes are quite small and very difficult to detect. This is precisely the focus of our interest when we consider the solar gravitational telescope scenario.

To investigate the intensity changes that result from the gravitational amplification of light, we need to develop a wave-theoretical treatment of light propagation in gravity. Equations (18)–(19) could be used for this purpose. These are derived from the wave equations (9)–(10) and provide a complete description of an EM wave propagating in a weak and static gravitational field (which, according to Fock Fock (1959), acts as a variable index of refraction). Specifically, (18) determines the change in the intensity of the EM field, while (19) describes the changes in the direction of propagation of the EM wave and its polarization.

We can solve Eqs. (18)–(19) iteratively to first order in $G$. This can be done along the path of wave propagation, which is established by relying on the geodesic equation (see Appendix B).

II.3 Solving the wave equations

To find the formal solution for the EM field, we begin with (18). This equation is well known: it is nearly identical to the time-independent Schrödinger equation that describes the scattering problem in a Coulomb potential in nuclear physics555A choice of constants $r_{g}=-\gamma/k$ makes (18) identical to the time-independent Schrödinger equation describing the scattering problem in a Coulomb potential Messiah (1968), where $\gamma=Z_{1}Z_{2}e^{2}/\hbar v$ with $Z_{1}e$, $Z_{2}e$ are the charges of the two particles and $v$ is their relative velocity. The first analytical solution to (18) was given by Mott in 1928 Mott (1928). A more elegant form was found a few months later by Gordon Gordon (1928), using the ansatz $\psi({\boldsymbol{\mathrm{r}}})=e^{ikz}f(r-z)$. The complex-valued function $f$ describes the perturbation of the incoming plane wave and transforms (18) into a solvable differential equation for $f$ Messiah (1968). The same equation appears in other problems of modern physics: for instance, in problems describing photothermal single-particle Rutherford scattering microscopy that involves the scattering of waves by a $1/r$ refractive index profile formed by the presence of a point-like heat source in a homogeneous medium (e.g., Selmke and Cichos (2013a, b); Selmke (2015)).. We take a spherical coordinate system $(r,\theta,\phi)$ and also use Cartesian coordinates such that the $x$ and $z$ axes, $r$ and angle $\theta$ are given as in Fig. 2.

We consider the propagation of a monochromatic EM wave along the $z$-axis coming from a source at infinity. As is known from textbooks on quantum mechanics (e.g., Morse and Feshbach (1953); Messiah (1968); Landau and Lifshitz (1989); Schiff (1968); Burke (2011)), (18) has a solution that is regular at the origin, which can be given as

[TABLE]

where $z$ is the projection of $\boldsymbol{\mathrm{r}}$ onto the optical axis (i.e., a coordinate along that axis; see Fig. 2), $\psi_{0}$ is an integration constant and ${}_{1}F_{1}$ is the confluent hypergeometric function Abramowitz and Stegaun (1965) (also known as Kummer’s function of the first kind, $M[\alpha,\beta,w]$; see Appendix C for more details and useful relations).

The solution (20), also shown in Fig. 2, describes a wave coming from a large distance along the $z$ axis (for the relevant geometry, see Fig. 2) and generalizes the incoming plane wave solution $\psi_{0}(\boldsymbol{\mathrm{r}})=e^{ikz}$, which is familiar from studying wave propagation in Euclidean spacetime. In fact, Eq. (20) reduces to $e^{ikz}$ when $r_{g}\rightarrow 0$. Thus, one may use the solution (20) to describe the incident “plane wave” that is sourced at infinity, in the presence of a gravitational monopole with a $1/r$ potential. All the important contributions to $\psi_{0}(\boldsymbol{\mathrm{r}})$ from gravitation are contained in the function ${}_{1}F_{1}[\alpha,\beta,w]$.

Given the asymptotic properties of ${}_{1}F_{1}[\alpha,\beta,w]$ from (233) (see details of derivation in Appendix C.2), we obtain the asymptotic form of Eq. (20) as

[TABLE]

This approximation is valid for large values of the argument $k(r-z)\gg 1$ and for angles $\theta$ satisfying $\theta\gtrsim\sqrt{2r_{g}/r}$ (see Fig. 3). This region is relatively far from the optical axis; light refraction here is well described by geometric optics. This solution offers a good starting point for the development of the wave-theoretical treatment of the SGL.

Typically, one normalizes the solution at large distances from the deflecting center by requiring that the function $\psi$ behaves as $\lim_{k(r-z)\to\infty}\psi\psi^{*}=1$ (a.k.a. Gamow normalization Messiah (1968); Schiff (1968)), which results in $\psi_{0}=e^{\frac{\pi}{2}kr_{g}}{\Gamma(1-ikr_{g})}$. However, in our case, we require that at larger distances from the deflector the intensity of the EM field, $\psi$, is to be equal to that at the source, namely $\lim_{k(r-z)\to\infty}\psi\psi^{*}=E_{0}^{2}$ (in the vacuum $E_{0}=H_{0}$). This results in the following choice for the constant $\psi_{0}$:

[TABLE]

As a result, at large distances from the deflector, the incident wave ($\psi_{\tt inc}(\boldsymbol{\mathrm{r}})$, given by the first term in (21)) and the scattered wave ($\psi_{\tt s}(\boldsymbol{\mathrm{r}})$, given by the second term in (21)) take the following asymptotic forms:

[TABLE]

The solution provided in the form of Eqs. (23)–(24) is well known from the Coulomb scattering problem in nuclear physics. What is its meaning in general relativity? Equations (23)–(24) do not exhibit the familiar geodesic behavior that is characteristic of rays of light. Nonetheless, with some algebra (see Appendix B.2, Eq. (189)), we can show that (23) is consistent with a solution for the phase of an EM wave propagating in the background of a weak and static gravitational field. For a wave moving from a remote source along the $z$-axis, $({\boldsymbol{\mathrm{k}}}\cdot{\boldsymbol{\mathrm{r}}})=z$, where $\boldsymbol{\mathrm{k}}$ is the unperturbed unit vector of the photon’s trajectory (see Sec. B.1 for details). Therefore, from (189) and (170), for a wave moving along a geodesic, we obtain, for the change of phase along the path, \delta\varphi=k\big{(}({\boldsymbol{\mathrm{k}}}\cdot({\boldsymbol{\mathrm{r}}}-{\boldsymbol{\mathrm{r}}}_{0}))+r_{g}\ln({r+({\boldsymbol{\mathrm{k}}}\cdot{\boldsymbol{\mathrm{r}}})})/({r_{0}+({\boldsymbol{\mathrm{k}}}\cdot{\boldsymbol{\mathrm{r}}}_{0})})+{\cal O}(r_{g}^{2})\big{)}=k\big{(}z-z_{0}-r_{g}\ln({r-z})/({r_{0}-z_{0}})+{\cal O}(r_{g}^{2})\big{)}. Thus, the time-independent part of the phase of the incident wave has the form \varphi(r)=k\big{(}z-r_{g}\ln k(r-z)+{\cal O}(r_{g}^{2})\big{)}, given by (23), which is consistent with a geodesic solution.

To understand the meaning of Eq. (24), we rewrite it using $z=r\cos\theta$ as follows:

[TABLE]

with $f(\theta)$ being the scattering amplitude familiar from nuclear scattering.

One can see that the phase in the first expression in (25) is consistent with the phase of a radial geodesic or that of an outgoing spherical wave (see discussions in Appendix B.4, Eq. (211)). From (25), for the change of phase along a radial geodesic, we have \delta\varphi=k_{0}\big{(}r-r_{0}+r_{g}\ln{r}/{r_{0}}+{\cal O}(r_{g}^{2})\big{)}, which indicates that the time-independent part of the phase of a scattered wave is that of a spherical wave given by (211) as \varphi(r)=k(r+r_{g}\ln 2kr+{\cal O}(r_{g}^{2})\big{)} and is consistent with the phase of a radial geodesic (190). The quantity $f(\theta)$ in (25) is the scattering amplitude that was first derived by Rutherford for the electron scattering problem in nuclear physics Matzner (1968) and has been confirmed in many experiments. This amplitude modifies the outgoing spherical wave (211) (discussed in Appendix B.4).

Therefore, the two solutions to the time-independent wave equation (18) are both consistent with the familiar geodesic solutions in a weak and static gravitational field (as discussed in Secs. B.1 and B.2). The phase of the incident wave is consistent with the geodesic solution (189), while the scattered wave is consistent with a spherical wave solution (211) or, equivalently, with radial geodesics (190). With this knowledge we may already identify these features in (20). Solution to this equation is given in Fig. 2, which clearly shows the presence of both of these waves, namely the Coulomb-modified incident wave and the outgoing spherical wave modified by the scattering amplitude.

To interpret the approximate solutions (23)–(24), it helps to study the schematic geometry shown in Fig. 2. Light from a distant source reaches the point of observation (black dot on the right-hand side) via two paths. When the point of observation is a significant distance away from the focal line, these two paths are qualitatively different.

Perturbations to the path on the same side of the focal line as the point of observation (the “top” ray of light in Fig. 2) are dominated by deflection. Neighboring rays in a tight family of rays diverge (“spread out”) minimally. Therefore, this path is well approximated by Eq. (23), which describes a plane wave slightly perturbed by deflection. This wave is shown in the top left panel of Fig. 3. We call this part of the solution the (perturbed) “incident” wave.

In contrast, perturbations to the path on the side of the focal line opposite to the point of observation (the “bottom” ray of light in Fig. 2) are dominated by scattering. These rays reach the point of observation because they have a small impact parameter and a large angle of deflection. As a result, even neighboring rays, with only slightly different impact parameters, will suffer noticeably different deflections. The resulting wavefront is dominated by this divergent (“spreading out”) behavior, and thus it is well approximated as a spherical wave emanating from the gravitational lens itself; that is to say, Eq. (24), which is depicted in the top right panel of Fig. 3. We refer to this part of the solution as the “scattered” wave.

The combination of the deflected plane (incident) wave and the perturbed spherical (scattered) wave, shown in the bottom right panel of Fig. 3, offers a good approximation of the propagating wavefront everywhere except for the vicinity of the focal line. The bottom left panel of Fig. 3 compares this approximation to the original form of (20) (see Fig. 2). The conical region on the right-hand side of this figure is where the geometric optics approximation fails Turyshev and Andersson (2003). It is inconvenient, since the total solution for the wave function (20) gives the correct asymptotic expression for any angle. Technically, this is because there are no known approximations of the confluent hypergeometric function ${}_{1}F_{1}$ that are simultaneously valid both for large distances and also for small angles Deguchi and Watson (1986).

For a monopole lens, the very fact that any two rays intersect at the focal line means that these rays essentially have identical impact parameters. In this case, an observer will see a thin annulus around the lens representing the Einstein ring formed by the amplified intensity of the incident light coming from the direction of the source. At any given point outside the focal line, the rays will have different impact parameters. An observer will see two images of unequal brightness of a distant point source, one on each side of the lens. Far enough from the focal line, one ray will suffer minimal deflection due to its large impact parameter. Meanwhile, the other ray will not only be deflected but also dispersed as neighboring rays diverge. The approximation given by (23) describes the ray with minimal deflection, i.e., a slightly perturbed version of the incident wave. A weaker contribution dominated by the factor $r_{g}/(r-z)$ is given by (24), which approximates these diverging rays with a small impact parameter (passing close to the lens) as a perturbed spherical wave originating from the lens.

At a sufficient distance from the focal line, the impact parameter needed for one of the rays to reach these points will be smaller than the physical radius of the Sun. Therefore, these incident rays will be blocked by the Sun and no scattered rays will be produced. In these cases, an observer will see only one image described by (23).

II.4 Amplitude evolution of the incident wave

Given the solution (23) for the incident wave, we can now proceed with solving (19). This helps us determine the polarization changes of the EM wave. First, by defining $\varphi$ to be the phase of the incident wave $\psi_{i}$ in (23) and using the usual definition for the wave number, $K^{m}=dx^{m}/d\lambda=g^{mn}\partial_{n}\varphi$ or $K_{m}=\partial_{m}\varphi$, we have

[TABLE]

where ${\boldsymbol{\mathrm{\kappa}}}$ is the unit vector along the direction of the wave vector, such that $K^{\alpha}=K^{0}\kappa^{\alpha}$. Note that to ${\cal O}(r_{g}^{2})$, ${\boldsymbol{\mathrm{\kappa}}}={\boldsymbol{\mathrm{K}}}/|{\boldsymbol{\mathrm{K}}}|$ has the form ${\boldsymbol{\mathrm{\kappa}}}={\boldsymbol{\mathrm{k}}}+{\boldsymbol{\mathrm{\kappa}}}_{\tt G}+{\cal O}(r_{g}^{2})$, with $\boldsymbol{\mathrm{k}}$ being the unperturbed part and ${\boldsymbol{\mathrm{\kappa}}}_{\tt G}$ being the post-Newtonian term, with both of them given explicitly by (158).

It is convenient to introduce a parameter $\ell$, which is defined along the path of the photon’s trajectory as $\ell=({\boldsymbol{\mathrm{k}}}\cdot{\boldsymbol{\mathrm{r}}})=({\boldsymbol{\mathrm{k}}}\cdot{\boldsymbol{\mathrm{r}}}_{0})+c(t-t_{0})$ (see (168) and discussion in Appendix B.1.) Given $K^{0}=dx^{0}/d\lambda$, we have $d\ell=K^{0}d\lambda$, and, thus:

[TABLE]

Substituting (26) and (27) in (19), we obtain the following equation that can be used to study the post-Newtonian evolution of $\boldsymbol{\mathrm{d}}$:

[TABLE]

Given the two linearly independent unit 3-vectors ${\boldsymbol{\mathrm{n}}}={\boldsymbol{\mathrm{r}}}/r$ and ${\boldsymbol{\mathrm{\kappa}}}$, we can define a triplet of unit vectors, ${\boldsymbol{\mathrm{\kappa}}}$, ${\boldsymbol{\mathrm{\pi}}}=[{\boldsymbol{\mathrm{\kappa}}}\times{\boldsymbol{\mathrm{n}}}]/|[{\boldsymbol{\mathrm{\kappa}}}\times{\boldsymbol{\mathrm{n}}}]|$, and ${\boldsymbol{\mathrm{\epsilon}}}=[{\boldsymbol{\mathrm{\pi}}}\times{\boldsymbol{\mathrm{\kappa}}}]$, forming a local right-handed orthonormal basis: $({\boldsymbol{\mathrm{\kappa}}}\cdot{\boldsymbol{\mathrm{\pi}}})=({\boldsymbol{\mathrm{\kappa}}}\cdot{\boldsymbol{\mathrm{\epsilon}}})=({\boldsymbol{\mathrm{\pi}}}\cdot{\boldsymbol{\mathrm{\epsilon}}})=0$ (see discussion in Appendix B.3). Then, we can write Callen (2006) the vector ${\boldsymbol{\mathrm{r}}}$ in this basis as

[TABLE]

where ${\boldsymbol{\mathrm{r}}}_{\tt G}\sim{\cal O}(r_{g})$ is the post-Newtonian part of ${\boldsymbol{\mathrm{r}}}$ (derived in (177)) and we used (167) to write $({\boldsymbol{\mathrm{r}}}\cdot{\boldsymbol{\mathrm{k}}})=\ell$ and ${\boldsymbol{\mathrm{b}}}_{0}=[{\boldsymbol{\mathrm{k}}}\times[{\boldsymbol{\mathrm{r}}}\times{\boldsymbol{\mathrm{k}}}]]+{\cal O}(r_{g})$, is the impact parameter (169). Similarly, we can write $\boldsymbol{\mathrm{d}}$ as

[TABLE]

where $d_{||0}=({\boldsymbol{\mathrm{d}}}\cdot{\boldsymbol{\mathrm{k}}})+{\cal O}(r_{g})$ and ${\boldsymbol{\mathrm{d}}}_{\perp 0}=[{\boldsymbol{\mathrm{k}}}\times[{\boldsymbol{\mathrm{d}}}\times{\boldsymbol{\mathrm{k}}}]]+{\cal O}(r_{g})$ are the components of ${\boldsymbol{\mathrm{d}}}$ in the directions parallel and orthogonal to ${\boldsymbol{\mathrm{k}}}$, correspondingly, and ${\boldsymbol{\mathrm{d}}}_{\tt G}$ is the post-Newtonian part of vector ${\boldsymbol{\mathrm{d}}}$. Next, we have:

[TABLE]

As a result, Eq. (28) takes the form

[TABLE]

Taking into account that $d_{||0}$ and ${\boldsymbol{\mathrm{d}}}_{\perp 0}$ are constant, we integrate (32) with respect to $\ell$ from $-\infty$ to $\ell$ and obtain a solution for the components of ${\boldsymbol{\mathrm{d}}}={\boldsymbol{\mathrm{d}}}_{0}+{\boldsymbol{\mathrm{d}}}_{\tt G}+{\cal O}(r_{g}^{2})$ in the local basis. The ${\boldsymbol{\mathrm{B}}}$ field will evolve in a similar manner. As a result, the solutions for ${\boldsymbol{\mathrm{d}}}$ and ${\boldsymbol{\mathrm{b}}}$ are both real and have the following form:

[TABLE]

We introduce a right-handed Cartesian coordinate system $(x,y,z)$ with corresponding unit vectors $({\boldsymbol{\mathrm{e}}}_{x},{\boldsymbol{\mathrm{e}}}_{y},{\boldsymbol{\mathrm{e}}}_{z})$ and with the origin at the center of mass of the Sun. We take the $z$-axis to be directed along the unperturbed direction of the light ray, i.e., along the vector ${\boldsymbol{\mathrm{k}}}$, while the $x$ and $y$ axes will be directed along the unperturbed directions set by the vectors ${\boldsymbol{\mathrm{\epsilon}}}$ and ${\boldsymbol{\mathrm{\pi}}}$, correspondingly. In Appendix B.3 we show that in this coordinate system the vectors ${\boldsymbol{\mathrm{\kappa}}},{\boldsymbol{\mathrm{\pi}}},{\boldsymbol{\mathrm{\epsilon}}}$ take the following form (see (202)–(204)):

[TABLE]

Also, in this coordinate system, $d_{||0}=d_{z0}$ and ${\boldsymbol{\mathrm{d}}}_{\perp 0}=(d_{x0},d_{y0},0)$, ${\boldsymbol{\mathrm{b}}}_{0}=[{\boldsymbol{\mathrm{k}}}\times[{\boldsymbol{\mathrm{r}}}\times{\boldsymbol{\mathrm{k}}}]]+{\cal O}(r_{g})=(x,y,0)+{\cal O}(r_{g})$, and, thus, $({\boldsymbol{\mathrm{d}}}_{\perp 0}\cdot{\boldsymbol{\mathrm{b}}}_{0})=d_{x0}x+d_{y0}y+{\cal O}(r_{g})$. Similarly, we have $({\boldsymbol{\mathrm{b}}}_{\perp 0}\cdot{\boldsymbol{\mathrm{b}}}_{0})=b_{x0}x+b_{y0}y+{\cal O}(r_{g}).$ We choose the components of the incident wave so that it represents a transverse electric and transverse magnetic (TEM) wave, namely we require: $d_{z0}=d_{y0}=b_{z0}=b_{x0}=0$. Based on (33)–(34), the directional vectors of this EM field evolve as

[TABLE]

As a result, substituting (36)–(37) and (23) into (11), accounting for the fact that $\ell=({\boldsymbol{\mathrm{k}}}\cdot{\boldsymbol{\mathrm{r}}})=z$ and using (170) in the second term in (36) and (37), and also taking the amplitudes of the unit vectors of the EM field at the source to be $d_{x0}=b_{y0}=1$, we present the EM field of the incident wave as

[TABLE]

which, with the help of (35), indicates that ${\boldsymbol{\mathrm{D}}}_{\tt inc}\propto{\boldsymbol{\mathrm{\epsilon}}}$ and ${\boldsymbol{\mathrm{B}}}_{\tt inc}\propto{\boldsymbol{\mathrm{\pi}}}$. As the local base vectors ${\boldsymbol{\mathrm{\epsilon}}},{\boldsymbol{\mathrm{\pi}}}$ and ${\boldsymbol{\mathrm{\kappa}}}$ are forming a triplet of orthonormal vectors, the three vectors ${\boldsymbol{\mathrm{D}}}_{\tt inc}$, ${\boldsymbol{\mathrm{B}}}_{\tt inc}$ and $\boldsymbol{\mathrm{\kappa}}$ that characterize the incident wave (38)–(39) are also orthogonal to each other, namely from (35) one can verify that $({\boldsymbol{\mathrm{D}}}_{\tt inc}\cdot{\boldsymbol{\mathrm{B}}}_{\tt inc})=({\boldsymbol{\mathrm{D}}}_{\tt inc}\cdot{\boldsymbol{\mathrm{\kappa}}})=({\boldsymbol{\mathrm{B}}}_{\tt inc}\cdot{\boldsymbol{\mathrm{\kappa}}})=0+{\cal O}(r_{g}^{2})$. So, as expected, the components orthogonal to the wave vector do not change as a photon moves along its trajectory, which, in the case of ${\boldsymbol{\mathrm{D}}}_{\tt inc}$, (38), is in the plane spanned by ${\boldsymbol{\mathrm{e}}}_{x}$ and ${\boldsymbol{\mathrm{e}}}_{z}$. Thus, the gravitational field of a static monopole does not change the polarization of an EM wave. At the same time, the component along the wave vector is mixed with the orthogonal component and rotates by a small angle $\delta\theta=({r_{g}}/{b_{0}})(1+({{\boldsymbol{\mathrm{r}}}\cdot{\boldsymbol{\mathrm{k}}}})/{r})$ as it moves along the trajectory with the entire EM wave being perpendicular to the wave vector (similar results were reported in Mo and Papas (1971)). Similar behavior is evident from (39) for ${\boldsymbol{\mathrm{B}}}_{\tt inc}$ in the plane formed by vectors ${\boldsymbol{\mathrm{e}}}_{y}$ and ${\boldsymbol{\mathrm{e}}}_{z}$.

To proceed with the solution of the scattering problem, we need to transform (38) and (39) from Cartesian into spherical coordinates. The curvilinear coordinates appropriate to represent the problem are the spherical polar coordinates $(r,\theta,\phi)$ defined as usual by the relationships $(x,y,z)=r(\sin\theta\cos\phi,\,\sin\theta\sin\phi,\,\cos\theta)$. Transforming (38) and (39) from the Cartesian system $(x,y,z)$ to this new system of spherical coordinates according to the usual rules of such coordinate transformations Korn and Korn (1968), we obtain the incident wave, ${\boldsymbol{\mathrm{D}}}_{\tt inc}$ and ${\boldsymbol{\mathrm{B}}}_{\tt inc}$, in the following form:

[TABLE]

where $u$ is given by (17) and \psi_{i}({\boldsymbol{\mathrm{r}}})=e^{ik\big{(}r\cos\theta-r_{g}\ln kr(1-\cos\theta)\big{)}} is the incident wave (23). As we can see, the phase and the directional vector of the incident wave are both Coulomb-modified. This reflects on the fact that the long-range $1/r$ field due to the gravitational monopole changes the incident wave even at large distances from the deflector.

II.5 Amplitude evolution of the scattered wave

We now consider the evolution of the amplitude of the scattered wave. Similarly to the incident wave, we may solve (19) for the scattered wave given the solution (24). This helps us determine the polarization changes of the scattered EM wave. First, we recognize the fact that the amplitude in (24) is a slowly varying function of distance while the phase varies rapidly. Thus, we may consider only the phase in finding solution for (19). Following the approach demonstrated in Sec. II.4, we may present (19) for radial geodesics, i.e., ${\boldsymbol{\mathrm{k}}}={\boldsymbol{\mathrm{n}}}$, as

[TABLE]

where the parameter $\ell$ now is $\ell=r=r_{0}+c(t-t_{0})$. Similarly to the discussion of the propagation of the incident wave amplitude, we present $\boldsymbol{\mathrm{d}}$ as

[TABLE]

where, in this case, $d_{||0}=({\boldsymbol{\mathrm{d}}}\cdot{\boldsymbol{\mathrm{n}}})+{\cal O}(r_{g})$ and ${\boldsymbol{\mathrm{d}}}_{\perp 0}=[{\boldsymbol{\mathrm{n}}}\times[{\boldsymbol{\mathrm{d}}}\times{\boldsymbol{\mathrm{n}}}]]+{\cal O}(r_{g})$ are the components of ${\boldsymbol{\mathrm{d}}}$ in the directions parallel and orthogonal to ${\boldsymbol{\mathrm{n}}}$, correspondingly, and ${\boldsymbol{\mathrm{d}}}_{\tt G}$ is the post-Newtonian part of vector ${\boldsymbol{\mathrm{d}}}$. Then, taking into account that $d\ell=cdt$ and, thus, $d(d_{||0}{\boldsymbol{\mathrm{n}}})/d\ell=d{\boldsymbol{\mathrm{d}}}_{\perp 0}/d\ell=0$, we can present Eq. (42) in the following form:

[TABLE]

Taking into account that $d_{||0}$ and ${\boldsymbol{\mathrm{d}}}_{\perp 0}$ are constant, we integrate (44) with respect to $\ell$ from $-\infty$ to $\ell$ and obtain a solution for the components of ${\boldsymbol{\mathrm{d}}}={\boldsymbol{\mathrm{d}}}_{0}+{\boldsymbol{\mathrm{d}}}_{\tt G}+{\cal O}(r_{g}^{2})$ in the local basis along the radial path. The ${\boldsymbol{\mathrm{B}}}$ field will evolve in a similar manner. As a result, the solutions for ${\boldsymbol{\mathrm{d}}}$ and ${\boldsymbol{\mathrm{b}}}$ of the scattered wave have the following form:

[TABLE]

where we remember that for radial motion $\ell=r$. We again choose the TEM wave, thus, $d_{||0}=b_{||0}=0$ and write the solution (45) in the following form:

[TABLE]

As a result, using the entire solution (25) and normalizing $d_{\perp 0}=b_{\perp 0}=1$, the components of the scattered wave, ${\boldsymbol{\mathrm{D}}}_{\tt s}$ and ${\boldsymbol{\mathrm{B}}}_{\tt s}$, in the spherical coordinate system may be given in the following form:

[TABLE]

As expected, the scattered EM wave is proportional to the scattering amplitude $f(\theta)$ and multiplies the outgoing spherical wave as given by (25). Equations (47)–(48) may be presented in the form showing their explicit dependence on all the parameters involved:

[TABLE]

where $\sigma_{0}$ is the quantity known in nuclear physics as the Coulomb phase shift $\sigma_{0}=\arg\Gamma(1-ikr_{g})$. It is defined via the ratio of two gamma function terms in (25): ${\Gamma(1-ikr_{g})}/{\Gamma(1+ikr_{g})}=e^{2i\sigma_{0}}$.

These expressions complete our description of the scattering problem in the geometric optics. In the next section, we use these results to derive the Poynting vector that characterizes energy transmission in this situation.

II.6 Poynting vector in the geometric optics approximation

We may now compute the components of the Poynting vector in geometric approximation using the solutions for the incident and scattered waves. The components of the Poynting vector are computed as ususal Landau and Lifshitz (1988); Logunov (1987):

[TABLE]

where ${\boldsymbol{\mathrm{D}}}={\boldsymbol{\mathrm{D}}}_{\tt inc}+{\boldsymbol{\mathrm{D}}}_{\tt s}$ and ${\boldsymbol{\mathrm{B}}}={\boldsymbol{\mathrm{B}}}_{\tt inc}+{\boldsymbol{\mathrm{B}}}_{\tt s}$ are the total solutions for the EM field that includes incident and scattered waves. Thus, for (51) we have:

[TABLE]

where ${\boldsymbol{\mathrm{S}}}_{\tt inc}=(c/4\pi u)[{\rm Re}({\boldsymbol{\mathrm{D}}}_{\tt inc})\times{\rm Re}({\boldsymbol{\mathrm{B}}}_{\tt inc})]$ is the Poynting vector due to the incident wave, ${\boldsymbol{\mathrm{S}}}_{\tt s}=(c/4\pi u)[{\rm Re}({\boldsymbol{\mathrm{D}}}_{\tt s})\times{\rm Re}({\boldsymbol{\mathrm{B}}}_{\tt s})]$, is that due to the scattered wave, with {\boldsymbol{\mathrm{S}}}_{\times}=(c/4\pi u)\big{(}[{\rm Re}({\boldsymbol{\mathrm{D}}}_{\tt inc})\times{\rm Re}({\boldsymbol{\mathrm{B}}}_{\tt s})]+[{\rm Re}({\boldsymbol{\mathrm{D}}}_{\tt s})\times{\rm Re}({\boldsymbol{\mathrm{B}}}_{\tt inc})]\big{)} being the interferometric or mixed term. Using the expressions for the incident and scattered fields given by (40)–(41) and (49)–(50), correspondingly, we may compute all the terms on the right-hand side of (52). Then, after averaging (52) over time, we get the needed expressions. Thus, for the Poynting vector of the incident wave with (40)–(41) we have

[TABLE]

As expected, the incident wave propagates along the wave vector ${\boldsymbol{\mathrm{\kappa}}}$, which is given by (35). Using expressions (49)–(50), we compute the Poynting vector for the scattered EM wave as

[TABLE]

Note that this term is below our approximation threshold of ${\cal O}(r_{g}^{2})$ and thus it may be omitted. However, it provides information on the largest contribution from the scattered term alone. Note that if, for a particular value of $r$, the angle $\theta$ decreases to the point where the ratio ${r_{g}}/{2r\sin^{2}\frac{\theta}{2}}$ becomes 1, the term (54) is of the same size as (53). If $\theta$ continues to decrease, the interferometric term in (52) also becomes significant. We derive this term next.

Before we derive an expression for ${\boldsymbol{\mathrm{S}}}_{\times}$, it is instructive to represent $\sigma_{0}$ in (49)–(50) in terms of its functional dependence. For this, we need to evaluate the ratio of two gamma functions in (24). To do that, we will use Stirling’s formula that approximates the gamma function for large values of its argument $|\alpha|\rightarrow\infty$ (e.g., Abramowitz and Stegaun (1965)):

[TABLE]

As a result, we have

[TABLE]

Therefore, to a sufficient accuracy, for large values of the argument $|\alpha|=kr_{g}$ (i.e., when considering the propagation of high frequency EM waves), the quantity $\sigma_{0}$ may given as $2\sigma_{0}=-2kr_{g}\ln(kr_{g}/e)-\frac{\pi}{2}$. This allows us to compute the interference term and present it in the following form:

[TABLE]

where the second term in the argument of $\sin()$ comes both from $\sigma_{0}$ and from the argument of the exponent in the expression (25) for the scattering amplitude, $f(\theta)$. One can see that pretty much for every value of $r$ and $\theta$ the Poynting vector of the incident wave $\bar{\boldsymbol{\mathrm{S}}}_{\tt inc}$ (53) dominates the interference term $\bar{\boldsymbol{\mathrm{S}}}_{\times}$ (57). However, when $\theta$ becomes smaller, the interference term starts to grow. If, for a particular $r$, the ratio ${r_{g}}/{2r\sin^{2}\frac{\theta}{2}}$ approaches 1, the magnitude of $\bar{\boldsymbol{\mathrm{S}}}_{\times}$ becomes comparable to that of $\bar{\boldsymbol{\mathrm{S}}}_{\tt inc}$, reaching the value of \bar{\boldsymbol{\mathrm{S}}}_{\times}=(c/4\pi u^{2})E_{0}^{2}\sin 3kr_{g}\big{\{}\hat{i}_{r}-\hat{i}_{\theta}\sqrt{{2r_{g}}/{r}}\big{\}}+{\cal O}(r_{g}^{2}). If $\theta$ continues to decrease, i.e., when $\theta\rightarrow 0$, the terms representing the scattered (54) and interferometric (57) terms continue to grow, and ultimately diverge on the optical axis, where $\theta=0$. This is precisely the area where geometric optics breaks down, necessitating a wave-theoretical treatment. We develop that treatment in Sec. III.

II.7 Boundary conditions in the geometric optics approximation

Lastly, we note that to develop a solution to a diffraction problem, we need to introduce a set of boundary conditions. These conditions are necessary to select specific values for the arbitrary integration constants that are appropriate for a particular problem under consideration. Considering the case of diffraction of the EM wave by the gravitational field of a large star (i.e., an idealized spherical sun with no luminosity and no corona), we need to consider only two of such conditions:

(i) the asymptotic boundary conditions and

(ii) the physical boundary conditions

(as was done, for instance, in Herlt and Stephani (1975, 1976)).

As far as the asymptotic boundary condition is concerned, we already introduced such a condition when we selected the value for the constant $\psi_{0}$ in (20) in the form of (22). This choice was made to satisfy the condition that at large distances from the deflector the incident wave must resemble the Coulomb-modified plane wave with a unit magnitude (i.e., Gamow conditions), but scaled to match the field intensity at the source, namely $\lim_{k(r-z)\to\infty}\psi\psi^{*}=E_{0}^{2}$. This condition led to the solutions for both incident and scattered waves, given by (40)–(41) and (49)–(50), correspondingly.

However, the solutions for the EM waves that we established describe scattering on an object that is characterized only by its Schwarzschild radius, $r_{g}$. This may be sufficient for the problems describing scattering of massless scalar waves by black holes (e.g., Matzner (1968); Futterman et al. (1988); Andersson (1995); Andersson and Jensen (2001)), but is not enough to describe scattering by the Sun, whose physical size is much larger than its Schwarzschild radius, i.e., $R_{\odot}\gg r_{g}$. Therefore, following Herlt and Stephani (1976) we introduce another requirement that our solution must to satisfy: the fully absorbing boundary condition. This condition requires that for rays with impact parameter less than the solar radius, i.e., $b_{0}\leq R_{\odot}$, no wave propagates behind the Sun and no diffracted wave exists. In the geometric optics approximation this condition introduces the shadow behind the Sun, determines its shape, and moves the interference region to heliocentric distances beyond $z_{0}=547.8$ AU (i.e., the point where two gravitationally deflected rays of light that are just grazing the Sun on its opposite sides will intersect.)

Both of these boundary conditions are useful and will take an explicit analytical form in the case of the wave optics treatment of the scattering of an EM wave by the gravitational field of a large star that we discuss next.

III Electromagnetic wave in the field of a static monopole

In the previous section, we obtained all the tools that are required to investigate the EM field in the interference zone of the SGL. Our next goal is to find a solution to the EM field in that region. In this section, we accomplish this objective using the approach developed for classical diffraction theory, by finding the set of equations that determine the EM field via Debye potentials and then matching these equations with the incident wave.

III.1 Representing the field in terms of Debye potentials

It is known Born and Wolf (October 13, 1999); Kerker (1969); van de Hulst (1981) that Maxwell’s equations can be represented in terms of the electric Debye potential ${}^{e}{\hskip-1.0pt}\Pi$ and the magnetic Debye potential ${}^{m}{\hskip-1.0pt}\Pi$. This also applies to the case of an EM wave propagating in the static gravitational field of a Schwarzschild black hole or a large star Herlt and Stephani (1975, 1976); Philipp and Perlick (2015). In Appendix E we demonstrate how such a representation may be done for the EM wave propagating in the vacuum in the background of a weak and static gravitational field, represented by the metric (1)–(2), which is a good approximation for the gravitational field in the solar system. The complete solution for the EM field may be given as (see (305)–(310) for details):

[TABLE]

This solution can be derived from the two potentials ${}^{e}{\hskip-1.0pt}\Pi$ and ${}^{m}{\hskip-1.0pt}\Pi$, which both have to satisfy the same differential equation (300), which is just the wave equation (see (303)):

[TABLE]

where $\Pi$ can be either ${}^{e}{\hskip-1.0pt}\Pi$ or ${}^{m}{\hskip-1.0pt}\Pi$. Typically Born and Wolf (October 13, 1999), in spherical polar coordinates (see Fig. 2 for details), the solution of this equation is represented using an expansion, with terms in the form

[TABLE]

and with coefficients that are determined by boundary conditions. Direct substitution into (299) reveals that the functions $R,\Theta$ and $\Phi$ must satisfy the following ordinary differential equations:

[TABLE]

The solution to (65) is given as usual Born and Wolf (October 13, 1999):

[TABLE]

with $\beta=m^{2}$, with $m$ being an integer number and $a_{m}$ and $b_{m}$ are integration constants.

Equation (64) is well known for spherical harmonics. Single-valued solutions to this equation exist when $\alpha=l(l+1)$ with ($l>|m|,$ integer). With this condition, the solution to (64) becomes

[TABLE]

Now we focus on the equation for the radial function (63), which may be rewritten as

[TABLE]

This second-order differential equation has two well-known solutions that are linearly independent: the regular function $F_{\ell}(kr_{g},kr)$ and the irregular function $G_{\ell}(kr_{g},kr)$. A regular function is so named because it is zero at $r=0$. Any solution to (63) may be chosen as linear combination of these two functions Messiah (1968); Thompson and Nunes (2009):

[TABLE]

where $F_{\ell}$ and $G_{\ell}$ are the Coulomb functions (discussed in Appendix D) and $c_{\ell}$ and $d_{\ell}$ are arbitrary constants.

According to (62), a particular integral $\Pi_{i}$ is obtained by multiplying together the functions given by (66), (67) and (69); we then obtain a general solution to (299). Collecting results for $\Phi(\phi)$, $\Theta(\theta)$ and $R(r)$, given by (66), (67), and (69), to the order of ${\cal O}(r_{g}^{2})$, the Debye potential has the form

[TABLE]

where $c_{l},d_{l},a_{m},b_{m}$ are arbitrary and yet unknown constants.

We must now determine these constants in such a way as to satisfy the boundary conditions. For this to be possible, one must be able to express the potentials $\,{}^{e}{\hskip-1.0pt}\Pi^{(i)}$ and $\,{}^{m}{\hskip-1.0pt}\Pi^{(i)}$ of the incident wave in a series of the from (70).

To proceed with the solution of the scattering problem, we consider the incident wave given by (40)–(41). Its properties should give us the partial wave amplitudes $c_{\ell}$ and $d_{\ell}$ in (70). To do this may not be straightforward, because these fields are singular at $\theta=0$ and cannot be written in terms of Legendre polynomials $P_{n}^{1}(\cos\theta)$ at all.

To determine $\,{}^{e}{\hskip-1.0pt}\Pi$ or $\,{}^{m}{\hskip-1.0pt}\Pi$, we use Eqs. (40)–(41) that describe the incoming wave and substitute them into (305)–(310). For example, for $D_{r}$ Eq. (40) yields

[TABLE]

where $\psi_{i}({\boldsymbol{\mathrm{r}}})$ is the incident scalar wave (23). Together with (58) (or the first part of (305)), after omitting the $e^{-i\omega t}$ factor, we obtain

[TABLE]

Our first problem, therefore, is to find an electromagnetic field, which for $r\rightarrow\infty,\theta\sim\pi$ has the same asymptotic behavior as the incident field given in (40), but which is regular everywhere, for all values of $\theta$ and $r$. Instead of using only a partial asymptotic solution representing the incident wave, $\psi_{i}({\boldsymbol{\mathrm{r}}})$, this field can be constructed using the full solution given by (20) and (22), for which (23) represents one of its asymptotic limits when $r\rightarrow\infty$:

[TABLE]

We may extend this to find the solution for the EM field in all regions by taking, instead of $\psi_{i}({\boldsymbol{\mathrm{r}}})$, the entire solution for $\psi$ from (73). Equation (72) indicates that

[TABLE]

is a suitable definition of the wanted regular field Herlt and Stephani (1975, 1976). The exact solution for $D_{r}$ based on (73) should differ from the incident wave (40) only for outgoing waves, the amplitudes of the incoming waves should be equal.

The function $\psi$ on the left-hand side of this equation may be expressed in the form of a differentiable series of Legendre polynomials Messiah (1968); Born and Wolf (October 13, 1999):

[TABLE]

where $F_{\ell}$ is the Coulomb function discussed in Appendix D. This representation is analogous to the following representation of a plane wave $\psi_{0}({\boldsymbol{\mathrm{r}}})=e^{ikz}$, given as

[TABLE]

where $j_{\ell}(kr)$ is the spherical Bessel function given by (277). Note, when $r_{g}\rightarrow 0$, one may see from (275) that function $\psi_{0}({\boldsymbol{\mathrm{r}}})$ is the limit of $\psi({\boldsymbol{\mathrm{r}}})$.

Using (76) and the identities

[TABLE]

we can write the left-hand side of (74) as

[TABLE]

This expression allows us to present a trial solution for $\,{}^{e}{\hskip-1.0pt}\Pi$ as a series of a form similar to (78), to order ${\cal O}(r_{g}^{2})$:

[TABLE]

Considering the asymptotic expansion of (79), we can substitute (78) and (79) into (74). Remembering that $F_{\ell}$ satisfies (68) and comparing coefficients, we obtain the relation

[TABLE]

The calculations for the magnetic potential, ${}^{m}{\hskip-1.0pt}\Pi$, are similar. In fact, in the vacuum, the solutions for the electric and magnetic potentials of the incident wave, ${}^{e}{\hskip-1.0pt}\Pi$ and ${}^{m}{\hskip-1.0pt}\Pi$, may be given in terms of a single potential $\Pi(r,\theta)$ as

[TABLE]

Therefore, by matching the general form for the Debye potentials (70) to the incident EM wave (75), we see that Maxwell’s equations (58)–(60) can only be satisfied by selecting $c_{\ell}=1$ and $d_{\ell}=0$, and also by choosing $m=1$, with $a_{1}=0$ for the magnetic potential, and $b_{1}=0$ for the electric potential. Thus, we have expressed both Debye potentials of the incident wave, ${}^{e}{\hskip-1.0pt}\Pi$ and ${}^{m}{\hskip-1.0pt}\Pi$, in the form of the series (70) by determining all the unknown constants. As a result, (81) represents an exact vacuum solution via Debye potentials for the EM field scattered by a gravitational monopole.

In the background of the metric (1), with $u$ from (17), the general solution of Maxwell’s equations (7)–(8) that corresponds to a monochromatic wave with the symmetry of a plane wave can be given in terms of a function $\Pi$, given by (81). Using this result in Eqs. (58)–(60) with the help of (305) we see that, in order to obtain the components of the EM field in a vacuum, we need to construct the following expressions Herlt and Stephani (1976):

[TABLE]

and insert them into

[TABLE]

This completes the solution for the EM field in a vacuum in the background of a spherically symmetric, static gravitational field represented by its Schwarzschild radius. However, the Sun has a physical boundary with a radius that is much larger than $r_{g}$. To account for this fact, we need to apply the fully absorbing boundary condition, as discussed in Sec. II.7.

III.2 Boundary conditions

As we discussed in Sec. II.7, the physical size of the Sun necessitates a proper treatment. Usually, this is done by selecting a form of the Debye potential for each of the regions in question, imposing the relevant boundary conditions, and matching the potentials on the boundary. We will follow a similar approach. First we note that, in order to match the potentials (81) to those of the incident and scattered waves, the latter must be expressed in a similar series form but with arbitrary coefficients. Only the function $F_{\ell}(kr_{g},kr)$ may be used in the expression for the potential, since $G_{\ell}(kr_{g},kr)$ is divergent at the origin. On the other hand, the scattered wave must vanish at infinity and the Hankel functions, $H^{+}_{\ell}(kr_{g},kr)$ (see Appendix D.2 for a discussion of the Hankel and Coulomb functions, their relationships and their relevant properties), will impart precisely this property. This function is suitable as a representation of the scattered wave. For large values of the argument $(kr)$, it behaves as $e^{ik(r+r_{g}\ln 2kr)}$ and the Debye potential will satisfy $\Pi\propto e^{ik(r+r_{g}\ln 2kr)}/r$ for large $r$. Thus, for distances $r\gg r_{g}$, the diffracted wave is spherical, with its center at the origin $r=0$. Accordingly, it will be used in the expression for the diffracted wave:

[TABLE]

To select the arbitrary coefficients $a_{\ell}$ we will use the fully absorbing boundary condition discussed in Sec. II.7. For this, we first consider the effective potential in Eq. (68) for the radial function $R_{\ell}$. We notice that a transition from small-$(kr)$ power law behavior to large-$(kr)$ oscillatory behavior occurs outside the classical turning point, which is the point where the effective potential in (68) vanishes, namely $1+2r_{g}/r-\ell(\ell+1)/(kr)^{2}={\cal O}(r_{g}^{2})$. Solving this quadratic equation, we determine the turning point

[TABLE]

As $r$ is positive, then with purely Newtonian (or, in nuclear scattering, Coulomb) and centrifugal potentials (68) there is only one turning point corresponding to the $+$ sign in (87). Classically, the turning point is at the distance of closest approach or at the impact parameter. These quantities are related in the same manner as the classical impact parameter $b_{0}$ is related to the quantum mechanical partial wave $\ell$ Thompson and Nunes (2009); Landau and Lifshitz (1989):

[TABLE]

To set the boundary conditions, we realize that rays with impact parameter $b_{0}\leq R_{\odot}$ are absorbed by the Sun. Thus, the fully absorbing boundary condition signifies that all the radiation intercepted by the body of the Sun is fully absorbed by it and no reflection or coherent reemission occurs. All intercepted radiation will be transformed into some other forms of energy, notably heat. Thus, we require that no scattered waves exist with impact parameter $b_{0}\ll R_{\odot}$ or, equivalently, for $\ell\leq kR_{\odot}$ It means that we need to subtract the scattered wave (86) from the incident wave for $\ell\leq kR_{\odot}$. In other words, to derive the solution for the Debye potential $\Pi^{\tt(I)}$ in the region outside the Sun (denoted by Latin superscript I), we set $a_{\ell}=-1$ in the expression for the scattering potential $\Pi^{({\tt s})}$ given by (86) and add to the expression for $\Pi_{\tt inc}$ from (81). This results in

[TABLE]

This is the second asymptotic boundary condition which is set on the “future infinity” light cone and deals with the fact that the physical boundary of the Sun is much larger than its Schwarzschild radius, $R_{\odot}\gg r_{g}$. This is in addition to the earlier condition that was established in “past infinity”, to fix the value for $\psi_{0}$ in (73).

We have thus obtained the Debye potential representing the total solution for the problem of diffraction of EM waves by a large spherical star. Solution (89) describes the EM field outside the Sun, which is our primary interest, and which we discuss next.

III.3 Exact solution for the Debye potentials

We observe that, in addition to the solution for the Debye potential in the form of the infinite series of partial waves (81), in a vacuum there exists an exact analytical solution for this quantity. To demonstrate this, we use the wave equation (300) written in the spherical coordinate system and present the expression for $D_{r}$ via derivatives with respect to $\theta$, as it was originally obtained in (298) and shown in (305), ultimately leading to (82). Then, from the two expressions for $D_{r}$ given by (85) and also by (74) with the $\exp(-\omega t)$ term reinstated, we obtain

[TABLE]

As a result, (90) yields the following equation to determine the Debye potential $\Pi$:

[TABLE]

We may now integrate this equation with respect to $\theta$ to obtain

[TABLE]

where $c(r)$ is the integrating constant. Integrating again from $\pi$ to $\theta$, we have

[TABLE]

Using (73) for $\psi$ and relying on the properties of the hypergeometric function from Appendix C, especially (215), we can evaluate the integral:

[TABLE]

By taking the integration constant to be

[TABLE]

we obtain the following expression for the Debye potential:

[TABLE]

which gives us the Debye potential of the incident wave in terms of the Coulomb wave function $\psi$, i.e., essentially in terms of the confluent hypergeometric function Herlt and Stephani (1975, 1976). This solution is always finite and is valid for any angle $\theta$.

As a result, the solution (96) for the Debye potential allows us to replace the first term in (89) and rewrite it as

[TABLE]

This is our main result. It contains all the information about the EM field around the Sun in all the regions of interest for the diffraction problem (see Fig. 4). We will evaluate the terms in this expression for each of these regions.

III.4 Solution to the diffraction problem and different regions

In order to understand the solution (97) that we obtained, we need more information on the second term in this expression. Considering the region outside the Sun, $r\gg r_{g}$, we may replace $H^{+}_{\ell}(kr_{g},kr)$ with its asymptotic expansion (271). Extending it to distances closer to the turning point, as derived in Appendix F and shown in (326), we obtain

[TABLE]

Next, we use the asymptotic representation for $P^{(1)}_{l}(\cos\theta)$ from Korn and Korn (1968):

[TABLE]

At this point, we may replace the sum in (98) with an integral:

[TABLE]

and evaluate this integral by the method of stationary phase. Note that the lower bound in this integral should be of the size of the Einstein radius of the lens. However, taking into account the physical dimensions of the Sun, such a detail is insignificant. Expression (100) shows that the $\ell$-dependent part of the phase has the structure:

[TABLE]

Therefore, the points of stationary phase where $d\varphi_{\pm}/d\ell=0$ are given by the following equation:

[TABLE]

with $\sigma_{\ell}$ taken from expression (265) where we formally replaced the sum with an integral, namely $\sum_{j=1}^{\ell}\rightarrow\int^{\ell}dj$.

If we take $\ell$ from the semiclassical approximation presented by (88), then for small angles $\theta$, Eq. (102) yields $\pm\sin\theta=2{r_{g}}/b_{0}-b_{0}/r+{\cal O}(r_{g}^{2})$. As a result, we see that the points of stationary phase satisfy the equation

[TABLE]

The potential $\delta(\Pi^{\tt(I)})$ from (100) contributes only if the points of stationary phase are within the interval $0\leq\theta\leq\pi$ and $1\leq\ell\leq kR_{\odot}$. As the largest impact parameter in (103) is set by the upper integration limit in (100), or $b_{0}^{\tt max}=R_{\odot}$, we see that this equation gives us the boundary of those regions influenced by $\delta(\Pi^{\tt(I)})$. This equation allows for a simple geometric and physical interpretation. We remember that the classical scattering orbit in a Newtonian potential is a hyperbola, described in polar coordinates $(\rho,\theta,\phi)$, starting at $\theta=\pi$, by Herlt and Stephani (1976); Thompson and Nunes (2009):

[TABLE]

which, based on the analysis in Appendix B.1, describes the geodesic path of the photon in the gravitational field of a monopole. From this, we see that the boundary in question coincides with the rays that are just grazing the Sun in the forward direction, $0\leq\theta\leq\frac{\pi}{2}$. Furthermore, for distances $z\leq z_{0}=R^{2}_{\odot}/2r_{g}$ (derived from (103) with $\theta=0$), one needs to take the plus sign in (103) and for distances beyond that point, $z\geq z_{0}$, the minus sign should be taken.

As a result, we established the boundary that separates three regions of interest (see Fig. 4 for details), namely:

i) For impact parameters $b_{0}\leq R_{\odot}$, the boundary conditions establish the shadow behind the Sun where no light from the source may appear;

ii) Impact parameters larger that the solar radius, $b_{0}>R_{\odot}$, correspond to regions of geometric optics where only one ray from a point source could pass through each point. The solution for the EM field in this region is given by the incident and scattered waves (40)–(41) and (49)–(50), correspondingly. However, as we discussed in Sec. II.6, the scattered wave is negligibly small everywhere in this region and offers practically no contribution;

iii) For distances beyond $z_{0}=R^{2}_{\odot}/(2r_{g})$, as we approach the optical axis, $\theta\rightarrow 0$, we enter the interference region where, in the immediate vicinity of the optical axis, the beam of extreme intensity is present.

Proper description of the EM field in this region requires a wave-theoretical treatment, which we develop next.

III.5 The electromagnetic field in the region of interference

We now consider the region of interference, i.e., the region in the immediate vicinity of the optical axis, $\theta\approx 0$, and at distances beyond $z\geq z_{0}=R_{\odot}^{2}/(2r_{g})$, so that the argument in (96) is small, namely $kr(1-\cos\theta)\ll 1$. We realize that in this region the second term in (89) produces no contribution and the EM field can be derived in its entirety from (96) Herlt and Stephani (1976). In addition, it can be shown by direct computation that the second term enclosed in round brackets in (96) can be neglected. The EM field and the Poynting vector due to it are orders of magnitude (factor of $(kr_{g})^{-1/2}$) smaller than those originating from the first term. The second term is important only near the axis $\theta=\pi$ where it serves to avoid a singularity. Thus the task that remains is the derivation of the Poynting vector of the field given by

[TABLE]

where, for convenience, and again following the logic of Herlt and Stephani (1976), we introduced the notation

[TABLE]

As we remember, $F[1]$ was first seen in (20) as a part of the solution of the time-independent Schrödinger equation for the scalar intensity of the EM wave, $\psi$. From (105) we see that $F[2]$ determines the properties of the Debye potential that corresponds to that solution.

In Appendices C.2 and C.3 we discuss the properties of these two functions and their behavior at small angles $\theta$ and also at large distances. Using the asymptotic behavior of $F[2]$ at large values of argument $k(r-z)\gg 1$ and $\psi_{0}$ from (73) and expressing $z=r\cos\theta$, we compute the asymptotic behavior of the Debye potential $\Pi$ from (105) as

[TABLE]

We can verify that the first term in (107) is the Debye potential corresponding to the incident wave, while the second term corresponds to the scattered wave. In fact, by substituting (107) into (82)–(85), after some algebra, we can see that the solution given by (107) yields results that are identical to the expressions for the incident and scattered fields given by (40)–(41) and (49)–(50), obtained earlier using different approach. Therefore, the exact solution for the Debye potential (105) may be used for any region describing the EM field.

Using the solution for the Debye potential $\Pi$ given by (105), we may now compute all the quantities in (82)–(84):

[TABLE]

By taking the asymptotic behavior of $F[1]$ and $F[2]$ from (233) and (237), correspondingly, together with $\psi_{0}$ from (73), substituting these into (108)–(110), and using the results in (85), we can verify that at large distances our solution gives the correct expression for each component of the incident (40)–(41) and scattered (49)–(50) EM waves. We can use the quantities (108)–(110) to compute the resultant EM field.

The solution (108)–(110) is valid for any angle and distance from the lens. However, for practical purposes, we are interested only in the small region on the optical axis just after the point where grazing rays intersect (see Fig. 4). We established earlier that, in the post-Newtonian approximation, the trajectories of light rays are governed by geodesic equations. These equations tell us that the focal line along which rays of light grazing the Sun intersect begins at $z_{0}=547.8$ AU. As was discussed in Herlt and Stephani (1976), beyond that point, the solar gravitational monopole forms a folded caustic (Fig. 5) that is characterized by a very high density of the EM field along the focal line, or optical axis. In the immediate vicinity of the optical axis $\rho\ll r_{g}$, the caustic is in the shape of a pencil-sharp beam. This region of the caustic, characterized by $0\lesssim\theta\ll\sqrt{2r_{g}/r}$, is where we direct our attention next.

III.6 Transformation to cylindrical coordinates

As argued in Herlt and Stephani (1976), for practical purposes it is convenient to introduce a cylindrical coordinate system $(\rho,\phi,z)$ instead of the spherical coordinates $(r,\theta,\phi)$. In the far field, $r\gg r_{g}$, this can be done by defining $R=ur=r+{r_{g}}/{2}+{\cal O}(r_{g}^{2})$ and introducing the coordinate transformations $\rho=R\sin\theta,z=R\cos\theta$, which, from (1), yield the line element:

[TABLE]

As a result, taking into account (85) and using the rules of vector transformations between curvilinear coordinates given by (152), for the metric (111) we have the following components of the EM field in cylindrical coordinates:

[TABLE]

where

[TABLE]

Using (108)–(110) for $\alpha$ and $\beta$, for a high-frequency EM wave (i.e., neglecting ${\cal O}((kr)^{-1})$ terms), we obtain

[TABLE]

We will use these results to study the properties of the EM field characterizing the diffraction of light by the SGL.

III.7 The electromagnetic field in the image plane

The components of the EM field in the cylindrical coordinate system $(\rho,\phi,z)$ are given by (112)–(114) with amplitudes given by (115)–(117). We note that at large distances from the Sun, we may neglect the terms $\sim r_{g}/r$ leading, in particular, to ${\boldsymbol{\mathrm{D}}}\simeq{\boldsymbol{\mathrm{E}}}+{\cal O}(r_{g}/r)$ and ${\boldsymbol{\mathrm{B}}}\simeq{\boldsymbol{\mathrm{H}}}+{\cal O}(r_{g}/r)$. Together with (115)–(117) and neglecting ${\cal O}((kr)^{-1})$ and ${\cal O}(r_{g}/r)$ terms (i.e., keeping only the largest terms), the physical components of the electric field take the form

[TABLE]

Similar expressions may be derived for the magnetic field ${\boldsymbol{\mathrm{H}}}$. Furthermore, in the immediate vicinity of the optical axis, $\rho\lesssim r_{g}$, we may use approximations for the functions $F[1]$ and $F[2]$ given by (252)–(253). For all practical applications, we may neglect terms containing $\theta^{2}$, not only because in the immediate vicinity of the optical axis $\rho\lesssim r_{g}$ and, thus, $\theta$ is very small; furthermore, the Bessel functions at those distances $\rho$ are also small. We are then left with the following solution for the EM field in the image plane:

[TABLE]

with $x=k^{2}rr_{g}(1-\cos\theta)$. Expressing $x$ in terms of cylindrical coordinates of (111) yields

[TABLE]

Using this result and $\theta=\rho/z+{\cal O}(\rho^{2}/z^{2})$, we can express the ratio in the second term of (121) as:

[TABLE]

These results allow us to present (121) in the form showing explicit dependence on all variables involved:

[TABLE]

Clearly, at the focal region of the SGL, when $z\geq z_{0}=R_{\odot}^{2}/2r_{g}=547.8$ AU, the factor in front of the $z$-components of the EM field, ${\hat{E}}_{z}$ and ${\hat{H}}_{z}$, is negligibly small. Thus, both of these components may be neglected, leaving only transverse components of the EM field on the image plane.

Solution (124) offers a good approximation for the EM field within a pencil-sharp beam in the very narrow vicinity of the optical axis, $\rho\lesssim r_{g}$; it is also quite accurate even for larger distances $\rho\sim 10^{2}\,r_{g}$. It shows that the EM field is distributed narrowly in the immediate region of the optical axis and falls off sharply as one moves away from it.

III.8 The Poynting vector in cylindrical coordinates

To consider the imaging properties of the SGL, we need to know the energy flux at the image plane, which is given by the Poynting vector. Components of the Poynting vector Landau and Lifshitz (1988); Logunov (1987) are given by (52). To compute ${\boldsymbol{\mathrm{S}}}$ in the cylindrical coordinate system, we use (112)–(114) and (115)–(117), and express the components of the Poynting vector as

[TABLE]

Averaging (125) over time and considering only high-frequency EM waves (i.e., neglecting ${\cal O}((kr)^{-1})$ terms), we get

[TABLE]

where the asterisk $({}^{*})$ denotes the complex conjugate. All properties of the diffraction field are encoded in these formulae (126)–(128). As noted in Herlt and Stephani (1976), extracting these properties is challenging because of the number of parameters that must be considered: the heliocentric distance $z$, the distance $\rho=z\theta$ from the axis $\theta=0$ in the image plane, the frequency of the wave $\omega$ and the telescope aperture.

Equations (254)–(255) allow us to present (126)–(128) up to the terms of $\propto\theta^{2}$:

[TABLE]

Using again the result (122) and $\theta=\rho/z+{\cal O}(\rho^{2}/z^{2})$, we can express the ratio in the second term of (131) as

[TABLE]

When a practical SGL is considered, this ratio is negligible. Therefore, the second term in (131) may be omitted.

Next, we consider the constant $\psi_{0}$ given by (73), for which the following is valid: $\psi^{2}_{0}=E_{0}^{2}\,e^{\pi kr_{g}}\Gamma(1-ikr_{g})\Gamma(1+ikr_{g}).$ Using the properties of the gamma function Abramowitz and Stegaun (1965), we have $\Gamma(1-ikr_{g})\Gamma(1+ikr_{g})={\pi kr_{g}}/{\sinh\pi kr_{g}},$ which for $\psi_{0}^{2}$ results in the following expression:

[TABLE]

Given the fact that in the focal region of the SGL, the ratio $r_{g}/r\ll 1$ is very small, the terms in (129)–(131) that include this ratio may also be omitted. As a result, using (122) for the argument of the Bessel function, we can present the components of the Poynting vector (129)–(131) in the following most relevant form:

[TABLE]

with ${\bar{S}}_{\rho}={\bar{S}}_{\phi}=0$ for any practical purposes. Note that in the case when $r_{g}\rightarrow 0$, the Poynting vector reduces to its Euclidean spacetime vacuum value, namely ${\bar{S}}\rightarrow{\bar{\boldsymbol{\mathrm{S}}}}_{0}=(0,0,(c/8\pi)E_{0}^{2})$, which may de deduced from (53) by taking $r_{g}=0$. Note that in the limit $\lambda/r_{g}\rightarrow 0$, (134) corresponds to the geometric optics approximation which yields a divergent intensity of light on the caustic.

Result (134) completes our derivation of the wave-theoretical description of light propagation in the background of a gravitational monopole. The result that we obtained extends previous derivations that are valid only on the optical axis (e.g., Eshleman (1979)) to the neighborhood of the focal line and establishes the structure of the EM field in this region. As such, it presents a useful wave-theoretical treatment of focusing light by a spherically symmetric mass, which is of relevance not only for the SGL discussed here but also for microlensing by objects other than the Sun.

IV Towards a Solar Gravitational Telescope

We now have all the tools necessary to establish the optical properties of the SGL in the region of interference, i.e., at heliocentric distances $z\geq z_{0}=R^{2}_{\odot}/2r_{g}=547.8$ AU on the optical axis. First, given the knowledge of the Poynting vector in the image plane (134), we may define the monochromatic light amplification of the lens, ${\mu}$, as the ratio of the magnitude of the time-averaged Poynting vector of the lensed EM wave to that of the wave propagating in empty spacetime ${\boldsymbol{\mathrm{\mu}}}={\boldsymbol{\mathrm{\bar{S}}}}/|{\boldsymbol{\mathrm{\bar{S}}}}_{0}|$, with $|{\bar{\boldsymbol{\mathrm{S}}}}_{0}|=(c/8\pi)E_{0}^{2}$. The value of this quantity is then given by

[TABLE]

As evident from (134), we see that the largest amplification of the SGL occurs along the $z$ axis. The other components of the Poynting vector are negligible.

We now consider the light amplification of the SGL in the focal region. Figure 6 shows the resulting Airy pattern (i.e., the point spread function or PSF) of the SGL from (135). Due to the presence of the Bessel function of the zeroth order, $J_{0}^{2}(2\sqrt{x})$, the PSF falls off more slowly than traditional PSFs, which are proportional to $J^{2}_{1}(2\sqrt{x})/x^{2}$, as seen in Fig. 7. Thus, a non-negligible fraction of the total energy received at the image plane of the SGL is present in the side lobes of its PSF. This indicates that for image processing purposes, one may have to develop special deconvolution techniques beyond those that are presently available (e.g., Koechlin et al. (2005, 2006)), which are used in modern microlensing surveys. Most of these techniques rely on raytracing analysis and typically are based on geometric optics approximation.

Furthermore, the light amplification $\mu$ weakly depends on the distance from the Sun. For practical purposes, it is easier to show this property by plotting the gain of the SGL, $g$, which is related to light amplification as $g(\lambda,z)=10\log_{10}\mu(\lambda,z)$. Figure 8 plots the gain of the SGL at two heliocentric distances $z=600$ AU and $1,000$ AU for two wavelengths $\lambda=1.0~{}\mu$m and 2.0 $\mu$m.

We may express the argument of the Bessel function in (135) in terms of the quantities of interest, namely the heliocentric distance along the optical axis $z$, the distance in the image plane $\rho$ (as measured from the optical axis), and the impact parameter $b_{0}$. With the help of (122) we have:

[TABLE]

where $\alpha_{0}={2r_{g}}/{R_{\odot}}=8.490\times 10^{-6}~{}{\rm rad}=1.751^{\prime\prime}$ is the angle of deflection by the SGL for the light rays just grazing the Sun. Given numerical values of various quantities involved, we obtain

[TABLE]

or,  equivalently,

[TABLE]

This result clearly shows the dependence of the SGL’s light amplification on the observing wavelength, $\lambda$, the distance along the focal line, ${z}$, and the distance from the focal line in the image plane, $\rho$. The value of maximum amplification of the SGL, $\mu_{0}=4\pi^{2}r_{g}/\lambda$, is independent of ${z}$. For optical wavelengths, this amounts to $\mu_{0}\sim 1.2\times 10^{11}$, giving the SGL its enormous light amplification. For small deviations from the optical axis, the light amplification (135) drops sharply, as seen in Fig. 6, but the overall envelope decreases more slowly than that of a traditional PSF (Fig. 7).

The ability of a lens to resolve detail is ultimately limited by diffraction. Light coming from a point source diffracts through the lens aperture, forming a diffraction pattern in the image plane known as an Airy pattern (see Fig. 6). The angular radius of the central bright lobe, called the Airy disk, is measured from the center to the first null. Therefore, we define the resolution of the SGL using the location where $J_{0}(2\sqrt{x})=0$, which is satisfied for the value of the argument of $2\sqrt{x}\approx 2.40483$. We can then solve (136) for $\theta_{\tt SGL}=\rho/z$:

[TABLE]

where $D_{\odot}=2R_{\odot}$ is the solar diameter. For the wavelength $\lambda=1~{}\mu$m, the resolution of the SGL at ${z}_{0}=547.8~{}{\rm AU}$ is $\theta_{0}\approx 5.50\times 10^{-16}~{}{\rm rad}=0.11~{}{\rm nas}$. The resolution increases with ${z}$ as $\theta_{0}\sqrt{{z}_{0}/{z}}$ as

[TABLE]

For an exoplanet situated at the distance $z_{\tt p}$ from the Sun, the angular resolution (139) translates into resolvable surface features of $\delta\rho_{\tt SGL}=\theta_{\tt SGL}z_{\tt p}$, which improves with heliocentric distance as

[TABLE]

Depending on the impact parameter, the deflection angle of the SGL is given as $\alpha=2r_{g}/b_{0}=\alpha_{0}(R_{\odot}/b_{0})$. Rays with impact parameter $b_{0}$ will intersect the optical axis at the distance of $z=b_{0}/\alpha=547.8\,(b_{0}/R_{\odot})^{2}$ AU. In the pencil-sharp region along the focal line the amplification (135) of the SGL stays nearly constant well beyond 2,500 AU, while its angular resolution (140) increases by a factor of $\sim 1/\sqrt{5}$ in the same range of heliocentric distances.

Across the image plane, the amplification oscillates quite rapidly. For small deviations from the optical axis, $\theta\approx\rho/{z}$. Using this relation in (140), we see that the first zero occurs quite close to the optical axis:

[TABLE]

(Note in (142) the inverse ratio of ${z}$ vs. ${z}_{0}$ and $b_{0}$ vs $R_{\odot}$.) Equation (142) favors larger wavelengths and larger heliocentric distances or, similarly, impact parameters.

Thus, we have established the basic optical properties of the solar gravitational lens. By achromatically focusing light from a distant source Kraus (1986); Turyshev and Andersson (2003), the SGL provides a major brightness amplification and extreme angular resolution. Specifically, from (135) for $\lambda=1~{}\mu$m, we get a light amplification of the SGL of $\mu\simeq 1.2\times 10^{11}$, corresponding to a brightness increase by $\delta{\rm mag}=2.5\ln\mu=27.67$ stellar magnitudes in case of perfect alignment. Furthermore, (140) gives us the angular resolution of the SGL of $\theta_{\tt SGL}\simeq 1.1\times 10^{-10}$ arc seconds.

We note that if the diameter of the telescope $d_{0}$ is larger than the diffraction limit of the SGL (i.e., larger than the diameter of the first zero of the Airy pattern), it would average the light amplification over the full aperture. Such an averaging will result in the reduction of the total light amplification. To estimate the impact of the large aperture on light amplification, we average the result (135) over the aperture of the telescope:

[TABLE]

For an aperture of $d_{0}=1$ m at $z=600$ AU, this results in the reduction in light amplification by a factor of $0.025$, leading to the effective light amplification of $\bar{\mu}_{z}=2.87\times 10^{9}$ (i.e., 23.65 mag), which is still quite significant. The effect of the large aperture is captured in Fig. 9, where we plot the behavior of each of the two terms in curly braces in (143) and also their sum. Although each term oscillates and reaches zero, their sum never becomes zero.

As seen from a telescope at the SGL, light from a distant target fills an annulus at the edge of the Sun, forming the Einstein ring. At a distance $z$ on the focal line, an observer looking back at the Sun will see the Einstein ring with an angular size that is given by $\alpha_{\tt ER}=2b_{0}/z=4r_{g}/b_{0}$. Using this equation, we determine the angular size of the ring as

[TABLE]

A telescope with aperture $d_{0}$, placed at the heliocentric distance $z$ on the optical axis, receives light from a family of rays with different impact parameters with respect to the Sun, ranging from $b_{0}$ to $b_{0}+\delta b_{0}$. Using (144), these rays are deflected by different amounts given as $\alpha_{1}=(b_{0}+{\textstyle\frac{1}{2}}d_{0})/z=\alpha_{0}R_{\odot}/(b_{0}+{\textstyle\frac{1}{2}}d_{0})$, for one edge of the aperture, where $\alpha_{0}=2r_{g}/R_{\odot}$, and $\alpha_{2}=(b_{0}+\delta b_{0}-{\textstyle\frac{1}{2}}d_{0})/z=\alpha_{0}R_{\odot}/(b_{0}+\delta b_{0}-{\textstyle\frac{1}{2}}d_{0})$, for the other edge. Taking the ratio of $\alpha_{2}/\alpha_{1},$ we can determine the relation between $\delta b_{0}$ and the telescope diameter, $d_{0}$, which, to first order, is given as $\delta b_{0}=d_{0}$.

As a result, the area of the Einstein ring that is seen by the telescope with aperture $d_{0}$, to first order, is given by $A_{\tt ER}=\pi((b_{0}+\delta b_{0})^{2}-b_{0}^{2})\simeq 2\pi b_{0}d_{0}$. For different impact parameters the area behaves as

[TABLE]

Therefore, the magnifying power of a 1 m telescope placed at heliocentric distance $z$ on the focal line of the SGL is equivalent to a telescope with diameter of $D=2\sqrt{2b_{0}d_{0}}=74.6\,({b_{0}/R_{\odot}})^{\frac{1}{2}}$ km or, in terms of the heliocentric distances, it is given as $D=74.6\,(z/z_{0})^{\frac{1}{4}}$ km, which is a weak function of the observer’s position on the focal line.

To image an exoplanet, observing this annulus with thickness of $\delta b_{0}=d_{0}$ is, of course, the primary objective. A diffraction-limited 1 m telescope would have a resolution of $\delta\theta=\lambda/d_{0}=0.21^{\prime\prime}$ at $\lambda=1~{}\mu{\rm m}$. The thickness of the Einstein ring from the heliocentric distance of $z=600$ AU is $d_{0}/z=2.30$ nas. Thus, although the thickness of the Einstein ring is unresolved by the telescope at the SGL, the ring itself is well-resolved and can be used for imaging purposes. In fact, the entire circumference of the ring at the same distance of $z=600$ AU has the length of $L_{\tt ER}=2\pi b_{0}/z=10.05^{\prime\prime}\,(b_{0}/R_{\odot})$, and it is resolved with $L_{\tt ER}/\delta\theta\sim 48.7\,(b_{0}/R_{\odot})$ resolution elements. Thus, the ring could be used to provide information on a particular surface area on the target exoplanet. By sampling various parts of the ring, we will be able to collect data relevant to that particular surface area on the exoplanet.

Considering the plate scale: an Earth-sized exoplanet at $z_{\tt ep}=30$ pc away from the Sun, when imaged from the focal region of the SGL at heliocentric distance of $z\sim$600 AU, has the the image size of $2R_{\oplus}z/z_{\tt ep}\sim 1,238~{}{\rm m}\,(z/600\,{\rm AU})(30\,{\rm pc}/z_{\tt ep})$. A single telescope would have to traverse this area in the immediate vicinity of the focal line to scan the image of the exoplanet. Such a scaling law suggests that to image this object with $\sim 10^{3}\times 10^{3}$ pixels, the telescope would need to move in the image plane from pixel to pixel, each of which has the size of $\sim 1.2$ m. Each surface element resolved on the surface of the exoplanet would form its own Einstein ring around the Sun. However, because of the properties of the PSF of the SGL (which has prominent side lobes, as seen Fig. 7), the total flux within each Einstein ring corresponding to a particular surface element would also have contributions (in the form of Einstein arcs) from adjacent surface elements. Therefore, to form a reliable image of an exoplanet’s topography, multiple such images must be deconvoluted. This can be accomplished as the properties of the Sun and, thus, of the SGL are well understood.

Considering a realistic mission to the SGL to image a preselected target, one would have to consider the effects of the proper motion of the host star with respect to the Sun, as well as orbital dynamics of the target exoplanet and its diurnal rotation. Even if these factors are accounted for by a trajectory design and raster scan in the image plane, the exoplanet may also change as it is being scanned, due to changes in illumination, seasonal changes, cloud cover, the presence of one or more natural satellites and other factors; therefore, image deconvolution must also take place in the temporal dimension, possibly aided with reasonable models of periodic changes in appearance.

This interesting problem set must be addressed before exoplanet imaging using the SGL can become reality. Nonetheless, the potential benefits of a solar gravitational telescope (SGT) are well considered in comparison with the parameters of a comparable diffraction-limited optical telescope. Given the very small angular diameter ($\sim 1.4\times 10^{-11}$ rad) of an Earth-like planet at 30 pc, obtaining a single-pixel image would require a diffraction-limited telescope with an aperture of $\sim$74.6 km. To match the magnifying power of the SGL and obtain an image at a resolution of a thousand linear pixels, a telescope aperture of $4\times 10^{5}$ km ($\sim 16R_{\oplus}$) would be needed. Building an optical imaging interferometer with such a set of baselines is not feasible. At the same time, a mission to the SGL offers access to unique conditions needed for direct imaging of an exoplanet. Perhaps, it is the time we start taking the SGL seriously.

V Discussion and Conclusions

In this paper, we considered the propagation of EM waves in the gravitational field of the Sun, which is represented by the Schwarzschild monopole taken within the first post-Newtonian approximation of the general theory of relativity. We have developed a wave-theoretical treatment for light diffraction in the field of a static gravitational mass monopole and considered the case of a monochromatic EM wave coming from a point source at a large distance from the monopole. We obtained a solution for the EM field everywhere around the lens and especially in the immediate vicinity of its focal line, where the geometric optics leads to diverging results. As anticipated, because of wave effects in the focal region, our wave-optical treatment is immune to singularities, allowing us to describe the optics of the SGL and understand its image formation properties. As such, in contrast to models based purely on geometric optics, our approach allows us to consider practical questions related to the design of a SGT, in part by permitting the use of traditional tools of telescope design. The results that we obtained allow us to compute the PSF, resolution and FOV, as well as the evolution of these quantities at various heliocentric distances along the focal line. These will help improve our understanding of the unique properties of the SGL for imaging and spectroscopic investigations.

Our presentation is streamlined, taking full advantage of the weak-field gravity in the solar system. We also benefit from the tools and techniques borrowed from nuclear physics, specifically from the physics of scattering in a Coulomb field. Our approach can be extended to include higher-order solar gravity multipoles, if needed. We find that the formalism for Coulomb scattering from the nuclear physics literature is directly applicable. However, whereas nuclear particle physics studies focus on the scattering of scalar particles, we were able to develop the formalism required to describe the scattering of a vector EM field in the post-Newtonian approximation of the solar gravitational field.

Our results represent the first step towards developing a comprehensive theory of image formation by the SGL and the tools needed for mission design, data collection and processing, and ultimate image deconvolution Gaudi and Petters (2002a, b); Gould (2001); Wambsganss (1998). Several effects of gravitational and dynamical areas will require further analysis. In particular:

i) distinguishing the bright solar disk from the annulus of an Einstein ring, and the constraints it places on the performance of the SGT;

ii) effects due to the solar corona and solar plasma on light propagation;

iii) effects due to solar oblateness and solar rotation on the spatial and temporal properties of the caustic formed by the SGL;

iv) effects of reflex motion of the Sun with respect to the solar system’s barycentric coordinate reference system due to the presence of the giant gaseous planets in the solar system;

v) effects of proper motion of the exoplanet’s parent star, orbital motion of the planet around the barycenter of its planetary system, diurnal rotation of the planet, orientation of its axis of rotation, precision and nutation;

vi) temporal changes in the targeted planet’s appearance due to changing illumination, varying cloud cover, changes in atmospheric chemistry, varying surface features (ice cover, vegetation), varying illumination by its host star, and eclipses due to any satellites.

Some of these aspects will be addressed in the upcoming study of a mission to the SGL that is to be conducted at JPL Turyshev et al. (2017). The results of this study will be available elsewhere.

Concluding, we emphasize that our present understanding of the properties of the SGL and its value for imaging and spectroscopy is about at the same level as we knew gravitational waves back in the 1970s. At that time, the physics of gravitational waves was already well understood, but the technology needed for their detection was a long way in the future. That “future” for the research in gravitational waves came at the centennial of general relativity with the results of the first direct detection of the gravitational waves reported by the LIGO team Abbott et al. (2016). It is our hope and desire that by the theory’s sesquicentennial, we will be in possession of a fully developed set of technologies as well as the spacecraft, instruments, and data analysis tools required to collect data and present us with high-resolution imaging and spectroscopy of habitable exoplanets, relying on the physics of the solar gravitational lens.

Acknowledgements.

We thank L.D. Friedman for suggesting us to consider the problem and continuing interest through the entire research process. We are also thankful to M.V. Sazhin and M. Shao for their interest, support and encouragement during the work and preparation of this manuscript. This work was performed at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration.

Appendix A Three-dimensional metric and $(3+1)$ decomposition

We summarize basic rules for vector transformations and differential operators in curvilinear coordinates, for convenience and also to introduce the notations used throughout the present paper.

Following Landau and Lifshitz (1988) (see §84), we consider a generic interval and its $3+1$ decomposition:

[TABLE]

where the three-dimensional metric $\kappa_{\alpha\beta}$ is given as:

[TABLE]

If $g_{mn}$ is diagonal, so is $\kappa_{\alpha\beta}$. In the following, we assume a diagonal metric. We consider the standard basis Korn and Korn (1968) with unit basis vectors $\boldsymbol{\mathrm{i}}_{1}(x^{1},x^{2},x^{3})$, $\boldsymbol{\mathrm{i}}_{2}(x^{1},x^{2},x^{3})$, $\boldsymbol{\mathrm{i}}_{3}(x^{1},x^{2},x^{3})$ respectively directed along the coordinates $x^{1},x^{2},x^{3}$. When $\kappa_{\alpha\beta}$ is diagonal, these basis vectors form an orthonormal basis ($\boldsymbol{\mathrm{i}}_{\alpha}\cdot\boldsymbol{\mathrm{i}}_{\beta}=\delta_{\alpha\beta}$).

Components of a vector $\boldsymbol{\mathrm{F}}$ in this basis are defined by $\hat{F}_{\alpha}=(\boldsymbol{\mathrm{F}}\cdot\boldsymbol{\mathrm{i}}_{\alpha})$ (no summation), such that

[TABLE]

We now form the covariant basis as

[TABLE]

and the corresponding contravariant basis as

[TABLE]

where $[\boldsymbol{\mathrm{a}}\boldsymbol{\mathrm{b}}\boldsymbol{\mathrm{c}}]=(\boldsymbol{\mathrm{a}}\cdot[\boldsymbol{\mathrm{b}}\times\boldsymbol{\mathrm{c}}])$ represents the vector triple product.

We obtain the covariant components of a vector $\boldsymbol{\mathrm{F}}$ as $F_{\alpha}=({\bf F}\cdot{\boldsymbol{\mathrm{e}}}_{\alpha})$ and the contravariant components as $F^{\alpha}={\bf F}\cdot{\boldsymbol{\mathrm{e}}}^{\alpha}$. Consequently,

[TABLE]

The expressions to transform $F_{\alpha}$ from coordinates $\xi_{m}$ with Lamé coefficients $h_{m}$ to $\xi_{n}^{\prime}$ with Lamé coefficients $h_{n}^{\prime}$ are given as (see Chapter 1.3 in Morse and Feshbach (1953)):

[TABLE]

where for an orthonormal coordinate systems endowed with the diagonal 3-metric $\kappa_{\alpha\beta}$ (147), we have $h_{\alpha}=\sqrt{\kappa_{\alpha\alpha}}$.

The differential operators ${\rm grad}_{\kappa}\psi=\boldsymbol{\mathrm{\nabla}}_{\kappa}\psi,~{}{\rm div}_{\kappa}\boldsymbol{\mathrm{F}}=(\boldsymbol{\mathrm{\nabla}}_{\kappa}\cdot\boldsymbol{\mathrm{F}}),~{}{\rm curl}_{\kappa}\boldsymbol{\mathrm{F}}=[\boldsymbol{\mathrm{\nabla}}_{\kappa}\times\boldsymbol{\mathrm{F}}],$ and $\Delta_{\kappa}\psi=(\boldsymbol{\mathrm{\nabla}}_{\kappa}\cdot\boldsymbol{\mathrm{\nabla}}_{\kappa})\psi$ in orthonormal coordinate systems endowed with the diagonal 3-metric $\kappa_{\alpha\beta}$, Eq. (147), are given as Morse and Feshbach (1953); Korn and Korn (1968); Landau and Lifshitz (1988):

[TABLE]

Appendix B Light propagation in weak and static gravity

B.1 Geodesics in weak and static gravity

To investigate the propagation of light in the vicinity of the Sun, we consider the metric (1). We represent the trajectory of a photon as

[TABLE]

where $k^{\alpha}$ is the unit vector in the unperturbed direction of photon’s propagation and $x^{\alpha}_{\tt G}(t)$ is the post-Newtonian term. We define the four-dimensional wave vector in a curved space-time as usual:

[TABLE]

where $\lambda$ is the parameter along the ray’s path and $\kappa^{\alpha}={dx^{\alpha}}/{dx^{0}}$ is the unit vector in that direction, i.e., $\kappa_{\epsilon}\kappa^{\epsilon}=-1$ (do not confuse $\kappa^{\alpha}$ with the three-dimensional metric $\kappa_{\alpha\beta}$ in 147). From (157) we see that the unit vector $\kappa^{\alpha}$ may be represented as $\kappa^{\alpha}=k^{\alpha}+k^{\alpha}_{\tt G}(t)+{\cal O}(G^{2}),$ where $k^{\alpha}_{\tt G}(t)=dx^{\alpha}_{\tt G}/{dx^{0}}$ is the post-Newtonian perturbation. The wave vector obeys the geodesic equation: ${dK^{m}}/{d\lambda}+\Gamma^{m}_{kl}K^{m}K^{l}=0,$ which yields

[TABLE]

Equation (159) is an integral of motion due to energy conservation. Indeed, we can present it as

[TABLE]

Therefore, in the static field energy is conserved, and we have the following integral of motion:

[TABLE]

where $x^{0}_{\tt G}(\lambda)$ is the post-Newtonian correction. We recall that the wave vector $K^{m}$ is a null vector, which, to first order in $G$ and with $K^{0}=k^{0}+{\cal O}(G)$ yields K_{m}K^{m}=0=(k^{0})^{2}\big{(}1+\gamma_{\epsilon\beta}k^{\epsilon}k^{\beta}+{\cal O}(G)\big{)}. Then, Eq. (160) becomes

[TABLE]

We can now represent (163) in terms of derivatives with respect to time $x^{0}$. First we have

[TABLE]

Substituting (164) into (163) and using (159), we have

[TABLE]

Remember that for light $ds^{2}=0$. Then, from the fact that it moves along the light cones, the following expression is valid $g_{mn}({dx^{m}}/{dx^{0}})({dx^{n}}/{dx^{0}})=0=1+k_{\epsilon}k^{\epsilon}+{\cal O}(G)$, which for (165) yields

[TABLE]

We begin by examining the Newtonian part of (157) and representing it as

[TABLE]

Following Kopeikin (1997); Kopeikin et al. (2011), we define $b_{0}^{\alpha}\equiv{\boldsymbol{\mathrm{b}}}_{0}=[[{\boldsymbol{\mathrm{k}}}\times{\boldsymbol{\mathrm{x}}}_{0}]\times{\boldsymbol{\mathrm{k}}}]+{\cal O}(G)$ to be the impact parameter of the unperturbed trajectory of the light ray. The vector ${\boldsymbol{\mathrm{b}}}_{0}$ is directed from the origin of the coordinate system toward the point of the closest approach of the unperturbed path of light ray to that origin. We also introduce the parameter $\ell=\ell(t)$ as follows:

[TABLE]

These quantities allow us to rewrite (167) as

[TABLE]

The following relations hold:

[TABLE]

They are useful for presenting the results of integration of the light ray equations in different forms. Clearly, when the coordinate system oriented along the initial direction of the ray’s path, then $\ell=({\boldsymbol{\mathrm{k}}}\cdot{\boldsymbol{\mathrm{x}}})=z$.

Below, we focus our discussion on the largest contribution to the gravitational deflection of light: that due to the field produced by a monopole. In this case, the Newtonian potential may be given by $c^{-2}U({\boldsymbol{\mathrm{r}}})={r_{g}}/{2r}+{\cal O}(r^{-3},c^{-4}),$ where $r_{g}=2GM/c^{2}$ is the Schwarzschild radius of the source. Therefore, the quantity $u$ in (1) has the form

[TABLE]

If needed, one can account for the contribution of the higher-order multipoles using the tools developed in Kopeikin et al. (2011); Turyshev et al. (2014).

Limiting our discussion to the monopole given by (171), we have $c^{-2}\partial^{\alpha}U=-(r_{g}/2r^{2})\partial^{\alpha}r+{\cal O}(G^{2},r^{-4})$. We recall that $\partial^{\alpha}r=\partial^{\alpha}\sqrt{-x_{\epsilon}x^{\epsilon}}=-x^{\alpha}/r.$ Then, $c^{-2}\partial^{\alpha}U=(r_{g}/2r^{3})x^{\alpha}+{\cal O}(G^{2},r^{-4})$. In this case, Eq. (166) takes the form:

[TABLE]

Making the substitution $d/dx^{0}=d/d\ell$, we have the following equation:

[TABLE]

We integrate (173) from $-\infty$ to $\ell$ to get the following result:

[TABLE]

or, equivalently, with the help of (168)–(169) we have the following expression for the wave vector $\kappa^{\alpha}$ from (158):

[TABLE]

We may now integrate (174) from $\ell_{0}$ to $\ell$ to obtain

[TABLE]

which results in

[TABLE]

or, equivalently, substituting $\ell$ and $r$ from (168)–(169), we have

[TABLE]

Therefore, the trajectory of a photon in a static weak gravitational field is described by (177), while the direction of its wave vector $\kappa^{\alpha}={dx^{\alpha}}/{dx^{0}}$ is given by (175). For a radial light ray given by $k^{\alpha}=x^{\alpha}_{0}/r_{0}=n^{\alpha}_{0}$ and $b_{0}=0$, then Eqs. (175) and (178) become

[TABLE]

The solutions given by Eqs. (178) and (180) describe the motion of a photon along a geodesic in the post-Newtonian approximation in the static spacetime of a monopole. While Eq. (178) describes the motion along an arbitrary geodesic, Eq. (180) deals only with radial propagation of light.

B.2 Geometric optics approximation for the wave propagation in the vicinity of a massive body

In geometric optics, the phase $\varphi$ is a scalar function, a solution to the eikonal equation Fock (1959); Landau and Lifshitz (1988); Kopeikin (2009); Kopeikin et al. (2011):

[TABLE]

Given the wave vector $K_{m}=\partial_{m}\varphi$, and its tangent $K^{m}=dx^{m}/d\lambda=g^{mn}\partial_{n}\varphi$ where $\lambda$ is an affine parameter, we note that (181) states that $K^{m}$ is null ($g_{mn}K^{m}K^{n}=0$), thus

[TABLE]

Equation (181) can be solved by assuming an unperturbed solution that is a plane wave:

[TABLE]

where $\varphi_{0}$ is an integration constant and, to Newtonian order, $\underline{k}^{m}=(k^{0},k^{\alpha})=k_{0}(1,{\boldsymbol{\mathrm{k}}})$, where $k_{0}=\omega/c$, is a constant null vector of the unperturbed photon trajectroy, $\gamma_{mn}\underline{k}^{m}\underline{k}^{n}={\cal O}(G)$; $\varphi_{G}$ is the post-Newtonian perturbation of the eikonal. The wave vector $K^{m}(t,{\boldsymbol{\mathrm{x}}})$ then also admits a series expansion in the form

[TABLE]

where $k^{m}_{\tt G}(t,{\boldsymbol{\mathrm{x}}})=\gamma^{mn}\partial_{n}\varphi_{\tt G}(t,{\boldsymbol{\mathrm{x}}})$ is the first order perturbation of the wave vector. Substituting (183) into (181) and defining $h^{mn}=g^{mn}-\gamma^{mn}$ with $g_{mn}$, we obtain an ordinary differential equation to for $\varphi_{\tt G}$:

[TABLE]

where ${d\varphi_{\tt G}}/{d\lambda}=K_{m}\partial^{m}\varphi$. Similarly to (157), to Newtonian order, we represent the light ray’s trajectory as

[TABLE]

and substituting a monopole potential characterized by the Schwarzschild radius $r_{g}$ for $U$, we obtain

[TABLE]

Representation (186) allows us to express the Newtonian part of the wave vector $K^{m}$, as given by (184), as K^{m}={dx^{m}}/{d\lambda}=k^{0}\big{(}1,{\boldsymbol{\mathrm{k}}}\big{)}+{\cal O}(G), where $k^{0}$ is immediately derived as $k^{0}={cdt}/{d\lambda}+{\cal O}(G)$ and $|{\boldsymbol{\mathrm{k}}}|=1$. Keeping in mind that $\underline{k}^{m}$ is constant, we establish an important relationship:

[TABLE]

which we use to integrate (187). As a result, in the body’s proper reference frame Turyshev et al. (2013); Turyshev and Toth (2015), we then obtain

[TABLE]

which, for a radial light ray characterized by $k^{\alpha}=x^{\alpha}_{0}/r_{0}=n^{\alpha}_{0}$ (similarly to (180)), yields

[TABLE]

It is worth pointing out that the results obtained here for the phase of an EM wave (189) and (190) are equivalent to those obtained in the preceding section obtained for the geodesic trajectory of a photon (178) and (180).

B.3 Local basis vectors

In Sec. II.4 we introduced the local basis vectors ${\boldsymbol{\mathrm{\kappa}}}={\boldsymbol{\mathrm{K}}}/|{\boldsymbol{\mathrm{K}}}|,{\boldsymbol{\mathrm{\pi}}}=[{\boldsymbol{\mathrm{\kappa}}}\times{\boldsymbol{\mathrm{n}}}]/|[{\boldsymbol{\mathrm{\kappa}}}\times{\boldsymbol{\mathrm{n}}}]|$ and ${\boldsymbol{\mathrm{\epsilon}}}=[{\boldsymbol{\mathrm{\pi}}}\times{\boldsymbol{\mathrm{\kappa}}}]$. These vectors are very convenient to develop the results in this paper. In this appendix, we express these vectors in various coordinates with accuracy to the order of ${\cal O}(r_{g}^{2})$. We do that by using an expression for the trajectory of the photon (178) and its phase $\varphi$, (189) or, similarly, (23). We recognize from (23) that for a wave coming from $-\infty$ along the $z$-axis, the time-independent part of the phase with $({\boldsymbol{\mathrm{k}}}\cdot{\boldsymbol{\mathrm{r}}})=z$ has the form:

[TABLE]

From the definition for the wave vector, $K_{\alpha}=\partial_{\alpha}\varphi$, and with the help of (170), we have

[TABLE]

The covariant wave vector $K^{\alpha}$ is given as:

[TABLE]

From (175), we have

[TABLE]

Also, defining $K_{0}=k_{0}=\omega_{0}/c$ (see Landau and Lifshitz (1988)) in a static field, we have

[TABLE]

Therefore, collecting all the terms we have the following expression for $K^{\alpha}$:

[TABLE]

We can verify that the following relations hold:

[TABLE]

Next, we use expression (178) for the position vector of a photon on its trajectory, written as

[TABLE]

where $\ell$ and $r$ are given by (168) and (169), correspondingly. Expressions (196) and (198) allow us to compute all local vectors for a ray moving in the plane formed by ${\boldsymbol{\mathrm{k}}}$ and ${\boldsymbol{\mathrm{r}}}$ vectors:

[TABLE]

In the Cartesian coordinate system $(x,y,z)$ used to develop (36)–(37), remembering that the impact parameter has the form ${\boldsymbol{\mathrm{b}}}_{0}=[{\boldsymbol{\mathrm{k}}}\times[{\boldsymbol{\mathrm{r}}}\times{\boldsymbol{\mathrm{k}}}]]+{\cal O}(r_{g})=(x,y,0)+{\cal O}(r_{g})$, we present (199)–(201) in the following convenient form:

[TABLE]

The local basis vectors (202)–(204) represent the right-handed set of orthonormal unit vectors, that is the following relationships exist $[{\boldsymbol{\mathrm{\epsilon}}}\times{\boldsymbol{\mathrm{\pi}}}]={\boldsymbol{\mathrm{\kappa}}}+{\cal O}(r_{g}^{2}),$ $[{\boldsymbol{\mathrm{\pi}}}\times{\boldsymbol{\mathrm{\kappa}}}]={\boldsymbol{\mathrm{\epsilon}}}+{\cal O}(r_{g}^{2}),$ $[{\boldsymbol{\mathrm{\kappa}}}\times{\boldsymbol{\mathrm{\epsilon}}}]={\boldsymbol{\mathrm{\pi}}}+{\cal O}(r_{g}^{2})$, thus, $({\boldsymbol{\mathrm{\epsilon}}}\cdot{\boldsymbol{\mathrm{\pi}}})=({\boldsymbol{\mathrm{\epsilon}}}\cdot{\boldsymbol{\mathrm{\kappa}}})=({\boldsymbol{\mathrm{\pi}}}\cdot{\boldsymbol{\mathrm{\kappa}}})=0+{\cal O}(r_{g}^{2})$. One can also verify that ${\boldsymbol{\mathrm{\epsilon}}}^{2}={\boldsymbol{\mathrm{\pi}}}^{2}={\boldsymbol{\mathrm{\kappa}}}^{2}=1+{\cal O}(r_{g}^{2}).$

B.4 Spherical waves in the weak and static gravity

We know from quantum mechanics that spherical waves are important for the scattering problem. To study spherical waves in a weak and static gravitational field, we need to find solutions to the EM field by solving (18), namely:

[TABLE]

We seek a spherically symmetric solution with the following properties:

[TABLE]

In this case the d’Alembertian $\Delta\psi$ reduces to

[TABLE]

Therefore, (205) takes the form

[TABLE]

A formal solution to (208) may be given in the terms of confluent hypergeometric function Abramowitz and Stegaun (1965):

[TABLE]

where ${}_{1}F_{1}$ is the confluent hypergeometric function of the first kind (213) and $A$ is arbitrary constant.

Following the same approach that was demonstrated in Sec. II.3, we studied the asymptotic behavior of the solution (209). It turned out that such a solution may be given as follows:

[TABLE]

By choosing the constant $A=\pm ike^{\frac{\pi}{2}kr_{g}}{\Gamma(1\pm ikr_{g})}$, we may present the solution for a spherical wave in the weak and static gravity in the following form:

[TABLE]

representing both incoming and outgoing radial waves, with $c_{1},c_{2}$ being arbitrary constants. Note that the spherical wave solution (211) that we obtained is consistent with the solution for the phase of a radially propagating beam of light (i.e., radial geodesic) given by (190). Equation (211) establishes the functional dependence of the logarithmic term, which is important for the discussions of the scattering problem in Sec. II.3.

Appendix C The confluent hypergeometric function

C.1 Mathematical properties of the confluent hypergeometric function

We present some of the properties of the confluent hypergeometric function, denoted here as $F[\alpha|\beta|w]$, which are useful to derive our results. As defined (e.g., Abramowitz and Stegaun (1965)), $F[\alpha|\beta|w]$ is the regular solution of

[TABLE]

It is given by Abramowitz and Stegaun (1965):

[TABLE]

The function ${}_{1}F_{1}$ satisfies the following identities:

[TABLE]

Specifically,

[TABLE]

In Sec. III.3 we introduced two useful functions (106):

[TABLE]

Equation (216) leads to the following useful relation between $F[1]$ and $F[2]$:

[TABLE]

We will use this property when evaluating various contributions to the EM field on the image plane and the relevant Poynting vector, discussed in Secs. III.7 and III.8, correspondingly.

C.2 Asymptotic behavior of $F[1]$ and $F[2]$ at large values of argument

The asymptotic form of ${}_{1}F_{1}$ for large $|w|$, fixed $\alpha$, $\beta$ can be obtained by writing Abramowitz and Stegaun (1965); Burke (2011):

[TABLE]

where the functions ${\tt W}_{1}$ and ${\tt W}_{2}$ have the following asymptotic behavior Messiah (1968); Deguchi and Watson (1986):

[TABLE]

with the function $\tt G$ given Landau and Lifshitz (1989) in the form

[TABLE]

where the integration path $C_{1}$ goes from minus infinity around the origin ($t=0$) counterclockwise and back to minus infinity. Integrating by parts, we obtain the asymptotic series

[TABLE]

This is an asymptotic expansion. For arbitrary values of $\alpha,\beta$ and $w$, successive terms may eventually grow in size beyond limit. However, it is true that there exist functions $|\vartheta_{n}(\alpha,\beta,w)|<1$ such that

[TABLE]

i.e., when the series is truncated after $(n-1)$ terms, the error is no greater than the $n$-th term Szász (1951); Carrier et al. (1966).

Given the asymptotic properties of ${}_{1}F_{1}[\alpha,\beta,w]$ from (219), we take the solution to equation (18) which is given by (20) as \psi(\boldsymbol{\mathrm{r}})=\psi_{0}e^{ikz}{}_{1}F_{1}\big{(}ikr_{g},1,ik(r-z)\big{)} and split it in the form of $\psi(\boldsymbol{\mathrm{r}})=\psi_{\tt inc}(\boldsymbol{\mathrm{r}})+\psi_{s}(\boldsymbol{\mathrm{r}})$, where $\psi(\boldsymbol{\mathrm{r}})_{\tt inc}$ is the incoming and $\psi_{\tt s}(\boldsymbol{\mathrm{r}})$ is the scattered waves, correspondingly, which are given as

[TABLE]

Using the asymptotic forms (220) and (221), for large values of the argument $k(r-z)\rightarrow\infty$, functions ${\tt W}_{1}$ and ${\tt W}_{2}$ have the following asymptotic behavior:

[TABLE]

From the asymptotic expansion of ${\tt G}$ given by (223), we find that

[TABLE]

where, in (230), we used the fact that for the large values of the argument $k(r-z)\rightarrow\infty$ and for the high-frequency EM waves $kr_{g}\gg 1$. This allows us to write

[TABLE]

where in (232) we neglected the term ${\cal O}(r_{g}^{2}/(r-z)^{2})$, as being beyond the first post-Newtonian approximation taken in (1)–(2). Also, examining the order terms in these approximations, we note that although their absolute magnitudes are large, they are small compared to the logarithmic term $kr_{g}\ln k(r-z)$ present in the series expansion of the preceding exponential. That is to say that the order term contributes to the Shapiro delay (which is present in the phase of (231) in the form of $\delta d_{\tt shap}=-r_{g}\ln k(r-z)+{\cal O}(r_{g}^{2})$) at the second post-Newtonian order, namely it is of the order of ${\cal O}(r_{g}^{2}/(r-z))$, which is beyond the first post-Newtonian approximation accepted in this paper.

Collecting the terms, we may now present the asymptotic behavior F[1]=F[ikr_{g},1,ik(r-z)\big{]} at large values of the argument $k(r-z)\rightarrow\infty$, which, to post-Newtonian order, is given as below:

[TABLE]

The approximations given by (231)–(232) and by the resulting expression (233) are good so long as $r_{g}/(r-z)\lesssim 1$, which, together with $z=r\cos\theta$, yields a constraint

[TABLE]

Similarly, we study the behavior of the function $F[2]=F[1+ikr_{g},2,ik(r-z)]$. First, we present $F[2]={\tt W}_{3}+{\tt W}_{4}$, where for large values of $k(r-z)$, functions ${\tt W}_{3}$ and ${\tt W}_{4}$ have the following asymptotic behavior:

[TABLE]

As a result, the asymptotic behavior of F[2]=F[1+ikr_{g},2,ik(r-z)\big{]} at large values of the argument $|w|=k(r-z)\rightarrow\infty$ and angles $\theta$ outside the immediate vicinity of the optical axis, i.e., satisfying (234), is given as

[TABLE]

C.3 Asymptotic behavior of $F[1]$ and $F[2]$ at small angles

To understand the properties of the SGL near the optical axis, we need to the investigate the behavior of the solution at small angles. Based on the properties of the hypergeometric function (213), here we consider the behavior of $F[1]$ and $F[2]$ from (217) when $\theta$ is small. Using $\alpha=ikr_{g}$ and $w=ikr(1-\cos\theta)$, we define

[TABLE]

We next rearrange (213) as

[TABLE]

where $()_{n}$ denotes Pochhammer’s symbol666http://mathworld.wolfram.com/PochhammerSymbol.html with $()_{0}=1$, and $s(n,k)$ is the Stirling number of the first kind Abramowitz and Stegaun (1965); $s(0,0)=1$. Reversing the order of summation yields

[TABLE]

with $A_{n}=(-1)^{n}\sum\limits_{k=0}^{\infty}s(n+k,k)\dfrac{(\alpha w)^{k}}{[(n+k)!]^{2}}$. The Stirling number of the first kind can be evaluated Abramowitz and Stegaun (1965) in terms of the Stirling number of the second kind, which, in turn, also has a closed form sum:

[TABLE]

This can be evaluated for specific values of $n$:

[TABLE]

We also note that the Bessel functions are given by $J_{n}(2\sqrt{x})=(\sqrt{x})^{n}\sum\limits_{k=0}^{\infty}\dfrac{(-x)^{k}}{k!(n+k)!}$. Given $x=-\alpha w$, we have $A_{n}=(-1)^{n}\sum\limits_{k=0}^{\infty}s(n+k,k)\dfrac{(-x)^{k}}{[(n+k)!]^{2}}$, therefore

[TABLE]

Substituting (242)–(245) into (240), we obtain a very useful expression for the confluent hypergeometric function $F[1]={}_{1}F_{1}[\alpha,1,w]$ in terms of Bessel functions:

[TABLE]

This result is also consistent with (13.3.8) in Abramowitz and Stegaun (1965). Using the properties of the Bessel functions Abramowitz and Stegaun (1965), namely that

[TABLE]

we can present $J_{4}(2\sqrt{x})$ in (246) as $J_{4}(2\sqrt{x})=({3}/{\sqrt{x}})J_{3}(2\sqrt{x})-J_{2}(2\sqrt{x})$, which allows us to write (246) in a slightly different form as

[TABLE]

Following the same approach, we may obtain a relation for the function $F[2]$:

[TABLE]

In a small angle approximation we use $\theta=\rho/z+{\cal O}((\rho/z)^{2})$ and noting that $r\sim z+{\cal O}(r_{g})$, we present $w$ as $|w|=kr(1-\cos\theta)\approx\pi\rho^{2}/z\lambda+{\cal O}(r_{g})$. Thus, in the immediate vicinity of the optical axis $|w|<1$. As a result, the functions $F[1]$ and $F[2]$ may be presented as

[TABLE]

When $|w|$ is small enough such that terms containing $w^{2}$ may also be omitted, we can keep only the leading terms in these expressions:

[TABLE]

Based on these expressions, we may compute the following combinations:

[TABLE]

where $A^{*}$ denotes a complex conjugate of $A$ and $x$ is given by (238).

Appendix D Properties of Coulomb functions

D.1 Differential equation

In spherical coordinates, the problem of scattering of an EM wave on a gravitational monopole for each value of partial momentum $\ell$ leads to the following radial equation (we follow very closely the discussion in Messiah (1968)):

[TABLE]

Partial solutions to this equation may be obtained in terms of spherical Coulomb functions. These are the functions of $\rho=kr$. They depend on the wavenumber, $k$, distance to the deflector, $r$, and its Schwarzschild radius $r_{g}$. There exists a regular solution ($\sim r^{\ell+1}$) at the coordinate origin, $F_{\ell}(kr_{g},kr)$ and irregular solutions $G_{\ell}(kr_{g},kr)$ together with $H^{+}_{\ell}(kr_{g},kr),H^{-}_{\ell}(kr_{g},kr)$ that are singular ($\sim 1/r^{\ell}$) at the coordinate origin.

With a substitution

[TABLE]

equation (256) may be presented in the form of the Laplace equation:

[TABLE]

where $\alpha=\ell+1-ikr_{g}$, $\beta=2\ell+2$ are complex constant coefficients. The solution to (258) is the confluent hypergeometric function ${}_{1}F_{1}$ given in (213). Equation (258) has a regular solution ${}_{1}F_{1}[\ell+1-ikr_{g},2\ell+2,z]$ and two irregular solutions ${\tt W}_{1}[\ell+1-ikr_{g},2\ell+2,z]$ and ${\tt W}_{2}[\ell+1-ikr_{g},2\ell+2,z]$. Based on these functions we can construct the solutions that we discuss below.

D.2 Relationships between the Coulomb functions and their asymptotic properties

Given ${}_{1}F_{1}[\ell+1-ikr_{g},2\ell+2,z]$ and ${\tt W}_{1,2}[\ell+1-ikr_{g},2\ell+2,z]$, we introduce the following functions (see Messiah (1968)):

[TABLE]

Alternatively, we can define a different, but equivalent, set of solutions (156) with $F_{\ell}(kr_{g},kr)$ given by (259), but also defining $G_{\ell}(kr_{g},kr)$ and the Coulomb–Hankel functions $H^{(\pm)}_{\ell}(kr_{g},kr)$ as

[TABLE]

where $U(\alpha,\beta,z)$ is the corresponding irregular confluent hypergeometric function defined in Abramowitz and Stegaun (1965).

The quantities $c_{\ell}$ and $\sigma_{\ell}$ (i.e., Coulomb phase shift) are the following functions of $r_{g}$:

[TABLE]

For $\ell=0$, (263) takes the form

[TABLE]

or, for $\ell\not=0$, (263) takes the form

[TABLE]

Both $F_{\ell}$ and $G_{\ell}$ are real-valued functions:

[TABLE]

The asymptotic behavior of the Coulomb functions outside the turning point defined by (87), when $r\rightarrow\infty$ and $r\gg r_{\tt t}=-r_{g}+\sqrt{r^{2}_{g}+\ell(\ell+1)/k^{2}}$, is well known and given as

[TABLE]

Their behavior near the origin of the coordinate system, when $r\rightarrow 0$, is

[TABLE]

In the case when $r_{g}=0$, then up to a factor of $kr$ one obtains spherical Bessel functions $j_{\ell},n_{\ell},h^{(\pm)}_{\ell}$:

[TABLE]

where $j_{\ell},n_{\ell},h^{(\pm)}_{\ell}$ are

[TABLE]

Appendix E Representation of the field in terms of Debye potentials

To represent the EM field equations in terms of Debye potentials, we start with (7)–(8). Assuming, as usual (we follow closely the discussion presented in Born and Wolf (October 13, 1999), adapted for the gravitational lens), the time dependence of the field in the form $\exp(-i\omega t)$ where $k=\omega/c$, the time-independent parts of the electric and magnetic vectors satisfy Maxwell’s equations: Eq. (7)–(8) in their time-free form:

[TABLE]

In spherical coordinates (Fig. 2), the field equations (278)–(279) with the help of (154)-(155) to order ${\cal O}(G^{2})$ become

[TABLE]

Our goal is to find a general solution to these equations in the form of a superposition of two linearly independent solutions \big{(}{}^{e}{\boldsymbol{\mathrm{D}}},{}^{e}{\boldsymbol{\mathrm{B}}}\big{)} and \big{(}{}^{m}{\boldsymbol{\mathrm{D}}},{}^{m}{\boldsymbol{\mathrm{B}}}\big{)} that satisfy the following relationships:

[TABLE]

With $\hat{B}_{r}={}^{e}{\hskip-2.0pt}\hat{B}_{r}=0$, (281) and (282) become

[TABLE]

Substituting these relationships into (284) and (285) we obtain

[TABLE]

From ${\rm div}(u^{2}{}^{e}{\boldsymbol{\mathrm{B}}})=0$ given by Eq. (8) and using our assumption that $\,{}^{e}{\hskip-2.0pt}\hat{B}_{r}=0$ we have

[TABLE]

which ensures that (283) is also satisfied, since it becomes, after the substitution from (288), (289),

[TABLE]

which is satisfied because of (292). The complementary case with ${}^{m}{\hskip-2.0pt}\hat{D}_{r}=0$ is treated identically, in accord with (287).

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 electric wave). These can both be derived from the corresponding Debye scalar potentials ${}^{e}{\hskip-1.0pt}\Pi$ and ${}^{m}{\hskip-1.0pt}\Pi$.

Given ${}^{e}{\hskip-2.0pt}\hat{B}_{r}=0$, ${}^{e}{\hskip-2.0pt}\hat{D}_{\phi}$ and ${}^{e}{\hskip-2.0pt}\hat{D}_{\theta}$ in (283) can be represented as a scalar field’s gradient:

[TABLE]

Using

[TABLE]

in (294), we obtain

[TABLE]

It can be seen that (288) and (289) are satisfied by

[TABLE]

If we substitute both of (297) into (280) we obtain

[TABLE]

Substituting expressions (297) into (290)–(291) yields \dfrac{-ik}{\sin\theta}\dfrac{\partial}{\partial\phi}\Big{\{}\dfrac{\partial}{\partial r}\Big{[}\dfrac{1}{u^{2}}\dfrac{\partial}{\partial r}\big{(}r\,{}^{e}{\hskip-1.0pt}\Pi\big{)}\Big{]}+k^{2}u^{2}(r\,{}^{e}{\hskip-1.0pt}\Pi)-{}^{e}{\hskip-2.0pt}\hat{D}_{r}\Big{\}}=0 and ik\dfrac{\partial}{\partial\theta}\Big{\{}\dfrac{\partial}{\partial r}\Big{[}\dfrac{1}{u^{2}}\dfrac{\partial}{\partial r}\big{(}r\,{}^{e}{\hskip-1.0pt}\Pi\big{)}\Big{]}+k^{2}u^{2}(r{}^{e}{\hskip-1.0pt}\Pi)-{}^{e}{\hskip-2.0pt}\hat{D}_{r}\Big{\}}=0, i.e., the derivative of the same expression with respect to both $\phi$ and $\theta$ vanishes. This is clearly satisfied if we set the expression itself to 0. Dividing by $u^{2}$ and using (298) leads to

[TABLE]

Defining $u^{\prime}=\partial u/\partial r$, this equation may be rewritten as

[TABLE]

which is the wave equation for the quantity ${\,{}^{e}{\hskip-1.0pt}\Pi}/{u}$:

[TABLE]

We are concerned only with the field produced by the gravitational monopole, thus the quantity $u$ has the from $u({\boldsymbol{\mathrm{r}}})=1+{r_{g}}/{2r}+{\cal O}(r^{-3},c^{-4}),$ as given by (171). With this, we can rewrite (301) as

[TABLE]

Equation (302) is similar to the Schrödinger equation of quantum mechanics, used to describe scattering on the Coulomb potential. However, this equation has an extra potential of $r_{g}/r^{3}$. It is known Messiah (1968) (and also shown in Appendix F) that the presence of potentials of $\propto 1/r^{3}$ in (302) does not alter the asymptotic behavior of the solutions. Reference Matzner (1968), discusses justification for neglecting the $r^{-3}$ term in (302), which reduces this equation to the time-independent Schrödinger equation that describes scattering in a Coulomb potential:

[TABLE]

In the case of the SGL, we will always be at the distances which are much larger than the Sun’s Schwarzschild radius, thus, we may neglect the term ${r_{g}}/{r^{3}}$ in (302). We will use (303) for the purposes of establishing the properties of the EM wave diffraction by the solar gravitational lens. An identical equation may be obtained for ${}^{m}{\hskip-1.0pt}\Pi$.

By means of (299), Eq. (298) may be written as

[TABLE]

It can be verified by substituting (296)–(299) and (304) into (280)–(285) that we have obtained a solution of our set of equations. In a similar way we may consider the magnetic wave. We find that this wave can be derived from a potential ${}^{m}{\hskip-1.0pt}\Pi$ which satisfies the same differential equation (299) as ${}^{e}{\hskip-1.0pt}\Pi$.

The complete solution of the EM field equations is obtained by adding the two fields (as discussed in Mie (1908); Born and Wolf (October 13, 1999); Kerker (1969)), namely ${\boldsymbol{\mathrm{D}}}={}^{e}{\hskip-1.0pt}{\boldsymbol{\mathrm{D}}}+{}^{m}{\hskip-1.0pt}{\boldsymbol{\mathrm{D}}}$; and ${\boldsymbol{\mathrm{B}}}={}^{e}{\hskip-1.0pt}{\boldsymbol{\mathrm{B}}}+{}^{m}{\hskip-1.0pt}{\boldsymbol{\mathrm{B}}}$; this gives

[TABLE]

Both potentials $\,{}^{e}{\hskip-1.0pt}\Pi$ and $\,{}^{m}{\hskip-1.0pt}\Pi$ are solutions of the differential equation (301), which, in the case of the weak gravity characteristic for the SGL, is given by (303).

Appendix F Solution for the radial equation in the WKB approximation

We consider Eq. (302) for the Debye potentials. Using the representation given by (62) and remembering $\alpha=\ell(\ell+1)$, we obtain the following equation for the radial function $R$:

[TABLE]

Following an approach similar to that presented in Herlt and Stephani (1976), we explore an approximate solution to (311) using the methods of stationary phase (i.e., the Wentzel–Kramers–Brillouin, or WKB approximation). As we are interested in the case when $k$ is rather large (for optical wavelengths $k=2\pi/\lambda=6.28\cdot 10^{6}\,{\rm m}^{-1}$), we will be looking for an asymptotic solution as $k\rightarrow\infty$. In fact, we will be looking for a solution in the form of

[TABLE]

Technically, however, it is more convenient to search for a solution to (311) in an exponential form:

[TABLE]

Defining $R^{\prime}=dR/dr$, with the help of a substitution of $R^{\prime}/R=w$, for the function $w$ we obtain the following equation:

[TABLE]

Using this substitution we have

[TABLE]

Substituting (315) into (314) we obtain

[TABLE]

Now, if we equate the terms with respect to the same powers of $k$, we get

[TABLE]

These equations may be solved as

[TABLE]

Using this approach we can identify $\alpha_{1}(\rho),\alpha_{2}(\rho)$, … Substituting (318) into (313), we have

[TABLE]

As a result, we obtain two approximate solutions for the partial radial function $R_{\ell}$:

[TABLE]

where $c_{\ell}$ and $d_{\ell}$ account for all the constants relevant to the point $r_{0}$ in (319)–(321).

We note that (311) is similar to the radial solution of the Schrödinger equation of quantum mechanics, which is used to describe scattering in a Coulomb potential. However, this equation has an extra potential in the form of $r_{g}/r^{3}$. It is known Messiah (1968) that the presence in (311) of potentials of $\sim 1/r^{3}$ does not alter the asymptotic behavior of the solutions. One can verify that neglecting $r_{g}/r^{3}$ in (311) leads to the following form of the radial function $R_{\ell}$:

[TABLE]

We see that the omission of the ${r_{g}}/{r^{3}}$ term in (311) leads to appearance of an “uncompensated” term ${r_{g}}/{4kr^{2}}=(1/8\pi)({r_{g}\lambda}/{r^{2}})$ in the exponent of (323). This term is extremely small; it decays fast as $r$ increases, and, thus, it may be neglected in the solution for the radial function for any practical purpose. A similar point was made in Matzner (1968), suggesting that one can neglect the $r^{-3}$ terms in (311) and reduce the problem to the case of the Schrödinger equation describing scattering in a Coulomb potential.

As a result, to describe the scattering of a plane EM wave by a gravitational monopole, we neglect the term ${r_{g}}/{r^{3}}$ in (311) and approximate it such that it becomes

[TABLE]

Finally, we may further improve the asymptotic expression for $R_{\ell}$ from (322) by accounting for the Coulomb phase shifts as given in (269)–(272). This can be done by simply redefining the constants $c_{\ell}$ and $d_{\ell}$ as

[TABLE]

This leads to the following expression for the radial function $R_{\ell}$:

[TABLE]

where the term $r_{g}/2r$ in the phase was neglected. As the asymptotic behavior of the Coulomb functions (269)–(272) was obtained for very larger distances from the turning point (87), or for $r\gg r_{\tt t}$, the solution (326) improves them by extending the argument of the Coulomb functions to shorter distances, closer to the turning point. (A similar result was obtained in Nambu and Noda (2016) using different approach developed to study image formation in gravitational lensing Nambu (2013).)