Integral Field Spectroscopy with the Solar Gravitational Lens

There are multiple, equivalent mathematical formulations for unbound hyperbolic trajectories. The formulas that follow are typically found in textbooks on astrodynamics (Schaub & Junkins 2003; McClain & Vallado 2001). The recognizable Cartesian parameterization

$$\begin{eqnarray}&&\displaystyle \frac{{x}^{2}}{{a}^{2}}-\displaystyle \frac{{y}^{2}}{{b}^{2}}=1\end{eqnarray} \tag{ 1 }$$

describes the shape of the hyperbola with two numbers (a, b) named the semimajor and semiminor axes. The graph has two asymptotes at $y=\pm \tfrac{b}{a}x$ and two branches opening along the x-axis that approach the asymptotes as the coordinates approach infinity.

A different formulation uses polar coordinates (r, θ)

$$\begin{eqnarray}&&r=\displaystyle \frac{l}{1-e\cos (\theta )}\end{eqnarray} \tag{ 2 }$$

where the shape is described by two new numbers, the semi-latus rectum $l=-\tfrac{{b}^{2}}{a}$ and eccentricity $e=\sqrt{1+{\left(\tfrac{b}{a}\right)}^{2}}$. An important note is that the polar coordinates have their origin as the focus of the right branch, and so to match the Cartesian grid, the (x, y) coordinates can be recovered with the relations $x=r\cos (\theta )-{ae}$ and $y=r\sin (\theta )$.

However, for the purposes of this investigation, we wish to parameterize the hyperbola with two different numbers, the radius at periastron r_p and the velocity at periastron v_p , so that one may think of the hyperbola shape phase space as a set of initial conditions for a rocket launching from an inner solar system before reaching the SGL. Everywhere along the hyperbola, the velocity associated with the trajectory is governed by the vis-viva equation

$$\begin{eqnarray}&&{v}^{2}=\mu \left(\displaystyle \frac{2}{r}-\displaystyle \frac{1}{a}\right)\end{eqnarray} \tag{ 3 }$$

where μ = GM_⊙ is the standard gravitational parameter of the Sun. So simply plugging in r_p , v_p , and solving for a allows determination of the semimajor axis

$$\begin{eqnarray}&&a=\displaystyle \frac{\mu {r}_{p}}{2\mu -{v}_{p}^{2}{r}_{p}}.\end{eqnarray} \tag{ 4 }$$

Additionally, the specific angular momentum $h={r}_{p}{v}_{p}=\sqrt{\mu l}$ of the orbit is conserved and can be used to find the semiminor axis

$$\begin{eqnarray}&&b=\displaystyle \frac{{r}_{p}^{3/2}{v}_{p}}{\sqrt{{v}_{p}^{2}{r}_{p}-2\mu }}.\end{eqnarray} \tag{ 5 }$$

However, this equation will return an imaginary number when $2\mu \gt {v}_{p}^{2}{r}_{p}$, in this case, the chosen r_p , v_p do not correspond to a hyperbola. An example plot of a basic hyperbola is shown in Figure 1.

Zoom In Zoom Out Reset image size

Using the (r_p , v_p ) basis for hyperbola shapes, we aim to answer a few basic questions about the properties of these orbits to inform a design of a craft that intends to visit the SGL. The questions are: "How long will it take to reach a heliocentric distance of r_SGL ≈ 600 au?" "What is the angular velocity of the craft with respect to the Sun at that point in time?" and "How much acceleration is needed to cancel the relative lateral motion to enter a radial hyperbolic trajectory?"

The first question can be answered using the relationship between mean anomaly M, eccentric (or hyperbolic) anomaly E, and time t. Specifically,

$$\begin{eqnarray}&&E=\mathrm{arcosh}\displaystyle \frac{\cos \theta +e}{1+e\cos \theta }\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&M=e\sinh E-E\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&t=\sqrt{\displaystyle \frac{-{a}^{3}}{\mu }}M.\end{eqnarray} \tag{ 8 }$$

Then it is just matter of finding θ corresponding to r_SGL using

$$\begin{eqnarray}&&{\theta }_{\mathrm{SGL}}=\arccos \displaystyle \frac{l/{r}_{\mathrm{SGL}}-1}{e}\end{eqnarray} \tag{ 9 }$$

to compute the time a craft would need to reach such an extreme heliocentric distance. Contours are plotted for a variety of cases in Figure 2(a).

Zoom In Zoom Out Reset image size

The second question can be analytically derived using the conservation of specific angular momentum. Since $h=r\times \left(r\tfrac{d\theta }{{dt}}\right)$,

$$\begin{eqnarray}&&\displaystyle \frac{d\theta }{{dt}}=\displaystyle \frac{h}{{r}_{\mathrm{SGL}}^{2}}.\end{eqnarray} \tag{ 10 }$$

which is equivalent to Kepler's second law. We convert these angular velocities into planet-crossing and system-crossing timescales ${T}_{\mathrm{planet}}={\theta }_{\mathrm{planet}}/\tfrac{d\theta }{{dt}}$ and ${T}_{\mathrm{system}}={\theta }_{\mathrm{system}}/\tfrac{d\theta }{{dt}}$ using θ_planet = 2R_Earth/d_s , and θ_system = 2a_sys/d_s where d_s = 40.54 ly is the system distance, a_sys = 0.062 au is the target planet semimajor axis, and R_Earth is the Earth's radii. These estimates correspond to the TRAPPIST-1 system (Gillon et al. 2017) and are plotted in Figures 2(b) and (c).

The last question can be tackled by means of computing the flight path angle ϕ_fpa, which gives the angle between the velocity vector and the direction perpendicular to the radial vector. The flight path angle is computed with

$$\begin{eqnarray}&&{\phi }_{\mathrm{fpa}}=\arctan \displaystyle \frac{e\sin {\theta }_{\mathrm{SGL}}}{1+e\cos {\theta }_{\mathrm{SGL}}},\end{eqnarray} \tag{ 11 }$$

and the magnitude of the velocity vector ${v}_{\mathrm{SGL}}\,=\sqrt{\mu (2/{r}_{\mathrm{SGL}}-1/a)}$ is computed with the vis-viva equation. Thus the delta-v required to eliminate the transverse motion and enter a radial hyperbolic orbit is

$$\begin{eqnarray}&&{dv}={v}_{\mathrm{SGL}}\sin \left(\displaystyle \frac{\pi }{2}-{\phi }_{\mathrm{fpa}}\right).\end{eqnarray} \tag{ 12 }$$

Values of which are plotted in Figure 2(d).

Examining the plots in Figure 2, a few important trade-offs in the orbit parameter space become apparent. Orbits with large velocity at periastron correspondingly have the shortest time to reach great heliocentric distance, while the radius at periastron has a small effect, with a slight preference for starting farther out of the solar system. However, these wide-starting orbits have much greater angular velocity at r_SGL, and correspondingly smaller planet-crossing and system-crossing timescales. The orbits with the smallest angular velocity (and thus largest crossing timescales, as well as smallest delta-v to enter the radial trajectory) correspond to orbits with small radius at periastron and large velocity at periastron.

To make a concrete example, suppose a craft is launched in a similar manner to the New Horizons spacecraft, which was launched in 2006 to investigate Pluto (Guo & Farquhar 2008). Its trajectory departed Earth with a velocity estimate of 16.2 km s⁻¹, which was the highest recorded departure speed at the time, which corresponds to a heliocentric velocity of 46 km s⁻¹. Yet New Horizons did not fly by Pluto at 33 au until nearly a decade later (Stern et al. 2015). If a similar rocket would be used to launch a mission to the SGL, our calculations indicate it would take ∼130 yr to reach a distance of 600 au, only to cross the image of the target planet in ∼40 s, or the entire system in just ∼1000 minutes. However, the craft could correct the lateral motion with a delta-v of only ∼80 m s⁻¹ to enter a pure radial trajectory and extend its potential observation time further. Additionally, New Horizons used a Jupiter fly-by maneuver to obtain a velocity boost of ∼4 km s⁻¹, which can shorten the long wait but also causes the craft to have a wider heliocentric trajectory with a larger angular velocity, which would cost additional delta-v in the end.

Regardless of how elegant the analytic solutions are, they are only a two-body approximation to the dynamics of a craft exiting the solar system. The gravitational dynamics of all of the bodies in the solar system are known to cause motion of the Sun around the barycenter (Cionco & Pavlov 2018), due to the gravitational interactions of all of the massive bodies in orbit. The N-body problem has been known since Newton, and is essentially intractable except for special cases (Musielak & Quarles 2014) or statistical ensembles (Stone & Leigh 2019). However, numerical techniques can be quite effective to solve the equations of motion, and so we use REBOUND (Rein & Liu 2012) to investigate these effects. REBOUND has implemented a high-order adaptive time step integrator called IAS15 (Rein & Spiegel 2015), which enables machine-precision accuracy of complex N-body orbits at extremely cheap computational cost, and can be used to simulate and compute the motion of the Sun around the solar system barycenter extremely accurately into the future.

To evaluate the magnitude of this effect, we initialize a REBOUND simulation with N_particles = 10 particles, one for the Sun, each of the eight planets, and Pluto, using the built-in tool to load barycentric orbit elements from the JPL HORIZONS online ephemeris database (Giorgini et al. 1996). Then REBOUND integrates the equations of motion for the N-body system of particles forward in time for 100 yr, using N_timesteps = 10⁴ samples. This results in a sampling rate of two per week, enough to resolve effects from Mercury's short 88 day orbit. The results of the simulation consist of a data cube of size (N_particles, N_timesteps, 6), which contains the six x-, y-, and z- positions and velocities of each of the particles in the integration at every time step.

To evaluate the effect of the barycentric motion of the Sun on the pointing of the SGL, we investigate the worst-case scenario. Since the majority of the motion of the Sun occurs in the invariable plane (Souami & Souchay 2012), defined by the point of the barycenter and perpendicular to the total angular momentum, a target planet that is pointed directly along this angular momentum vector will suffer from the largest perturbations of the position of the Sun along the line of sight. Any other target pointing will be somewhat better from this perspective as the motion of the Sun is projected into the plane of observation. In our coordinate system, this plane is defined by the x- and y-coordinates, and so the telescope craft visiting the SGL has a position of (0, 0, r_SGL) once it has reached the observation location and maneuvered to enter a purely radial orbit trajectory. This next part of the analysis ignores the additional z-axis motion for simplicity. Although it is a small effect, this could be corrected by adding an additional term for a time-dependent z-coordinate, $z={r}_{\mathrm{SGL}}\,+(t/{T}_{\mathrm{SGL}})\sqrt{\mu (2/z-1/a)}$, which would require solving the vis-viva equation for v(z(t)).

In Figure 3(a), the (x, y) position of the Sun relative to the barycenter is plotted in units of solar radii, showing the reflex motion due to the perturbations from the planets causing a displacement of around ∼1–2 R_⊙. These values are used to estimate the angular velocity of the "SGL pointing" vector defined as starting at the stationary craft at (0, 0, r_SGL) and directed at the position of the Sun (x_⊙, y_⊙, 0). The discretely sampled position vector is the output of the REBOUND simulation. First, we compute the discrete x- and y- second-order position differences dx = (x_t+1 − x_t−1)/2 and dy =(y_t+1 − y_t−1)/2, using np.gradient, which handles the edge cases using a first-order finite difference, as well as the discrete time difference dt = 100 yr/10⁴ samples. The error on these estimates of the derivative is on the order of ${ \mathcal O }({{dt}}^{2})$. Finally, the angular velocity of the pointing vector is computed with

$$\begin{eqnarray}&&\displaystyle \frac{d\phi }{{dt}}=\displaystyle \frac{\sqrt{{{dx}}^{2}+{{dy}}^{2}}}{{r}_{\mathrm{SGL}}{dt}},\end{eqnarray} \tag{ 13 }$$

which can be similarly converted into a planet-crossing timescale as before, using ${T}_{\mathrm{cross}}={\theta }_{\mathrm{planet}}/\tfrac{d\phi }{{dt}}$. Both the angular velocity and planet-crossing timescale are plotted over time in Figures 3(b) and (c). From these estimates, one can see the average angular velocity is around 2–3 mas per day, and that the corresponding planet-crossing timescales are between 200–300 s. So to design a craft that obtains longer exposure times at ultra-precise alignments, additional navigation acceleration will be needed to counteract this effect.

Zoom In Zoom Out Reset image size

In order to compute the amount of acceleration necessary to counteract the barycentric motion of the Sun, we first compute the location the craft must station-keep at to maintain the SGL aimed at a fixed point on the sky. If d_s is the distance between the Sun and source or target planet, and d_l = r_SGL is the distance between the craft and the Sun or lens, a simple similar triangle argument shows the position (x, y, z = d_l ) required to station-keep is enhanced by a small factor

$$\begin{eqnarray}&&x={x}_{\odot }\Space{0ex}{0.25ex}{0ex}(1+\displaystyle \frac{{d}_{l}}{{d}_{s}}\Space{0ex}{0.25ex}{0ex})\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&y={y}_{\odot }\Space{0ex}{0.25ex}{0ex}(1+\displaystyle \frac{{d}_{l}}{{d}_{s}}\Space{0ex}{0.25ex}{0ex})\end{eqnarray} \tag{ 15 }$$

relative to the position of the Sun due to the additional distance. This enhancement is a small number that depends on the distance to the target being considered. For d_s = 40.54 ly and d_l = 600 au, $\tfrac{{d}_{l}}{{d}_{s}}\approx 2\times {10}^{-4}$. This position is shown as a function of time in Figure 4(a) in units of solar radii. Subsequently, the velocity ${v}_{x,y}=\tfrac{{dx},y}{{dt}}$ and acceleration ${a}_{x,y}=\tfrac{{{dv}}_{x,y}}{{dt}}$ required to maintain this position are computed numerically using the discrete second-order difference using np.gradient, and are plotted in Figures 4(b) and (c).

Zoom In Zoom Out Reset image size

Examining the figure, it is clear that a small, continuous acceleration of ∼2 × 10⁻⁷ m s⁻² in an oscillatory fashion is necessary to counteract the lateral motion of the Sun around the barycenter to maintain a fixed pointing of the SGL. Considering the possibilities of rocket designs, it seems that ion thrusters may be ideally suited to this task. While the specifics depend on a number of design choices, such a beam divergence, mass of the ion, and the electrical current and voltage used to accelerate the ions (Goebel & Katz 2008), a reasonable estimate of thrust from such an ion thruster is around ∼100 mN, which could give accelerations of 10⁻⁴ m s⁻² for a 1000 kg craft over very long durations, due to the high specific impulse of the ions, which is more than two orders of magnitude needed.

An alternative visualization of this acceleration trajectory is visualized in Figure 5. By converting to polar coordinates, the direction of the acceleration vector in the xy-plane is plotted as a function of time, which is now the radial coordinate. The log of the magnitude of the acceleration vector is represented by the color of the curve. One can see that acceleration curve spirals outward, as the vector rotates around with an approximate 10 yr period dominated by the effect of Jupiter. However, the perturbations of the other planets contribute as well, with consistent changes in the acceleration vector's direction and corresponding changes in magnitude.

Zoom In Zoom Out Reset image size

To continue our analysis, we estimate the total delta-v needed to maintain this acceleration over time as the integral ${{dv}}_{\mathrm{total}}(t)={\int }_{0}^{t}| {a}_{\mathrm{total}}(t^{\prime} )| {dt}^{\prime}$ of the magnitude of the acceleration vector $| {a}_{\mathrm{total}}(t)| =\sqrt{{a}_{x}{\left(t\right)}^{2}+{a}_{y}{\left(t\right)}^{2}}$. This total delta-v is a linear function of the length of the time you would like to the maintain the alignment for, with a slope corresponding to the average value $\langle | {a}_{\mathrm{total}}| \rangle =\tfrac{1}{t}{\int }_{0}^{t}| {a}_{\mathrm{total}}(t^{\prime} )| {dt}^{\prime}$. As a result, the approximate cost scaling to maintain alignment is

$$\begin{eqnarray}&&{{dv}}_{\mathrm{total}}(t)=\langle | {a}_{\mathrm{total}}| \rangle t\approx 6.7\displaystyle \frac{{\rm{m}}}{{\rm{s}}\cdot \mathrm{year}}t,\end{eqnarray} \tag{ 16 }$$

which demonstrates that the station-keeping cost to counteract the solar barycentric motion is sufficiently small as to be reasonably possible.

However, this entire analysis has so far ignored the motion of the target across the sky, composed of the relative proper motion of the system and the orbital motion of the target around its host star. A simple estimate of the acceleration that would be needed to counteract the orbital motion of a hypothetical planet with a face-on, circular orbit of semimajor axis a and period T is

$$\begin{eqnarray}&&{a}_{p}=\displaystyle \frac{{d}_{l}}{{d}_{s}}\displaystyle \frac{4{\pi }^{2}a}{{T}^{2}},\end{eqnarray} \tag{ 17 }$$

which can dominate over or be comparable to the effect of the solar barycentric motion depending on the assumed planet properties. If a = 1 au, T = 1 yr, d_l = 600 au, and d_s = 10 pc, then a_p ∼ 6 × 10^{− 6} m s⁻² ∼ 60 m s⁻¹ yr⁻¹, but if instead d_s = 100 pc, then a_p ∼ 0.2 × 10⁻⁶ m s⁻² ∼6 m s⁻¹ yr⁻¹. A more complete three-dimensional analysis such as that in Turyshev & Toth (2022) is necessary to evaluate the combined effect, as the orientation of the orbits is important to determine if this acceleration exacerbates or counteracts the solar barycentric motion. Ultimately, precise determination of the target ephemeris will be necessary to locate the exact position needed to send the craft, and predicting this location in advance may be rather difficult.

Different exoplanetary detection techniques are capable of constraining different planetary orbital elements. Transit photometry is able to determine the planet's period and to some extent its orbital inclination, as it must be nearly edge-on in order for a transit to occur (Deeg & Alonso 2018). Doppler spectroscopy or the radial velocity method measures the period and velocity semiamplitude, which can be used infer the orbit's eccentricity, mass-inclination product $m\sin i$ and, if the star's mass is known through another technique, the orbital semimajor axis (Eggenberger & Udry 2010). If the planet is amenable to both transit and radial velocity, the combination of both techniques can further constrain the orbit. The Rossiter-Mclaughlin effect can be used to measure the relative azimuthal angle on the sky between the planet's orbit and the stellar spin axis (Ohta et al. 2005; Masuda 2018), using asymmetric distortions in the line profiles of the stellar spectrum as the planet crosses and sequentially obscures different subsets of the stellar disk. However, absolute determination of the planet's orbital position angle on the sky would still require determination of the orientation of the stellar spin axis, which may be possible using spectro-interferometry to measure the astrometric motion associated with stellar rotation (Kraus et al. 2020), with spectro-astrometric analysis of high-resolution long-slit spectra (Lesage & Wiedemann 2014), or by combining projected spectroscopic rotation rates with photometric rotation periods (Bryan et al. 2020). It may also be possible to make these measurements during the early part of the mission to the SGL, by visiting the image of the star first to aid in navigation before moving to the image of the target planet. Turyshev & Toth (2022) discusses a strategy that uses astrometric observations of the host star using the SGL to help constrain the planet's position.

However, direct imaging of exoplanets can constrain planetary orbital elements by directly measuring the planet's astrometric position on the sky. Typical astrometric accuracies obtainable with current instruments on single telescopes are between 3 and 22 mas (Konopacky et al. 2016), while the most precise measurements currently obtained could be around ∼100 μas (GRAVITY Collaboration et al. 2019) using multiple telescopes combined into an interferometer. These uncertainties in measured precision are typically relative measurements, comparing the planet's position to the host star, and an absolute astrometric calibration would require precise determination of the stellar absolute astrometry, which could be measured using GAIA (Brown 2021; Lindegren et al. 2018). The uncertainty in angular position δ θ on the sky directly translates into a spatial, lateral uncertainty δ x = d_l δ θ of the image of the target in the SGL image plane. The large enhancement due to the great heliocentric distance means that accurate determination of the angular position on the sky of the target is critical to locating the position of the image of the target in the SGL focal plane. An angular uncertainty of δ θ ∼ 20 mas corresponds to an uncertainty in location of δ x ∼ 8800 km, while δ θ ∼ 100 μas corresponds to δ x ∼ 44 km. Since the image of the target planet is rescaled from its original size ${R}_{\mathrm{image}}=\tfrac{{d}_{l}}{{d}_{s}}{R}_{\mathrm{Earth}}$ by the "plate scale" of the SGL geometry, an Earth-sized planet at a distance of 40 ly has an image only ∼3 km in diameter. With current levels of uncertainty in relative astrometry, it could be missed even with navigation relative to the image of host star. A further caveat is that those measurements are typically for Jupiter-sized planets, since no Earth-like planets have been directly imaged, although mission concepts like HabEx and LUVOIR have been proposed to potentially achieve this (Gaudi et al. 2018; LUVOIR Collaboration 2019).

Additionally, the planet and star are continually in motion, due to the relative orbital motion of the system, and the proper motion and acceleration of the entire system within the galaxy. Typically stars have galactic velocities around ∼100 km s⁻¹ and accelerations around ∼10⁻¹⁰ m s⁻² (Silverwood & Easther 2019), which means determining the location needed to visit the image of either the star or planet would need to be predicted in advance, to plan an appropriate trajectory. Given that realistic timescales between launch and arrival at the SGL are of the order of ∼100 yr, and that typical proper motions are of the order of ∼.1'' yr⁻¹, the location of the image of a target star now and later would be displaced by millions of kilometers, with uncertainties in current proper motion further worsening the issue. Additionally, predicting the future orbital motion of the relative motion of the planet and star is subject to additional uncertainty.

Software such as orbitize! (Blunt et al. 2020) can be used to fit astrometric observations with three-dimensional Keplerian orbits, allowing estimates of position into the future, although the accuracy of prediction depends strongly on the orbital coverage of the data and its astrometric accuracy. To demonstrate this, Figure 6 was created. First, a series of simulated astrometric data are generated that correspond to a ground truth Keplerian orbit with a semimajor axis of 1 au, an eccentricity of 0.05, an inclination of 30°, for a system 100 pc away. Either 10 or 33 samples were generated over 0.3 or 1.0 yr with an observational uncertainty around ∼1 mas, which correspond to an orbital coverage of 33% and 100% for the partially and fully constrained cases respectively.

Zoom In Zoom Out Reset image size

The simulated data points are fit with the affine-invariant Markov Chain Monte Carlo algorithm developed by Foreman-Mackey et al. (2013), which samples the posterior distributions of all orbital parameters. After a burn-in time of 10⁴ steps and a total number of 10⁶ samples, the posterior distributions have mostly converged, although with less precision on the true orbital elements for the partially constrained orbit. For the partially constrained orbit, attempting to predict the location at a time corresponding to unconstrained and unmeasured locations in the orbit causes the posterior distribution of predicted locations to become large. Even with observational uncertainties on the order of ∼1 mas, the future location posterior has uncertainty an order of magnitude larger around ∼10 mas. However, a fully constrained orbit predicts the future location with accuracy comparable to the observational uncertainty on the data points themselves, which is a necessary precursor for navigation to the precise location of the planet's image at the SGL. This demonstrates the importance of complete astrometric orbital coverage on constraining the planetary orbital elements and predicting its future location.

The simplest model for the gravitational potential of the Sun is that of a point mass, ${\rm{\Phi }}=\tfrac{-{GM}}{r}$. This model has a very simple formula (Narayan & Bartelmann 1996)

$$\begin{eqnarray}&&{{\boldsymbol{\theta }}}_{{\boldsymbol{g}}}=\displaystyle \frac{2{R}_{S}}{b}{\hat{e}}_{b}\end{eqnarray} \tag{ 18 }$$

for the deflection of a photon as a function of impact parameter b. Here ${R}_{S}=\tfrac{2{{GM}}_{\odot }}{{c}^{2}}$ is the Sun's Schwarzschild radius and ${\hat{e}}_{b}$ is the unit vector in the direction of the impact parameter. Since a photon can travel from the source to the observer around the point mass and receive the same deflection for a path at any azimuthal angle, it is easy to see how gravitational lensing around a point source forms an Einstein Ring symmetrically around the point mass.

However, sophisticated helioseismological analyses (Mecheri et al. 2004) require the gravitational field to deviate from spherical symmetry. The rotation of the Sun creates an oblate gravitational potential that is described by a zonal spherical harmonic expansion

$$\begin{eqnarray}&&{\rm{\Phi }}=\displaystyle \frac{-{GM}}{r}\left[1+\sum _{n=2}^{\infty }{I}_{n}{\left(\displaystyle \frac{{R}_{\odot }}{r}\right)}^{n}{P}_{n}(\cos {\beta }_{S})\right],\end{eqnarray} \tag{ 19 }$$

where R_⊙ is the radius of the Sun, the P_n are the Legendre Polynomials (Abramowitz & Stegun 1968), β_S is the solar co-latitude, and I_n are the dimensionless multipole moments. Due to the North–South symmetry along the rotation axis, all terms with odd n are zero. The largest even n terms are I₂ ≈ 2 ×10⁻⁷, I₄ ≈ −4 × 10⁻⁹, I₆ ≈ −3 × 10⁻¹⁰, and I₈ ≈ 1 × 10⁻¹¹ (Roxburgh 2001).

The effect of this deviation from spherical symmetry on the resulting optical properties of the SGL was first investigated by Loutsenko (Loutsenko 2018). They demonstrated that the quadrupole moment I₂ results in the formation of an astroid-shaped caustic in the image plane with a diameter

$$\begin{eqnarray}&&{D}_{\mathrm{astroid}}=\displaystyle \frac{4{I}_{2}{R}_{\odot }^{2}}{\sqrt{2{R}_{S}r}}{\sin }^{2}{\beta }_{S}.\end{eqnarray} \tag{ 20 }$$

Turyshev and Toth (Turyshev & Toth 2021a) extended their work and computed the effect of any high-order multipole I_n on the resulting optical properties of the gravitational lens. In general, the caustic curves are given by higher-order hypocycloids (Turyshev & Toth 2021b), and the general case of photon deflection can be computed with

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{\theta }}}_{g} & = & \displaystyle \frac{2{R}_{S}}{b}\left[{\hat{e}}_{b}-\displaystyle \sum _{n=2}^{\infty }{I}_{n}{\left(\displaystyle \frac{{R}_{\odot }}{b}\right)}^{n}{\sin }^{n}({\beta }_{S})\right.\\ & & \times \left.\left({\hat{e}}_{b}\cos [n({\phi }_{\xi }-{\phi }_{S})]-{\hat{e}}_{{\phi }_{\xi }}\sin [n({\phi }_{\xi }-{\phi }_{S})]\right)\Space{0ex}{3.25ex}{0ex}\right],\end{array}\end{eqnarray} \tag{ 21 }$$

where ϕ_ξ and ϕ_S are the azimuthal orientations of the impact parameter and the solar rotation axis. When I_n = 0, this formula reduces to the well-known example of a point mass, but includes new terms for the deflection along the ${\hat{e}}_{{\phi }_{\xi }}$ direction, which is perpendicular to the impact parameter (see Figure 2 in Turyshev & Toth 2021a). This additional deflection gives the SGL entirely new behavior that is not symmetric with respect to azimuthal rotation, and explains geometrically the formation of the Einstein Cross for n = 2.

Furthermore, they go on to solve Maxwell's equations by treating the new gravitational potential as a perturbation of the Newtonian metric, and using the angular eikonal method, which takes inspiration from Mie scattering theory. In this formalism, the magnification of the lens μ, defined as the ratio of the Poynting vectors for a plane wave at the source and resulting image, is (Turyshev & Toth 2021c)

$$\begin{eqnarray}&&\mu ({\boldsymbol{x}})=2\pi {{kR}}_{S}| E({\boldsymbol{x}}){| }^{2},\end{eqnarray} \tag{ 22 }$$

where $k=\tfrac{2\pi }{\lambda }$ is the wavevector, and the quantity E( x ) is given by the diffraction integral (Turyshev & Toth 2021c)

$$\begin{eqnarray}&&\begin{array}{l}E({\boldsymbol{x}})=\displaystyle \frac{1}{2\pi }{\displaystyle \int }_{0}^{2\pi }d{\phi }_{\xi }\exp \left[-{ik}\left(\sqrt{\displaystyle \frac{2{R}_{S}}{r}}\rho \cos ({\phi }_{\xi }-\phi )\right.\right.\\ \quad +\,\left.\left.2{R}_{S}\displaystyle \sum _{n=2}^{\infty }\displaystyle \frac{{I}_{n}}{n}{\left(\displaystyle \frac{{R}_{\odot }}{\sqrt{2{R}_{S}r}}\right)}^{n}{\sin }^{n}({\beta }_{S})\cos [n({\phi }_{\xi }-{\phi }_{S})]\right)\right],\end{array}\end{eqnarray} \tag{ 23 }$$

where x = (ρ, ϕ) are the polar coordinates in the image plane. If all of the multipole moments I_n = 0, this diffraction integral can be analytically computed, and the result is that of the monopole-term only, $E({\boldsymbol{x}})={J}_{0}(k\sqrt{2{R}_{S}/r}\rho )$, where J₀ is the Bessel Function of the First Kind (Abramowitz & Stegun 1968). However, it is not analytically tractable for nonzero I_n , as would be the case for a realistic Sun, and so numerical methods must be relied upon to evaluate Equation (23).

In this paper, we investigate Equation (23) with a discrete approximation to the integral over the azimuthal angles ϕ_ξ . By replacing ${\int }_{0}^{2\pi }d{\phi }_{\xi }\exp [...]$ with a discrete sum indexed by j

$$\begin{eqnarray}&&\begin{array}{l}E({\boldsymbol{x}})\approx \displaystyle \frac{1}{2\pi }\displaystyle \sum _{j=1}^{{N}_{p}}{{\rm{\Delta }}}_{{\phi }_{\xi }}\exp \left[-{ik}\left(\sqrt{\displaystyle \frac{2{R}_{S}}{r}}\rho \cos ({\phi }_{\xi ,j}-\phi )\right.\right.\\ \,+\left.\left.2{R}_{S}\displaystyle \sum _{n=2}^{\infty }\displaystyle \frac{{I}_{n}}{n}{\left(\displaystyle \frac{{R}_{\odot }}{\sqrt{2{R}_{S}r}}\right)}^{n}{\sin }^{n}({\beta }_{S})\cos [n({\phi }_{\xi ,j}-{\phi }_{S})]\right)\right],\end{array}\end{eqnarray} \tag{ 24 }$$

it is possible to evaluate the integral numerically. Here ϕ_ξ,j ∈ [0, 2π] is a discrete approximation to the continuous ϕ_ξ and ${{\rm{\Delta }}}_{{\phi }_{\xi }}=\tfrac{2\pi }{{N}_{p}}$, where N_p is the number of points to use for the azimuthal sampling. N_p must be chosen carefully to be sufficiently large so that the approximation is valid. For "easy" cases where k, β_S , and ρ are relatively small, such as λ = 4 μm, β_S = 10°, and ρ = 0 m, we find that N_p ≳ 10⁴ is sufficient. "Hard" cases where k, β_S , and ρ are relatively large such as λ = 0.5 μm, β_S = 90°, and ρ = 10³ m can require significantly finer azimuthal sampling, up to N_p ≳ 10⁵ to ensure the discrete approximation to the integral does not break down.

The approximation in Equation (24) depends on many parameters, the location in the image plane x = (ρ, ϕ), the wavelength k = 2π/λ of light, the multiple moments I_n , the mass M_⊙ and radius R_⊙ of the Sun, and the orientation of the solar rotation axis (β_S , ϕ_S ). Furthermore, once one has computed the result of Equation (24), which gives the complex electromagnetic field in the SGL image plane rescaled by the magnification prefactors, it is rather straightforward to find the PSF that would result in the image plane of a telescope using the Fraunhofer approximation and the techniques of Fourier Optics (Hecht 2002; Goodman 1996). If P( x ) = ρ < D_tel/2 is the binary pupil transmission function of a circular aperture telescope, then the PSF in its image plane is

$$\begin{eqnarray}&&\mathrm{PSF}(\alpha ,\beta )={\left|{{ \mathcal F }}^{-1}(\alpha ,\beta )[P({\boldsymbol{x}})E({\boldsymbol{x}})]\right|}^{2},\end{eqnarray} \tag{ 25 }$$

where (α, β) are the angular coordinates in the conjugate plane, and ${{ \mathcal F }}^{-1}[...]$ is the inverse Fourier Transform. The inverse Fourier Transform, which has a positive sign in the exponent e^+ik... of the exponential basis, is used instead of the Forward Transform to remain consistent with observed orientation of asymmetric PSFs in high-contrast imaging from the ground (see Section 2.4 in Madurowicz et al. 2019). To showcase the effect of these parameters on the resulting images formed, Figure 7 was created. The figure contains three separate cases, with different wavelengths λ = (1, 4, 16) μm and solar co-latitude β_S = (11.25, 22.5, 45)°, with each case labeled in its title. Each case contains three plots, the magnitude of ∣E( x )∣, the phase inside the telescope pupil, and the resulting PSF in the telescope image plane.

Zoom In Zoom Out Reset image size

A few key observations from this Figure. First, the image formation in the SGL image plane is dominated by the quadrupole moment I₂, and the resulting ∣E∣ is shaped primarily like the four-pointed hypocycloid, the astroid. The higher-order multipole moments have a small effect, but are subdominant for realistic solar values. Second, for larger β_S , the size of the astroid caustic grows significantly. For the first case with β_S = 1125, the astroid diameter is ∼20 m, but for β_S = 45°, it is ∼200 m. At the solar equator, β_S = 90°, the astroid diameter reaches its maximum value ∼560 m. This has the implication that any given target planet, which has a unique location on the sky, will suffer from degradation to the average magnification over the pupil for a set telescope diameter, since the positioning and size of the telescope control the fraction of intercepted light. However, this only is true when the telescope diameter is much smaller than the astroid diameter. If the telescope is sufficiently large to capture the entire astroid, then all of the flux can be intercepted, but that would require telescopes much larger than is currently feasible.

Third, the phase of the electromagnetic field appears to have a grid-like pattern that is responsible for the formation of the Einstein Cross. The size scaling of this grid pattern is controlled mostly by the wavelength of the light, with shorter wavelengths having a narrower spacing between the zones, and longer wavelengths having larger zones. It is crucially important that the sampling of the phase inside the pupil is fine enough that the delicate grid pattern is not aliased, and that the pupil area is large enough so that at least one "grid element" is fully contained inside the pupil. If these conditions are not both met, the resulting PSF does not produce the characteristic Einstein Cross, seen on the right. Fourth, the angular size of the four points in the telescope's PSF grow linearly with wavelength, as is expected. The positioning of the four point images is just outside the angular extent of the solar disk, shown at a heliocentric distance of r_SGL = 650 au. When the telescope is precisely aligned with the optical axis of the point source and the gravitational lens, the locations of the four point sources is symmetric; however, a lateral misalignment of the telescope intercepts an offset phase pattern, which will cause the locations of the points to be perturbed around the edge of the ring; see Figure 8. When the telescope is positioned directly on the caustic edge, the images have moved sufficiently around the edge to merge and form into an arc. When the telescope has moved outside the astroid caustic, the result is two point sources similar to the monopole-only case. However, the effect of the multipoles is still present, and the intensity of each point is asymmetric.

Zoom In Zoom Out Reset image size

To further investigate the effect of the telescope diameter on the intercepted flux, we compute the average magnification across the telescope pupil 〈μ〉_P . Equation (24) is computed numerically for points (ρ, ϕ) in the image plane, and the average across the pupil is defined as the weighted sum of μ in the region inside the telescope pupil $P\in \{(x,y)| \sqrt{{\left(x-{c}_{x}\right)}^{2}+{\left(y-{c}_{y}\right)}^{2}}\leqslant {D}_{\mathrm{tel}}/2\}$.

$$\begin{eqnarray}&&\langle \mu ({\boldsymbol{x}}){\rangle }_{{\rm{P}}}=\displaystyle \frac{\mathop{\iint }\limits_{{\rm{P}}}\mu ({\boldsymbol{x}})d{\boldsymbol{x}}}{\mathop{\iint }\limits_{{\rm{P}}}d{\boldsymbol{x}}}\end{eqnarray} \tag{ 26 }$$

Here (c_x , c_y ) is the location of the center of the telescope in Cartesian coordinates, which has diameter D_tel = 2 m. This metric is used to compare the multiple scenarios plotted in Figure 9.

Zoom In Zoom Out Reset image size

Examining the behavior of the magnification in Figure 9, it is apparent that the effect of the higher-order multipole moments causes significant alteration to the magnification in the region of the astroid caustic. This regime is characterized by ρ ≲ D_astroid. Here the maximum of the magnification no longer occurs on the optical axis, but instead in the region of the caustic, and is reduced from the maximum of the monopole by a few orders of magnitude, of which the exact amount depends on the solar co-latitude β_S as well as the telescope diameter D_tel.

However, far outside the astroid caustic, the average magnification is consistent with the pure monopole case. In this section, we investigate the asymptotic behavior of these functions and demonstrate an analytic comparison between the wave and geometric theory, which results in a power-law approximation valid in this regime. Using the geometric theory, the magnification of a point-mass lens can be described by Narayan & Bartelmann (1996)

$$\begin{eqnarray}&&{\mu }_{\mathrm{geo}}=\displaystyle \frac{1}{1-{\left(\tfrac{{\theta }_{E}}{{\theta }_{\pm }}\right)}^{4}},\end{eqnarray} \tag{ 27 }$$

where θ_E is the angular scale of the Einstein Ring

$$\begin{eqnarray}&&{\theta }_{E}=\sqrt{2{R}_{S}\displaystyle \frac{{D}_{{LS}}}{{D}_{S}{D}_{L}}},\end{eqnarray} \tag{ 28 }$$

and the θ_± are the angular coordinates of the magnified images

$$\begin{eqnarray}&&{\theta }_{\pm }=\displaystyle \frac{1}{2}\left(\beta \pm \sqrt{{\beta }^{2}+4{\theta }_{E}}\right).\end{eqnarray} \tag{ 29 }$$

Here β is the angular coordinate of the point source and is geometrically related to the radial coordinate ρ in the image plane by

$$\begin{eqnarray}&&\beta =\arctan \displaystyle \frac{{D}_{S}}{\rho }-\arctan \displaystyle \frac{{D}_{L}}{\rho }\approx \rho \left(\displaystyle \frac{1}{{D}_{L}}-\displaystyle \frac{1}{{D}_{S}}\right).\end{eqnarray} \tag{ 30 }$$

Taking a Taylor expansion of Equation (27) in the regime where β ≪ θ_E , results in

$$\begin{eqnarray}&&{\mu }_{\mathrm{geo}}\approx \displaystyle \frac{1}{2}+\displaystyle \frac{{\theta }_{E}}{2\beta },\end{eqnarray} \tag{ 31 }$$

so evidently the dominant scaling of the geometric magnification is ${\mu }_{\mathrm{geo}}\propto \tfrac{1}{\beta }\propto \tfrac{1}{\rho }$. Taking a Taylor expansion in the regime where β ≫ θ_E results in

$$\begin{eqnarray}&&{\mu }_{\mathrm{geo}}\approx 1+{\left(\displaystyle \frac{{\theta }_{E}}{\beta }\right)}^{4},\end{eqnarray} \tag{ 32 }$$

and the dominant term is the constant, which is succinctly summarized in the limiting behavior ${\mathrm{lim}}_{\beta \to \infty }{\mu }_{\mathrm{geo}}=1$. This limit indicates that far enough away from the lens, the point source simply exists as it would be unperturbed by anything. To compare and contrast this geometric scaling argument with the wave-theoretic treatment described previously, we take an asymptotic expansion of the magnification for a pure monopole SGL. Starting with the expression for the magnification

$$\begin{eqnarray}&&{\mu }_{\mathrm{wave}}=2\pi {{kR}}_{S}{\left|{J}_{0}\left(k\sqrt{\displaystyle \frac{2{R}_{S}}{r}}\rho \right)\right|}^{2},\end{eqnarray} \tag{ 33 }$$

and since the Bessel Functions of the First Kind have asymptotic expansions (Abramowitz & Stegun 1968)

$$\begin{eqnarray}&&{J}_{n}(z)=\sqrt{\displaystyle \frac{2\pi }{z}}\left[\cos \left(z-\displaystyle \frac{n\pi }{2}-\displaystyle \frac{\pi }{2}\right)+{e}^{| \mathrm{Im}(z)| }O\left(\displaystyle \frac{1}{| z| }\right)\right],\end{eqnarray} \tag{ 34 }$$

we can conclude, for ρ ≳ D_astroid, that

$$\begin{eqnarray}&&{\mu }_{\mathrm{wave}}\approx 4{\pi }^{2}\sqrt{\displaystyle \frac{{{rR}}_{S}}{2}}{\rho }^{-1}{\cos }^{2}\left(k\sqrt{\displaystyle \frac{2{R}_{S}}{r}}\rho -\displaystyle \frac{\pi }{2}\right)\propto \displaystyle \frac{1}{\rho },\end{eqnarray} \tag{ 35 }$$

which shares the same power-law scaling with the geometric theory. Incidentally, the limiting behavior of this function is discrepant as ${\mathrm{lim}}_{\rho \to \infty }{\mu }_{\mathrm{wave}}=0$. This likely is a result of the formalism being used to describe the electromagnetic field in the geometric shadow of the Sun, where the direct path component of the field is obscured by the solar disk. So if ρ > ρ_critical, such that β > θ_E , then this description is no longer valid. Indeed, in this case, it would be simpler to imagine that the lens does not exist at all, and to treat the direct path on its own. To summarize, the monopole approximation is valid and consistent with the geometric power-law scaling ∝ ρ⁻¹ in the regime where D_astroid ≲ ρ < ρ_critical. Inside the astroid, one must consider the perturbations of the higher-order multipoles, and outside the critical distance, one may ignore the lensing entirely. This approximation is useful to quickly evaluate the response of a realistic SGL at large lateral separations, since numerical evaluation of Equation (24) becomes increasingly challenging at large ρ.

With formalism described in Sections 3.1 and 3.2, it is possible to calculate the PSF of the SGL as well as the PSF as it would appear in the image plane of a telescope in the focal plane of the SGL. In this section, we extend our analysis to the imaging of extended sources under the assumptions of incoherent emission and small angular extent of the source. The former assumption allows one to treat the extended source as a superposition of independent point sources, and the latter allows one to ignore minor perturbations to the orientation of the solar rotation axis (β_S , ϕ_S ) that would begin to matter if the angular extent of the extended source was sufficiently large. For justification, an Earth-sized planet at a target distance of d_S = 100 pc occupies a solid angle of ${{\rm{\Omega }}}_{E}=\pi {\left(\tfrac{{R}_{E}}{{d}_{S}}\right)}^{2}\approx 1.4\,\times {10}^{-23}$ steradians.

If the Cartesian coordinates in the source plane are $(x^{\prime} ,y^{\prime} )$, then we can approximate the continuous intensity distribution of the target planet $I(x^{\prime} ,y^{\prime} )$ with

$$\begin{eqnarray}&&{I}_{{n}_{x},{n}_{y}}\approx \langle {{\rm{III}}}_{{n}_{x}}(x^{\prime} ){{\rm{III}}}_{{n}_{y}}(y^{\prime} )I(x^{\prime} ,y^{\prime} )\rangle ,\end{eqnarray} \tag{ 36 }$$

where ${{\rm{III}}}_{{n}_{x}}(x^{\prime} )=\delta (x^{\prime} -{{\rm{\Delta }}}_{x}({n}_{x}-{N}_{x}/2))$ is a comb function, δ is the Dirac delta function, n_x ∈ [0, 1, 2, ..., N_x ] is the integer index of the pixel, Δ_x is the spatial sampling, and the N_x /2 is a centering offset, so the origin of $x^{\prime} =0$ occurs at the halfway index n_x = N_x /2. Additionally, 〈...〉 represents the average of the intensity distribution in the vicinity of the sampling point, and at the end, ${I}_{{n}_{x},{n}_{y}}$ is a vector (or matrix), indexed by pixel numbers (n_x , n_y ), which represents a discrete approximation to the continuous intensity distribution of the source. In our analysis, we use uniform sampling for both x and y and set Δ_x = 2R_Earth/N_x , N_x = 64 pixels. As a consequence,

$$\begin{eqnarray}&&x{{\prime} }_{{n}_{x}}={R}_{\mathrm{Earth}}\left(2\displaystyle \frac{{n}_{x}}{{N}_{x}}-1\right)\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&y{{\prime} }_{{n}_{y}}={R}_{\mathrm{Earth}}\left(2\displaystyle \frac{{n}_{y}}{{N}_{y}}-1\right),\end{eqnarray} \tag{ 38 }$$

and both $x^{\prime} ,y^{\prime} \in [-{R}_{\mathrm{Earth}},{R}_{\mathrm{Earth}}],$ since $\tfrac{{n}_{x,y}}{{N}_{x,y}}\in [0,1]$.

For our analysis, each pixel in this approximation acts as an independent, incoherent source in the imaging problem. Each point on the source object is reimaged by the SGL and appears as an astroid-shaped PSF in the SGL image plane, with the center of the astroid offset by its location in the source plane multiplied by the SGL plate scale. If the Cartesian coordinates in the SGL image plane are $({x}_{{n}_{x}},{y}_{{n}_{y}})$, then the location of the astroid center for a point source located at $(x{{\prime} }_{{n}_{x}},y{{\prime} }_{{n}_{y}})$ is

$$\begin{eqnarray}&&{x}_{{n}_{x}}=-x{{\prime} }_{{n}_{x}}\displaystyle \frac{{d}_{L}}{{d}_{S}},\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&{y}_{{n}_{y}}=-y{{\prime} }_{{n}_{y}}\displaystyle \frac{{d}_{L}}{{d}_{S}},\end{eqnarray} \tag{ 40 }$$

where the negative sign appears due to the inversion of the image. Figure 10 depicts the geometry of the problem for convenience.

Zoom In Zoom Out Reset image size

In this diagram, which is not to scale, the size of the astroid PSF in the SGL image plane is depicted to be a fraction of the size of the extended object, a diagram of the Earth. One can visualize how a telescope situated directly on the optical axis will see unique PSFs for each of the three points highlighted in the diagram. For the central point on the optical axis, the telescope will sit squarely inside the center of the astroid, where the PSF will form an Einstein cross in the telescope image plane. For the two points that are offset from the optical axis, the telescope centered on the optical axis will instead intercept light from outside the edge of the astroid, which will result in simply two point sources in the telescope's image plane. It may be helpful to again reference Figure 8, which shows the PSFs that form in the telescope image plane for telescope locations offset relative to the astroid.

To obtain the final intensity distribution in the telescope image plane, we simply sum the contributions of each point in the source plane, using the unique PSF that a telescope would see corresponding to that point. To do this, we compute the PSF over a grid of telescope positions offset relative to the optical axis. It is equivalent to consider a point source located at $(x^{\prime} ,y^{\prime} )$ in the source plane and a telescope centered on the optical axis (x = 0, y = 0) as it would be to consider a point source on the optical axis $(x^{\prime} =0,y^{\prime} =0)$, and have the telescope center offset by negated and plate-scale-reduced coordinates $({c}_{x}=x{{\prime} }_{{n}_{x}}\tfrac{{d}_{L}}{{d}_{S}},{c}_{y}=y{{\prime} }_{{n}_{y}}\tfrac{{d}_{L}}{{d}_{S}})$. Thus the relationship between source plane pixel indices and corresponding telescope centering offsets is

$$\begin{eqnarray}&&{c}_{x}=\displaystyle \frac{{d}_{L}}{{d}_{S}}{R}_{\mathrm{Earth}}\left(2\displaystyle \frac{{n}_{x}}{{N}_{x}}-1\right)\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&{c}_{y}=\displaystyle \frac{{d}_{L}}{{d}_{S}}{R}_{\mathrm{Earth}}\left(2\displaystyle \frac{{n}_{y}}{{N}_{y}}-1\right).\end{eqnarray} \tag{ 42 }$$

With the vectors of centering offsets forming a grid of telescope locations (c_x , c_y ), we compute the PSF in the telescope image plane using Equation (25), where the pupil transmission function is modified so that $P({\boldsymbol{x}})\,\in \{(x,y)| \sqrt{{\left(x-{c}_{x}\right)}^{2}+{\left(y-{c}_{y}\right)}^{2}}\leqslant {D}_{\mathrm{tel}}/2\}$ and the telescope is properly centered around (c_x , c_y ). As a computational shortcut, E( x ) in Equation (24) is only computed for points x inside the pupil where P ≠ 0 to save on computational cost. Additional computational cost is saved by reverting to the analytic monopole expression for $E({\boldsymbol{x}})={J}_{0}(k\sqrt{2{R}_{S}/r}\rho )$ when the telescope is sufficiently far from the optical axis. Since Figure 9 demonstrates that the monopole approximation starts to be very accurate for $\rho =\sqrt{{c}_{x}^{2}+{c}_{y}^{2}}\gtrsim 4{D}_{\mathrm{astroid}}$, this is the condition we check before falling back on the fast analytic computation for E( x ). The final intensity distribution I_extended is computed as a matrix product of the PSF grid ${\mathrm{PSF}}_{{c}_{x},{c}_{y}}$ and the source image vector ${I}_{{n}_{x},{n}_{y}}$.

$$\begin{eqnarray}&&{I}_{\mathrm{extended}}(\alpha ,\beta )={\mathrm{PSF}}_{{c}_{x},{c}_{y}}(\alpha ,\beta ){I}_{{n}_{x},{n}_{y}}\end{eqnarray} \tag{ 43 }$$

Thus with the PSF grid previously computed and in memory, the forward model of imaging an extended source can be represented as a linear transformation of the input intensity distribution. However, for this to be achieved with a single matrix operation, the two-dimensional indices of source coordinates indexed by pixel numbers (n_x , n_y ) and the point-source offsets corresponding to telescope centering offsets (c_x , c_y ) must be unraveled into a single monotonically increasing index. This choice of ordering is arbitrary and can be done in a combinatorially large number of unique ways, although we find unraveling indices one row at a time to be the simplest. Figure 11 was created to demonstrate imaging of extended sources.

Zoom In Zoom Out Reset image size

In general, the size of the astroid is controlled by the quadrupole moment I₂ and the solar co-latitude β_S , see Equation (20), while the size of the image D_image is controlled by the physical dimension of the source, in this case R_Earth, as well as the geometry of the problem, including the distances between the source and lens d_S , and the telescope and the lens d_L . The exact relation is

$$\begin{eqnarray}&&{D}_{\mathrm{image}}=2{R}_{\mathrm{Earth}}\displaystyle \frac{{d}_{L}}{{d}_{S}}.\end{eqnarray} \tag{ 44 }$$

For extended sources of different sizes, there are three regimes that are important to consider. One limit is that of very small sources, or D_image << D_astroid. In this case, one may think of the extended objects as essentially just a point, and the resulting image formed will be that of the PSF. The other limit is that of very large objects, or D_image >> D_astroid. In this case, the majority of point sources in the SGL image plane will be farther from the optical axis than the dimension of the astroid, and the resulting PSF in the telescope image plane for each location on the extended object will simply be two points. One can think of this limit as being equivalent to taking D_astroid → 0, or β_S → 0, where the telescope is situated directly above or below the solar rotational axis and the projection of the quadrupole moment into the plane of observation becomes negligible. In this limit, the higher-order effects from the multipole moments can be neglected, as nearly all points appear to be reimaged into two points in the telescope image plane, which is identical to the monopole-only analysis described in Madurowicz (2020).

The last regime is that where the dimension of the reimaged target and astroid have similar sizes, or D_image ∼ D_astroid. In this regime, one cannot ignore the higher-order effects of the multipoles, as many of the points on the target will be reimaged multiple times, some of the points may sit on the critical caustic, and some of the points may be outside the astroid entirely. Thus, the final intensity distribution formed in the telescope image plane is a superposition of the multifaceted PSF distribution and intensity distribution from the source.

Actual physical observations using the SGL will face an extreme difficulty dealing with contamination of source signal with noise from multiple sources along the line of sight, including the Sun and its corona, zodiacal dust in the solar system, and additional nontarget sources in the SGL field of view (FOV). This section investigates the expected brightness from the Sun itself as a primary source of noise, compared to the expected brightness of a potential, Earth-like, target source by collecting data and models built from observations of the Earth and Sun and scaling them to match the geometry of the SGL problem.

The first Earth data set used in this study is a series of images taken from the Advanced Baseline Imager on board the GOES-R weather satellite (GOES-R Calibration Working Group & GOES-R Series Program 2017). The data includes extremely high-resolution (either ∼5 k by ∼5 k or ∼10 k by ∼10 k, depending on wavelength) images of the Earth taken from geostationary orbit, across a wide range of wavelengths, from 0.47 μm in the visible light spectrum, to 13.27 μm in the mid-infrared. The specific data set used in this study is the level-1b radiances product, which has been processed and calibrated to represent authentic top-of-the-atmosphere radiances in units of [W m⁻² sr⁻¹ μm] or [mW m⁻² sr⁻¹ (cm⁻¹)⁻¹] depending on if the channel is dominated by reflection of sunlight or thermal emission from the Earth itself. The GOES-R data is plotted in Figure 12.

Zoom In Zoom Out Reset image size

While the GOES-R data will be useful to investigate the effect of nonuniform intensity distributions across different wavelengths, the wavelength coverage is not complete. With the goal of investigating higher spectral resolutions, we investigate the Earth model described in Robinson et al. (2011). The Earth spectral model describes disk-integrated specific intensities of the Earth across a massive range in wavelength space, from 0.1 μm to 200 μm, with spectral resolution $R=\tfrac{\lambda }{{\rm{\Delta }}\lambda }$ ranging from ∼100 at the longest wavelengths to ∼10⁵ at the shortest. The Earth spectral model also includes variations in brightness due to the inclination of the orbit and the orbital phase. Both the GOES-R data and the Earth spectral model are plotted in Figure 13 in units of spectral flux density or [W m⁻² μm] by multiplying by ${\int }_{0}^{2\pi }{\int }_{0}^{\pi }\sin \theta \cos \theta d\theta =\pi$ sr. These units will be convenient for comparison with a model of the spectral flux density of the Sun.

Zoom In Zoom Out Reset image size

To compare the brightness of the Earth and the Sun across many wavelengths, we use three separate solar models to act as validation and ensure consistency. The first and simplest model is that of a spherical blackbody with radius R_⊙ and temperature T = 5778 K, to match the effective temperature of the solar photosphere (Carroll & Ostlie 2007). The blackbody spectrum is governed by the Planck equation

$$\begin{eqnarray}&&{B}_{\lambda }(T)=\displaystyle \frac{2{{hc}}^{2}}{{\lambda }^{5}}{\left({e}^{\tfrac{{hc}}{\lambda {kT}}}-1\right)}^{-1},\end{eqnarray} \tag{ 45 }$$

where h is Planck's constant, c is the speed of light, k is the Boltzmann constant, and T is the blackbody temperature. Since the blackbody equation gives specific intensity at the photosphere, to calculate the flux density at the top of Earth's atmosphere, one must rescale by the solid angle of the Sun as seen from Earth

$$\begin{eqnarray}&&{F}_{\lambda }={B}_{\lambda }\pi {\left(\displaystyle \frac{{R}_{\odot }}{1\mathrm{au}}\right)}^{2}.\end{eqnarray} \tag{ 46 }$$

The blackbody model is quite a good description of the radiation from the Sun, due to its physical elegance and analytic form, but is incomplete due to absorption in the solar atmosphere and additional complications. However, the Sun is bright and easily observable, so spectral data at high resolution is plentiful and of high significance. We use two separate solar spectral models that are validated with real observations. The first model is the National Oceanic and Atmospheric Administration (NOAA) Climate Data Record (CDR) solar standard irradiance (Coddington et al. 2015), and the second is the American Society for Testing and Materials (ASTM) E-490 solar standard, which is compiled from multiple unique observations (Kurucz 1993; Fröhlich & Lean 1998; Woods et al. 1996; Smith & Gottlieb 1974) to serve as a reference spectrum across a large range of wavelengths. Both of these models are converted into identical units for comparison in Figure 13.

A few observations regarding the spectral models. The Sun is accurately described by a blackbody spectrum, with small perturbations due to absorption and emission lines that are mostly invisible at the plotted spectral resolution. The Earth model has two peaks, one in the visible spectrum corresponding to the reflection of sunlight, and one in the infrared corresponding to the thermal emission from the planet itself. Additionally, the Earth's spectrum is significantly altered by absorption in its atmosphere, with many valleys in different wavelength regions indicative of particular molecules such as H₂O and CH₄.

The previous figure depicts the spectral flux density that an observer would see for the Sun and the Earth if that observer was located at the top of Earth's atmosphere. However, this arrangement is not equivalent to the geometry of imaging a target planet through the SGL, and so these curves must be rescaled to the appropriate distance using the inverse square law. Additionally, the flux of the planet is magnified by imaging through the SGL, by an amount that depends on its location on sky as well as by the exact alignment between telescope, target, and the Sun. Figure 9 shows that reasonable values of the magnification can range from μ = 10⁷ to 10⁹. Additionally, one may obscure the light from the Sun using a starshade (Cash 2011) or coronagraph (Lyot 1939), which may need to be specifically designed to obscure the extended disk of the Sun (Ferrari et al. 2009), but could significantly reduced the observed flux from the Sun.

To make a fair comparison of the brightnesses with all of these factors considered, Figure 14 was created. In this Figure, there are three curves for the Sun and the Earth each. Each starts with the Robinson Earth spectral model or the ASTM-E solar spectral model and rescales the flux amplitude depending on different assumptions. From top to bottom, the solar curves are rescaled to d_L = 650 au without a coronagraph, d_L = 650 au with a coronagraph with a contrast of 10⁻¹⁰, and d_L = 10 pc behind a similar coronagraph. The rescaling factor is simply a multiplication by $\mathrm{contrast}\times {\left(\tfrac{{d}_{L}}{1\mathrm{au}}\right)}^{2}$. This level of coronagraphic sensitivity has not yet been reached but is projected to be possible by the HABEX mission concept (Gaudi et al. 2018). The three curves for the Earth, from top to bottom, all correspond to a target distance d_S = 10 pc, but with μ = 10⁹, 10⁷, or 0, corresponding to a possible magnified source behind the SGL and an unmagnified source that one would see without the SGL.

Zoom In Zoom Out Reset image size

The principal observation this Figure was designed to showcase is that in terms of flux contrast between the target source and the primary noise source, using the SGL is comparable to direct imaging. The spectral flux density essentially describes the total number of photons that would be captured by the telescope over some bandpass, but is agnostic to the distribution of that light in the focal plane of the telescope, which is discussed in 4.4. The dashed pair of curves correspond to a reasonable expectation of the flux for the standard SGL, while the dotted pair of curves correspond to a reasonable comparison for direct imaging. The ratio of the flux between these two objects gives the contrast as a function of wavelength, which is quite similar. However, the major advantage that the SGL arrangement has over direct imaging is in the absolute brightness of both the target and noise source, being approximately 10⁷ times brighter for both the target and the solar/stellar noise contribution. This advantage may prove useful to obtain extremely high-dispersion spectra of an exoplanetary atmosphere, since high dispersion implies correspondingly smaller signal to noise per bin, and high-dispersion instruments typically struggle with low throughput (Vogt et al. 1994).

Aside from the light from the target source and the Sun itself, the Einstein ring contains light from all astrophysical sources in the SGL's FOV. This includes primarily the light from the target planet's host star, but also includes the probability of faint distant sources that are coincidentally close to alignment with the optical axis. In Figure 9 in Section 3.3, we show the average magnification 〈μ〉_P over the telescope pupil as a function of the telescope's radial offset in the SGL image plane $\rho =\sqrt{{c}_{x}^{2}+{c}_{y}^{2}}$. For radial separations significantly greater than the dimension of the astroid, the average magnification is a power law with index −1, matching the expectation for a pure monopole SGL.

$$\begin{eqnarray}&&\mathrm{log}\langle \mu {\rangle }_{P}\approx 8.88-\mathrm{log}\left(\displaystyle \frac{\rho }{1{\rm{m}}}\right),\quad \rho \gg {D}_{\mathrm{astroid}}\end{eqnarray} \tag{ 47 }$$

Since a host star at d_S = 10 pc with a target planet at a = 1 au has an angular separation in the source plane $\delta \theta =\tfrac{a}{{d}_{S}}=0.1$ arcsec, this is equivalent to a telescope radial offset in the SGL image plane of $\rho ={d}_{L}\delta \theta =a\tfrac{{d}_{L}}{{d}_{S}}\sim 4.7\times {10}^{7}$ m and thus an average magnification of 〈μ〉_P ∼ 16. Figure 15 demonstrates the relative contribution across various wavelengths the host star will have to the observed flux of the target planet and the post-coronagraphic solar contribution.

Zoom In Zoom Out Reset image size

Examining Figure 15, it is apparent that the target planet's host star contributes a relatively significant amount of the total flux intercepted. In optical wavelengths, the total contrast between the target planet and its host star is about four orders of magnitude, while at thermal IR wavelengths, it is closer to one order of magnitude. This flux will be spatially separated in the telescope image plane, and as the host star is far off-axis, it will be reimaged into two points on opposite sides of the Einstein ring, similar to in the first panel in Figure 11, which could be separated from the light from the planet that forms the entire ring.

The target planet's host star forms the primary source of contamination inside the Einstein Ring, due to its relatively nearby separation from the target planet at ∼.1''. However, astrophysical false positives are known to appear in direct imaging, such as this one case of a white dwarf companion to a K-type star (Zurlo et al. 2013). Other astrophysical objects may appear in the FOV despite being gravitationally unbound to the host star, although this can be mitigated by multiepoch observations demonstrating common proper motion (Macintosh et al. 2015). Recent work investigating the probabilities of observing a galactic or extragalactic background source from the Hubble Source catalog, Gaia catalog, and Besancon galaxy simulations (Cracraft et al. 2021) demonstrate that the probability of such a background object is usually low, though it increases at lower galactic latitude.

A mission to the focus of the SGL naturally occurs inside of the Milky Way and thus must be prepared for the possibility that a background object appears inside the SGL's FOV and contaminates the Einstein Ring. To investigate this probability, we rescale star counts from the Besancon galaxy simulation (Czekaj et al. 2014) across the entire sky. Figure 16 shows an interpolated map of star counts of Johnson V-mag < 21 within 1 deg² converted into the probability of finding such a star within 1 arcsec².

Zoom In Zoom Out Reset image size

Examining the map in the Figure, the probability of some background object ranges by multiple orders of magnitude depending on sky location. In the direction of the galactic center, the probability is nearly certain that at least one such background star will be visible, while away from the galactic center, it may be as low as 1 in 10,000. The magnification of such a background object depends strongly on the angular separation between it and the target planet (Equation (47)), so it is difficult to quantify exactly how bright such a source will appear. If the background object has an angular separation of δ θ = 1'', this is equivalent to a telescope radial offset in the SGL image plane of ρ = d_L δ θ ∼ 4.7 × 10⁸ m and thus an average magnification of 〈μ〉_P ∼ 1.6, which is a factor of 10 less magnification than that of the host star at 0.1''. If the background object has an intrinsically low luminosity or is significantly farther away in distance than the target, the total flux may or may not be large compared to the light from the planet. There is always a small possibility that the background object coincidentally appears at a small angular separation and becomes further magnified by the SGL. However, the light from the background object will be necessarily restricted to two diametrically opposed points on the Einstein Ring, similarly to the host star, with an azimuthal angle corresponding to the position angle of the background object. It can be spatially mitigated in the telescope's image plane by separating the azimuthal elements of the Einstein ring.

Astrophysical backgrounds aside, there will be additional contamination of light from astrophysical foregrounds that must be mitigated to obtain precise measurements of the Einstein ring. We consider here primarily the light from the solar corona and zodiacal light from dust populations in the solar system. Since the dawn of civilization, total solar eclipses have provided a natural occultation of the solar disk and allowed people to observe in detail the solar corona (de Jong & van Soldt 1989; Dreyer 1877). Even in the modern day, with the development of the coronagraphic instruments (Lyot 1939), eclipse observations still provide some of the only detailed observations of the solar corona right up to the edge of the solar disk (Pasachoff et al. 2018; Judge et al. 2019). These eclipse observations demonstrate the azimuthally asymmetric nature of the coronal emission, which complicates the analysis of the azimuthal variations in the Einstein ring's structure.

Coronagraphic instruments are typically limited by an inner working angle, which prevents observations at small separations, and so cannot observe the corona near the edge of the solar disk. However, unlike eclipse observations, coronagraphic instruments on space-borne telescopes can provide continuous monitoring of the solar coronal environment. These space-borne observations enable investigation of time-dependent coronal phenomena, such as the relationship between coronal emission and changing magnetic field patterns related to the solar cycle (Antonucci et al. 2020).

Spectroscopic observations (Moore 1934) of the solar corona indicate that the light emanating from this region is the result of multiple distinct processes, the scattering of sunlight off of dust particles, scattering off of free electrons, and direct emission from the recombination of ionized atoms (Zessewitsch & Nikonow 1929; Mierla et al. 2008). These unique components have been termed the F-, K-, and E- coronas, respectively. F- is for Fraunhofer, who discovered the similarity of absorption lines between the Sun's photosphere and corona, indicating the light is reflected by large particles. K- is for "continuous" in German, as Doppler broadening of the photosphere absorption lines from the velocity of free electrons creates the appearance of a continuum spectrum. E- is for "emission." The F- and K- coronal contributions can be reliably separated using the difference in spatial dependence between dust and free electrons (Calbert & Beard 1972), while the E- corona may be spectroscopically mitigated due to the narrow emission lines.

The F- corona is largely considered to be a natural extension of the zodiacal light at smaller angular separations or elongations with regard to the Sun, and contains a wavelength dependence that increases toward long wavelengths due to thermal emission from the grains themselves in addition to the reflection of sunlight (Kimura & Mann 1998). However, interpretation of the F- coronal brightness and the corresponding dust distribution is complicated by inversion of the line-of-sight integral. Theoretical predictions for the near solar environment have speculated that the F- corona may decrease strength in the vicinity of the Sun (Lamy 1974; Mukai et al. 1974), due to orbital differentiation from Poynting–Robertson drag, pressure drag, radiation pressure, and sublimation of the grains. Observations from the WISPR instrument on board the Parker Solar Probe (Howard et al. 2019) have provided the first observation of such an F- coronal decrease, inferred from brightness measurements at various elongations at different heliocentric distances between 0.166 and 0.336 au.

Extensive brightness and color measurements of the solar corona have been taken and compiled previously, and for this work, we cross-reference multiple sources in the literature. By combining the data from Figure 2 in Mann (1992), Figure 1 in Kimura & Mann (1998), Figure 1 in Howard et al. (2019), and Figure 9 in Boe et al. (2021), we are able to analyze the solar coronal brightness across a wide range of elongations . The data are reproduced in Figure 17(a) for convenience. To build a model of the solar coronal brightness, we perform a joint analysis of the various data sets, by fitting for the maximum likelihood power-law model of the data

$$\begin{eqnarray}&&\mathrm{log}\left(\displaystyle \frac{{B}_{\mathrm{corona}}}{{B}_{\odot }}\right)=A+a\mathrm{log}\left(\displaystyle \frac{\epsilon }{{R}_{\odot }}-1\right),\end{eqnarray} \tag{ 48 }$$

which is linear in log coordinates, and by choosing the base of the exponent as $\tfrac{\epsilon }{{R}_{\odot }}-1$, the edge of the solar disk at = R_⊙ corresponds to the origin of the power law. The resulting posterior distributions on the power-law slope a and intercept A are shown in Figure 17(b). The covariances of the inferred parameters show constraints given by different individual data sets, while the joint analysis is restricted to the region of overlap between all data sets. The best-fit parameters as a result of this analysis are A = −7.836 ± 0.05 and a = −2.071 ±0.055. The resulting power-law curve is shown in black in Figure 17(a).

Zoom In Zoom Out Reset image size

Most of these measurements are taken with visible light wavelengths, although some effort has been made to make measurements in the near-infrared as well, and the results are compiled into Figures 5 in Boe et al. (2021) and Figure 2 in Kimura & Mann (1998). In Figure 2 of Kimura & Mann (1998), the scatter is large; there is a slight increase at longer wavelengths, going from roughly 10⁻⁹ B_⊙ at λ = 0.3 μm to 10⁻⁸ B_⊙ at λ = 3 μm for elongations of = 4R_⊙. Figure 5 of Boe et al. (2021) reports relative intensities over λ ∈ [0.5, 0.8] μm, and measures a spectral power-law index of B_corona ∝ λ^0.91±0.07. This reddening effect could be the result of thermal emission from the particles themselves, as well as preferential directional scattering of thermal emission from the Sun.

Figure 7 in Boe et al. (2021) also demonstrates the azimuthal variation in intensity, with isointensity curves fluctuating from ∈ [2.0, 2.4]R_⊙ at 10⁻⁸ B_⊙ to ∈ [1.3, 1.5]R_⊙ at 10⁻⁷ B_⊙. To show the azimuthal variation that can reasonably be expected here, we plot data from the HAO COSMO K-coronagraph (de Wijn et al. 2012; Hou et al. 2013) in Figure 18 with a black-hot color map.

Zoom In Zoom Out Reset image size

The data in Figure 18 are level-2 PB calibrated data downloaded directly from the HAO website. Some of the many calibration steps include camera nonlinearity, dark, and gain corrections, polarization demodulation, distortion correction, combining images from multiple cameras, removal of sky polarization, and correction for daily variation in sky transmission. From examining the image, it is clear that the observations are roughly in correspondence with the power-law analysis described previously, with the coronal brightness ranging from 10⁻⁶ B_⊙ at ∼ R_⊙ to 10⁻⁸ B_⊙ at ∼ 2R_⊙, although it is clear that the coronal brightness varies depending on azimuthal angle.

4.4.1. Coronagraph Design

Assessing the feasibility of imaging the Einstein ring at extremely high contrast is complicated by the redistribution of light in angular coordinates. The Einstein Ring has an approximate radius of ∼1'', while a host star/target planet separation for an orbit of a = 1 au at d_S = 10 pc is an order of magnitude smaller at 0.1''. Designing a coronagraph to reach high contrast at small inner working angle is challenging (Spaan & Greenaway 2006; Mawet et al. 2012), and is typically investigated in the context of direct imaging of point sources. This is complicated by the extended solar disk and the minimal separation between the edge of the disk and Einstein ring. The angular separation between the edge of the solar disk and the center of the Einstein ring for a source with d_S >> d_L is (Hippke 2020)

$$\begin{eqnarray}&&\delta \theta \approx \sqrt{\displaystyle \frac{2{R}_{s}}{r}}-\displaystyle \frac{{R}_{\odot }}{r},\end{eqnarray} \tag{ 49 }$$

which reaches a maximum value of ∼04 at a heliocentric distance of r ∼ 2300 au independently of wavelength. Specialized coronagraph designs for this particular case have only recently been investigated. An original suggestion from Turyshev and Toth was the use of an annular coronagraph (Turyshev & Toth 2020), which would only allow light in from the region of the Einstein ring while blocking light from the solar disk and outer corona. However, as we have seen, the Einstein Ring itself contains contamination from the target planet's host star, as well as potentially other objects in the field and so this may be insufficient. Loutsenko & Yermolayeva (2021) recently suggested a design that is a combination of an annular mask and a bar mask, with the intention of obstructing light from both the solar disk as well as the magnified off-axis contribution from the target planet's host star. Their design claims to reach contrasts of 10⁻⁸ between the solar disk and the region of the Einstein Ring, with similar performance for off-axis sources behind the bar mask. However, neither of these coronagraphs on their own can mitigate contamination from the region of the Einstein ring itself, which includes light from the corona and background sources. A spectroscopic mitigation strategy similar to the method described in Ruffio et al. (2019) could model all sources of contamination in the final image, using spectral differences between noise sources to obtain high contrast in the absence of or in combination with a coronagraph.

4.4.2. Background and Foreground Mitigation

In the process of designing a mission to the SGL, one should be prepared to deal with the possibility of a background object contaminating the Einstein ring. The target planet's host star will certainly contaminate two points along the Einstein Ring, and a possible background star could contaminate two additional points in a similar fashion. Even contamination from nontarget planets in the target solar system may be relevant if the angular separation between the target and contaminating planet is small enough. Selecting a target planet from a list of multiple potential target planets being considered should choose the safer option, which is farther from the galactic center. This will lower the probability of a background star, but the probability cannot be eliminated entirely. Early and deep observations of the near target region could provide advanced evidence for a clear background, but a robust design should include a contingency plan for mitigating this contamination if the observation is unlucky, as galactic proper motion could place a background object in the FOV sometime between launch and observation. By spatially separating the light in the Einstein ring, either a hardware or a software mask could be applied to the region around the points of contamination; although, scattered light from this unobstructable source will ultimately limit the possible signal-to-noise ratio (S/N) obtainable in uncontaminated pixels.

For proper mitigation of the corona and zodiacal light as astrophysical foregrounds, it will be necessary to model the temporal and azimuthal dependence of the coronal brightness to extract the full information contained in the Einstein ring. One observational strategy to achieve this is a differential measurement. A single telescope equipped with an integral field spectrograph could make spatially and spectrally resolved measurements of the Einstein ring and the corona as well as the region of the corona nearby the Einstein ring. By measuring the coronal brightness in the vicinity and extrapolating a model into the region of overlap, a differential measurement may be possible as a mitigation strategy. Another type of instrument design could achieve a similar measurement instead using two separate spacecraft: one being the primary observation telescope at the focal plane of the SGL and the other being an additional companion telescope observing the solar corona from a closer heliocentric distance, where the Einstein ring cannot form at d_L ≲ 550 au. Simultaneous observations could enable a differential measurement of the brightness of the coronal region with and without the Einstein ring present, and would not rely upon fitting a nearby coronal region and extrapolating. However, this concept would suffer from complexities involved with coordinating observations on different telescopes separated by a great distance.

Another alternative strategy for mitigation of the corona is by sending the observation craft to a significantly larger heliocentric distance of d_L = 2300 au for observing. Both the angular size of the Einstein ring and the solar disk shrink at larger d_L , but the Einstein ring shrinks proportionally to ${d}_{L}^{-1/2}$ while the solar disk shrinks proportionally to ${d}_{L}^{-1}$. This means at the optimal observation distance d_L = 2300 au, the angular separation between the edge of the solar disk and the center of the Einstein ring is maximal, and the foreground contamination will be minimal as the Einstein ring overlaps with the much fainter F-corona rather than the bright K-corona at small separations.

The solar foreground curves shown in Figure 19 are a combination of the coronal model fit in Section 4.3 and a fifth-order polynomial limb-darkening model for the brightness of the solar disk described in Hestroffer & Magnan (1998) and using data from Pierce & Slaughter (1977) and Neckel & Labs (1994). The data correspond to measurements at λ ∼ 0.58 μm, but additional measurements at wavelengths λ ∈ [0.78, 2.44] μm are reported in Figure 6 of Pierce (1954), which demonstrates the increasing effects of limb darkening at longer wavelengths. The specific intensity of the Earth B_Earth is taken from the model in Robinson et al. (2011), and the specific intensity of the Sun is taken from the blackbody model described in Section 4.1. The specific intensity ratio B_Earth/B_⊙ ∼ 10⁻⁶ is not the typical flux contrast ratio of the Earth and Sun ∼10⁻¹⁰, as it lacks the additional factor of ${({R}_{\mathrm{Earth}}/{R}_{\odot })}^{2}\sim {10}^{-4}$ to account for the difference in the solid angle of the bodies. A version of Liouville's theorem of conservation of phase space volume applied to gravitational lenses (Misner et al. 1973) proves that gravitational lenses preserves surface brightness or specific intensity up to a frequency shift due to gravitational redshift, while the magnification acts to modify the apparent solid angle of the object. This justifies our comparison of the specific intensities of the background and foreground along the line of sight, and demonstrates the advantage of observing at a farther heliocentric distance.

Zoom In Zoom Out Reset image size