∎ 11institutetext: Calvin Leung 22institutetext: MIT Kavli Institute for
Astrophysics and Space Research, Massachusetts Institute of Technology, 77
Massachusetts Ave, Cambridge, MA 02139, USA

# Wave Mechanics, Interference, and Decoherence in Strong Gravitational
Lensing

Calvin Leung  Dylan Jow  Prasenjit Saha  Liang Dai  Masamune Oguri  Léon V. E.
Koopmans

###### Abstract

Wave-mechanical effects in gravitational lensing have long been predicted, and
with the discovery of populations of compact transients such as gravitational
wave events and fast radio bursts, may soon be observed. We present an
observer’s review of the relevant theory underlying wave-mechanical effects in
gravitational lensing. Starting from the curved-spacetime scalar wave
equation, we derive the Fresnel-Kirchoff diffraction integral, and analyze it
in the eikonal and wave optics regimes. We answer the question of what makes
interference effects observable in some systems but not in others, and how
interference effects allow for complementary information to be extracted from
lensing systems as compared to traditional measurements. We end by discussing
how diffraction effects affect optical depth forecasts and lensing near
caustics, and how compact, low-frequency transients like gravitational waves
and fast radio bursts provide promising paths to open up the frontier of
interferometric gravitational lensing.

###### Keywords:

Gravitational lensing, wave optics, gravitational waves, transients

###### Contents

  1. 1 Introduction
  2. 2 The curved-spacetime scalar wave equation
  3. 3 Different Regimes in Wave Optical Gravitational Lensing
    1. 3.1 Beyond Scalar Wave Optics
  4. 4 Eikonal Optics
  5. 5 Diffractive optics ($\Omega\gtrsim 1$)
    1. 5.1 A useful analytical approximation
  6. 6 Angular broadening/finite-size effects
  7. 7 Exploiting Chromaticity in Lensing
  8. 8 Finite magnifications near the Einstein ring
  9. 9 Finite magnifications near fold caustics
  10. 10 Modified lensing probabilities in wave optics
  11. 11 Observational prospects and discussion

##  1 Introduction

Gravitational lensing in the geometric optics approximation has enjoyed
immense success as a unique and precise probe of cosmology and astrophysics.
However, electromagnetic and gravitational radiation both obey wave equations,
have a wave-like nature, and interfere and diffract as they propagate to Earth
from the cosmos. Just how the first theoretical predictions of gravitational
lensing (Chwolson, 1924) were made decades before the first observation of an
Einstein ring (Hewitt et al., 1988), wave effects in gravitational lensing
were predicted long before there was any hope of detecting them.

Now, there is hope on multiple fronts: Lensed supernovae (Kelly et al., 2015)
have already been detected. In principle, wave-mechanical effects should be
observable in the lensing of gravitational wave (GW) signals from binary black
holes (Abbott et al., 2016), binary neutron stars (Abbott et al., 2017), and
even black hole-neutron star mergers (Abbott et al., 2021a), as well as a
large population of fast radio bursts (FRBs) at low frequencies (CHIME/FRB
Collaboration et al., 2021). Growing observational efforts in these areas puts
astronomers well on their way to detecting interference effects in
gravitational lensing. Measuring these interference effects enables unique
tests of gravity: for example, testing the equivalence principle tests using
precise measurements of time lags between electromagnetic and gravitational-
wave images (Takahashi, 2017). It is also of interest from the perspective of
testing the foundations of quantum mechanics (Miller & Wheeler, 1997; Jacques
et al., 2007; Doyle & Carico, 2009; Leung et al., 2018).

Wave-mechanical effects in gravitational lensing also offer a significant
enhancement of gravitational lensing as an astrophysical tool. In traditional
lensing, only two numbers (the flux magnification and time delay) can be
measured, leading to difficult degeneracies in lensing-based inference: the
mass-sheet degeneracy is perhaps the best-known example (Falco et al., 1985;
Saha, 2000). In the wave-optics limit, much more information is present, since
the flux magnification, and potentially the phase of the wave, may both be
measured as a function of wavelength. In certain circumstances, time delays
may also be extracted interferometrically with dramatically-increased
precision. As an example, lensing of FRBs can allow for interferometric
measurements of lensing time delays. These time delays can be measured with a
precision of $\sim 10^{-9}-10^{-6}$ seconds: several orders of magnitude more
precise than any lensing delay ($\sim$ 1 millisecond) inferred from measured
FRB light curves alone (Wucknitz et al., 2021; Jow et al., 2020).

Wave-mechanical effects may also shed light on the emission physics of compact
objects (such as GWs and FRBs), by using gravitational lenses as
astrophysical-scale interferometers. This concept has been demonstrated in the
study of pulsar scattering in the ionized interstellar medium to resolve the
nanoarcsecond-scale emission region of the Crab pulsar (Lin et al., 2022). The
observation of interference effects in gravitational lensing is a new frontier
in using lensing to learn about our universe.

The theory of interference effects in gravitational lensing was developed
since as far back as 1973 (Ohanian, 1973, 1974; Peters, 1974). Further
developments have mainly focused on GW lensing (Deguchi & Watson, 1986a, b;
Nakamura, 1998; Nakamura & Deguchi, 1999; Baraldo et al., 1999; Takahashi &
Nakamura, 2003; Macquart, 2004) and because of the relative simplicity of the
gravitational waveforms produced by compact binary coalescences, and the low
(kilohertz) emission frequencies at which interference effects are most
pronounced. However, wave-mechanical effects in lensing are of growing
importance at frequencies far higher than gravitational waves. An early
example was in using “femtolensing” to probe extremely low mass scales by
analyzing the spectra of gamma ray bursts (GRBs) (Gould, 1992; Stanek et al.,
1993; Nemiroff et al., 1993) More recently, various other applications using
lensing to search for exotic phenomena have had to contend with wave-
mechanical effects, which limit the reach of the search. Some examples include
optical microlensing searches for primordial black holes (PBHs; see e.g.
Sugiyama et al. (2020), searches for PBHs using lensed FRBs (Katz et al.,
2020; Kader et al., 2022; Leung et al., 2022), searches for small-scale
density perturbations in minihalos or e.g. axion dark matter (Venumadhav et
al., 2017; Dai & Miralda-Escudé, 2020), and even searching for lensing by
cosmic strings (Xiao et al., 2022).

In this review, we aim to present the fundamental concepts underlying wave
interference in gravitational lensing. We derive the Fresnel-Kirchhoff
integral from first principles (Sec. 2), and analyze it in two opposite
regimes. We first analyze the stationary-phase (a.k.a. Eikonal) regime in Sec.
4. Then, we present some analytical solutions from the diffractive (a.k.a.
wave optics) regime (Sec. 5) and comment on numerical solution methods for
this regime, in which the saddle-point approximation fails. We discuss how the
finite angular sizes of sources pose a challenge to detecting interference
effects (Sec. 6).

In anticipation of detections of wave-mechanical effects, we discuss several
observationally-relevant features of wave mechanics in lensing. First, we
discuss how the chromaticity of lensing in the wave optics regime can be
observed and used to significantly increase the information extracted out of
systems which exhibit wave optics effects (Sec. 7). We then discuss the effect
of diffraction near caustics (Sec. 8, 9) and talk about how wave optics alters
calculations of the lensing optical depth (Sec. 10) in searches for exotic
objects. We finally comment on how various compact transients offer
opportunities for detections in the future (Sec. 11).

##  2 The curved-spacetime scalar wave equation

Much of gravitational lensing is captured by the scalar wave equation in a
curved spacetime:

$\displaystyle\qquad\qquad\qquad\qquad 0$
$\displaystyle=g^{ab}\nabla_{a}\nabla_{b}\phi.$ (1)

Scalar wave theory captures interference and diffraction phenomena in
gravitational lensing 111but ignores e.g. polarization-dependent effects., and
is sufficient for analyzing most interference phenomena in gravitational
lensing. The starting point of scalar wave theory is the Fresnel-Kirchoff
integral, from which most lensing phenomena can be derived. It requires a
little bit of general relativity, but we will assume the reader has little
background in that area.

The formalism of gravitational lensing is similar to the Newtonian dynamics of
point masses orbiting in a potential. In the same way that the external
potential is fixed for point masses, we consider the dynamics of a scalar
field representing our wave, $\phi$, on a static background spacetime, whose
curvature is specified by the metric tensor $g^{ab}$. In a weak gravitational
field 222Confusingly, “strong lensing”, which refers to the regime in lensing
in which multiple images arise, occurs mostly in the weak-field limit; we will
discuss the strong-field limit briefly in Sec. 3.1 with only scalar degrees of
freedom, the metric tensor can be written in terms of the gravitational
potential as $g_{00}=-(1+2U)$ and $g_{ii}=1-2U$ for $i=1,2,3$, and where $U$
is the Newtonian gravitational potential, and where we use units in which
$c=1$ for now.

In Eq. 1, $\nabla_{a}$ refers to the covariant derivative. The covariant
derivative needs to be expanded into partial differential equations using
several identities from general relativity. We introduce the notation
$g=\det(g^{ab})=g^{00}g^{11}g^{22}g^{33}$ for our diagonal metric. Since
scalars are Lorentz-invariant, the covariant derivative of a scalar is the
same as its partial derivative:

$\displaystyle\qquad\qquad\qquad\qquad\nabla_{b}\phi$
$\displaystyle=\partial_{b}\phi.$ (2)

The analogous rule for (the covariant divergence of) a vector is333This in
turn arises from an identity for a contracted Christoffel symbol:
$\Gamma^{c}_{ca}=\dfrac{\partial_{a}(\sqrt{-g})}{\sqrt{-g}}$; see e.g.
Lightman et al. (1975) (Exercise 7.7j) for more details.

$\displaystyle\qquad\qquad\qquad\qquad\nabla_{a}V^{a}$
$\displaystyle=\dfrac{1}{\sqrt{-g}}\partial_{a}(\sqrt{-g}V^{a}).$ (3)

Eq. 3 implies

$\displaystyle g^{ab}\nabla_{a}\nabla_{b}\phi$
$\displaystyle=\nabla_{a}(g^{ab}\nabla_{b}\phi)=\nabla_{a}(g^{ab}\partial_{b}\phi)$
(4)
$\displaystyle=\dfrac{1}{\sqrt{-g}}\partial_{a}(\sqrt{-g}g^{ab}\partial_{b}\phi)=0.$
(5) Eq. 5 is the wave equation often used as a starting point in studying
gravitational lensing in wave optics, e.g. Nakamura & Deguchi (1999); Macquart
(2004). We now expand the four non-zero terms in $U$, since we are working the
weak-field limit, $|U|\ll 1$; this justifies neglecting terms of
$\mathcal{O}(U^{2})$. Since the potential is static, the 00 term evaluates to
$\displaystyle\dfrac{1}{\sqrt{-g}}\partial_{0}(\sqrt{-g}g^{00}\partial_{0}\phi)$
$\displaystyle=-\dfrac{1}{1+2U}\partial^{2}_{t}\phi\approx(-1+2U)\partial^{2}_{t}\phi.$
(6) The $a=b=1$ term is
$\displaystyle\dfrac{1}{\sqrt{-g}}\partial_{1}(\sqrt{-g}g^{11}\partial_{1}\phi)$
$\displaystyle=\phi^{\prime\prime}+(-4\phi^{\prime}U^{\prime}+2\phi^{\prime\prime})U.$
(7)

where the primes denote $\partial_{1}$ (the $a=b=2$ and $a=b=3$ terms are
similar). From dimensional arguments, $U^{\prime}\sim|U|/R_{s}$ (the curvature
scale of the spacetime) and the scale of $\phi^{\prime}\sim 1/\lambda$, and
$\phi^{\prime\prime}\sim 1/\lambda^{2}$. We are justified in neglecting the
$\phi^{\prime}U^{\prime}$ term of Eq. 7, which is smaller than the other two
by a factor of $\mathcal{O}(\lambda/R_{s})$444This is called the Eikonal
approximation, from which we derive the equations of diffractive lensing. This
should not be confused with the Eikonal / stationary-phase limit of the
Fresnel-Kirchhoff diffraction integral.. The full equation is

$\displaystyle 0$
$\displaystyle=(-1+2U)\partial_{0}^{2}\phi-(1+2U)\sum_{i=1,2,3}\partial_{i}^{2}\phi+\mathcal{O}(U^{2})$
(8) which can compactly be written as $\displaystyle 0$
$\displaystyle=\nabla^{2}\phi-(1-4U)\partial^{2}_{t}\phi$ (9) where the
$\nabla$ (without the indices) refers to the flat-space Laplacian. We now have
a standard PDE for which various solution approaches exist. For a point-mass
(or a point charge (Davydov, 2013)), one solution taken by Peters (1974);
Deguchi & Watson (1986b) is to substitute an ansatz of the form
$\phi(\vec{r},t)=\tilde{\phi}(\vec{r})\exp(i\omega t)$. This gives us a
Schrödinger-like equation for a free particle in a Coulomb potential
$V=4\omega U$ and whose “mass” is $\omega/2$:
$\displaystyle-i\partial_{t}\phi$
$\displaystyle=-\dfrac{1}{\omega}\nabla^{2}\phi+4\omega U\phi.$ (10) We take
an alternate approach using the Feynman path integral to derive the Fresnel-
Kirchhoff diffraction integral by summing over all possible paths through the
lens plane (Nakamura & Deguchi, 1999; Yamamoto, 2003). Further simplifying Eq.
10 gives another commonly-used starting point: $\displaystyle
4\omega^{2}U\tilde{\phi}$ $\displaystyle=(\nabla^{2}+\omega^{2})\tilde{\phi}.$
(11) To solve this, we use the ansatz
$\tilde{\phi}(\vec{r})=F(\vec{r})\exp(ikr)/r$ with Eq. 10, giving
$\displaystyle 4\omega^{2}UF$
$\displaystyle=-k^{2}F+2ik\partial_{r}F+\partial^{2}_{r}F+\dfrac{1}{r^{2}}\nabla_{\theta}^{2}F+\omega^{2}F.~{}$
(12) In this equation, we have decomposed the spherical Laplacian into its
radial and angular parts: $\displaystyle\nabla^{2}$
$\displaystyle=\dfrac{1}{r^{2}}\dfrac{\partial}{\partial
r}\left(r^{2}\dfrac{\partial}{\partial r}\right)+\nabla_{\theta}^{2}$ (13)
Note that in the un-lensed case ($U=0$), $F=1$ solves Eq. 12 (assuming
$\omega=k$). The ansatz also captures the free wave intensity decaying as
$r^{-2}$. Therefore, to an overall normalization constant which we define
later, we expect $F$ to capture the interference pattern generated by lensing.
$F(\vec{r})$ is called the “amplification factor” of the wave amplitude (not
the flux). The characteristic transverse scale over $F$ varies is estimated by
considering a spherical wave of wavelength $\lambda$ emanating from a point
towards a plane at some distance $D\gg\lambda$. In the tangential (transverse)
direction, $F$ varies over a characteristic scale $R_{F}=\sqrt{D\lambda}$ (the
Fresnel scale). Therefore $\partial_{r}^{2}F$ is smaller than
$2ik\partial_{r}F$ by a factor of $\mathcal{O}(\sqrt{\lambda/D})$ and can be
neglected. Finally, using $\omega=k$ yields Eq. 14, which can be thought of as
the “flattened” version of Eq. 10, where the propagation (radial) direction is
collapsed. In the flat-sky approximation, we can work in transverse
coordinates $\bm{q}=r\bm{\theta}$, with $\bm{\dot{q}}=r\bm{\dot{\theta}}$, and
$\nabla^{2}_{q}=\dfrac{1}{r^{2}}\nabla^{2}_{\theta}$. $\displaystyle
i\partial_{r}F$ $\displaystyle=-\dfrac{1}{2\omega}\nabla^{2}_{q}F+2\omega UF$
(14) Since Eq. 14 also resembles a Schrödinger equation, we may solve it using
the Feynman path integral. The classical Hamiltonian corresponding to Eq. 14
governs transverse deflections of rays evolving over the propagation path; the
radius $r$ from the waves’ origin plays the role of time in the classical
analogy. Hence, the evolution of the transverse trajectories of rays is
specified by $\bm{q}$ and $\bm{\dot{q}}$. The conjugate momentum is
$\bm{p}=mr\dot{\bm{\theta}}$, where the analogous “mass”555not to be confused
with the lens mass! $m=\omega$ and the classical potential is $V=2\omega U$.
The Lagrangian corresponding to this Hamiltonian can be obtained by a Legendre
transform: $\displaystyle L({\bm{q}},{\bm{\dot{q}}},r)$
$\displaystyle={\bm{\dot{q}}}\cdot(m\dot{\bm{q}})-H(\bm{q},\bm{p})=\omega\left[\dfrac{{\bm{\dot{q}}}}{2}-2U\right]$
(15) where $H=\dfrac{\bm{p}^{2}}{2m}+V(r,\bm{q})$. Solutions to Eq. 14 can be
written as an integral over all possible ray trajectories $\bm{q}(r)$:
$\displaystyle F(\omega,\vec{r})$
$\displaystyle=\int\mathcal{D}[\bm{q}(r)]\exp\left[i\int
L(\bm{q},\bm{\dot{q}},r)~{}dr\right]$ (16)

The kinetic and potential terms in the Lagrangian correspond to the geometric
time delay ($\tau_{\mathrm{geo}}$) and the Shapiro time delay
($\tau_{\mathrm{grav}}$) in optics(Blandford & Narayan, 1986). Putting the
factors of $c$ back in, the geometric term is

$\displaystyle\hat{\tau}_{\mathrm{geo}}$
$\displaystyle=\int\dfrac{r^{2}\bm{\theta}(r)^{2}}{2}~{}dr\approx\dfrac{D_{L}D_{S}}{2D_{LS}c}|\bm{\theta}_{l}-\bm{\beta}|^{2}.$
(17)

where $\bm{\theta}_{l}$ and $\bm{\beta}$ are angles from the observer’s
perspective, and where we have approximated the paths as straight rays bent
“instantly” once a ray hits the lens plane (see Fig. 1 in (Feldbrugge & Turok,
2020)). This is the same as assuming that far away from the lens plane, the
index of refraction of the medium is homogeneous enough such that rays travel
in approximately straight lines.
The gravitational term is

$\displaystyle\hat{\tau}_{\mathrm{grav}}$ $\displaystyle=-\dfrac{2}{c^{3}}\int
U(\bm{q},r)~{}dr.$ (18) This equation integrates the gravitational potential
along the line of sight in the source frame. In the lens frame, we approximate
this integral by collapsing the gravitational potential over the $z$ axis (the
thin-lens approximation). Defining the physical impact parameter $\bm{b}$
($=r_{l}\bm{\theta}$ from the source’s perspective, and
$\bm{b}=D_{L}\bm{\theta}_{l}$ from the observer’s perspective) we define
$\displaystyle\hat{\psi}(\bm{b})$ $\displaystyle=\dfrac{2}{c^{3}}\int
U(\bm{b}+z\bm{\hat{z}})~{}dz\approx-\hat{\tau}_{\mathrm{grav}}.$ (19)

The approximations in Eq. 17 and 19 reduce the infinite-dimensional Feynman
path integral to the Fresnel-Kirchhoff diffraction integral:

$F(\omega,\bm{\beta})\propto\int
d^{2}\bm{\theta_{l}}\exp\left\\{i\omega\left[\dfrac{D_{L}D_{S}}{2D_{LS}c}|\bm{\beta}-\bm{\theta}_{l}|^{2}-\hat{\psi}(D_{L}\bm{\theta}_{l})\right]\right\\}~{}$
(20)

where $\omega$ is the observing frequency and $\bm{\beta}$ is the source’s
angular position with respect to the lens. In a cosmological context, the
distances thus far should be replaced by angular diameter distances, and
$\omega\to\omega(1+z_{l})$.

At this point it is convenient to define a characteristic scale to measure
deflection angles. While any scale may be chosen, each lens model comes with a
different “natural” choice of the angular scale. For instance, the angular
scale is often chosen such that the dimensionless angular radius of the
Einstein ring is 1. For example, for a singular isothermal sphere
characterized by its velocity dispersion $\sigma$, the deflection angle is a
constant and defines the angular scale $\theta_{E}$ by which to normalize Eq.
20:

$\theta_{E}=4\pi\dfrac{\sigma^{2}}{c^{2}}\dfrac{D_{LS}}{D_{S}}$ (21)

While the singular isothermal sphere is common in geometric optics, its
solution in wave optics is quite complicated. Instead, in wave optics, the
point mass model is much simpler and more commonly analyzed. We write Eq. 20
in terms of dimensionless variables using a notional point-mass lens model.

$\displaystyle\theta_{E}$
$\displaystyle=\sqrt{\dfrac{4GM}{c^{2}}\dfrac{D_{LS}}{D_{L}D_{S}}}$ (22)
$\displaystyle\bm{\theta}_{l}$ $\displaystyle\to\mathbf{x}\theta_{E}$ (23)
$\displaystyle\bm{\beta}$ $\displaystyle\to\mathbf{y}\theta_{E}$ (24)
$\displaystyle\hat{\psi}(\bm{\theta}_{l})$
$\displaystyle\to\psi(\bm{x})\theta_{E}^{2}\dfrac{D_{L}D_{S}}{D_{LS}}$ (25)
$\displaystyle\omega$
$\displaystyle\to\Omega\dfrac{cD_{LS}}{D_{L}D_{S}\theta_{E}^{2}}$ (26)

Note that

$\Omega=4GM\omega/c^{3}=4\pi R_{s}/\lambda$ (27)

has a physical interpretation as the answer to the question: how large is a
wavelength relative to the Schwarzschild radius, $R_{s}=2GM/c^{2}$? The
dimensionless frequency, $\Omega$, therefore plays a crucial role in
determining whether geometric optics is a valid approximation. Eq. 20 becomes

$\displaystyle F(\Omega,\bm{y})$ $\displaystyle=\dfrac{\Omega}{2\pi i}\int
d^{2}\bm{x}\exp(i\Omega\tau(\bm{x},\bm{y}))$ (28) where we have also
introduced a normalization constant into $F(\Omega,\bm{y})$ such that when
$\psi\to 0$, $F\to 1$, and where the dimensionless delay is
$\displaystyle\tau(\bm{x},\bm{y})$
$\displaystyle=\dfrac{|\bm{x}-\bm{y}|^{2}}{2}-\psi(\bm{x}).$ (29)

Eqs. 28 and 29 are valid in the limit of weak gravity and small angular
deflections. However, this non-dimensionalization scheme is not restricted to
a point-mass lens. The constant $M$ can be anything, but it is understood to
be some characteristic mass scale of the lens. Similarly $\theta_{E}$ need not
be the Einstein radius, but is understood to be some characteristic scale of
deflections.

##  3 Different Regimes in Wave Optical Gravitational Lensing

Eq. 28 looks simple, but describes extremely diverse phenomena, including all
of standard (geometric optics-based) gravitational lensing phenomena. This is
because the complex-valued amplification $F$ (not to be confused for the
positive-valued flux magnification from geometric optics) carries both
amplitude and phase information, the latter of which can result in
constructive and destructive interference. For an amplification factor of $F$,
the flux of the source gets magnified by $|F|^{2}$, in analogy with the Born
rule from quantum mechanics. To extend the analogy, the overall phase of $F$
is arbitrary, as the choice of the origin of the source plane $\bm{y}$ was
arbitrary. However, we will see later that relative phase between the
different contributions to $F$ (the Morse phase) can alter the observables.

Just like in traditional geometric optics, we care about the image positions,
defined as stationary points of $\tau(\bm{x},\bm{y})$ obtained by solving

$\dfrac{\partial\tau}{\partial{\bm{x}}}=0.$ (30)

However, in wave optics, we sum over all possible paths through the lens
plane: the resulting amplification is proportional to the size of the coherent
“patches,” or Fresnel zones, surrounding the stationary points in the lens
plane. The Fresnel zone is defined here as the patch surrounding each
stationary point in the lens plane for which the phase has not changed by
2$\pi$.

Figure 1: Top row: The phase $\Omega T(\bm{x})$ corresponding to the lensing
potential $T(\bm{x},\bm{y})=\dfrac{1}{2}|\bm{x}-\bm{y}|^{2}-\ln(|x|)$, for
three values of $\Omega=1$ (heavily diffracted), $\Omega=5$, $\Omega=25$. We
fix the source at $\bm{y}=(1,0)$ and evaluate the phase as a function of
position in the lens plane, visualized with a $2\pi$-periodic colormap. To
guide the eye, we have drawn two circles centered on the loci of the image
positions $\bm{x}_{\pm}$ evaluated in the geometric optics limit. The radii of
the pairs of circles are proportional to the flux magnification ratios
$\mu_{\pm}$, and are scaled as $1/\Omega$. Top left: In the heavily-diffracted
regime ($\Omega=1$), the time delay between images is less than $2\pi$
radians. The images are superimposed in the time domain within one wavelength,
and are rendered indistinguishable due to diffraction. Top center: As the
frequency increases ($\Omega=5$), the images are separated by more than one
wavelength, and there are several oscillations between the two stationary
points of the lens. The images are therefore no longer blended by diffraction;
they become distinguishable in both the time domain and the angular domain.
The interference leads to constructive and destructive interference of the
phase of the integrand when integrated over the whole lens plane, and shows a
simple sinusoidal functional form when plotted against $\Omega$. Top right: In
the high-frequency limit (here shown as $\Omega=25$ for visualization
purposes), the two images in the time domain are very well-resolved. The total
magnification may be treated as a sum of discrete images – the discrete
stationary points of the lens potential – using geometric optics.
Middle/bottom rows: The saddle-point approximation to the “plus”/“minus” image
individually, for the same three values of $\Omega$. It can be seen that the
saddle-point approximation works better and better for higher values of
$\Omega$, as the relevant parts of the integrand become increasingly confined
around the stationary points. Note that a second saddle point appears on the
right side of the bottom right panel; this is a image rendering artifact.
Figure 2: The transition from geometric to diffractive optics, as quantified
by $\Omega=4\pi R_{s}/\lambda$, as function of lens mass and observing
wavelength. Right half: At longer wavelengths, advances in instrumentation
(e.g. wideband voltage recording in radio telescopes and gravitational wave
detectors) have enabled measurements of the amplitude (coherent detection)
instead of the flux (incoherent detection). Left half: At high $\Omega$, the
stationary-phase approximation holds, and geometric/Eikonal optics applies. In
most scenarios, this is an excellent approximation. However, the stationary-
phase approximation breaks down for long wavelengths or shallow local minima
in the time-delay potential, ($\Omega\ll 1$). In this regime lensing is
referred to as diffractive, and the unique frequency-dependence of diffractive
lensing may be observable. In the $\Omega\ll 1$ regime, the lenses can be
thought of as being “too small” to impart any phase on the passing wave.

In the top row of Fig. 1 we illustrate the phase delay potential
$\Omega\tau(\bm{x},\bm{y})$ modulo $2\pi$ over the lens plane for several
values of $\Omega$. We hold a source fixed at $\bm{y}_{0}=(1,0)$ (one Einstein
radius away from the origin of the source plane) for various choices of
$\Omega$, with a point-mass lens $(\psi(\bm{x})=\log(|\bm{x}|)$. To compare
the coherent patches with traditional geometric optics, we place dotted
circles centered at the two stationary points of $\tau$, at
$\bm{x}=(x_{\pm},0)$ for

$x_{\pm}=\dfrac{y\pm\sqrt{y^{2}+4}}{2}.$ (31)

The radii of the circles are chosen according to their geometric-optics
magnifications, and are scaled as $1/\Omega$ to ensure that at high $\Omega$,
due to the prefactor in Eq. 28, the flux magnification $|F|^{2}$ stays
constant.

The phase over the lens plane varies significantly as a function of $\Omega$,
suggesting that the amplification, and therefore the lens magnification, can
become frequency-dependent, unlike in traditional gravitational lensing, where
magnification ratios are completely achromatic.

This frequency dependence arises because Eq. 27 introduces a characteristic
length scale (the wave period), and therefore chromaticity, into gravitational
lensing. Lensing phenomenology can be classified depending on whether the
integral in Eq. 28 can be treated using the saddle-point approximation near
its stationary points.

The wisdom is that $\Omega$ is the most relevant quantity; however we shall
see later (in Eq. 33) that it is really the product $\Omega A$ which sets the
size of a coherent patch in the lens plane around the individual stationary
points of $\tau(\bm{x},\bm{y})$ (Jow et al., 2022; Berry, 2021), where $A$ is
the determinant of the magnification matrix:

$A_{ab}=\dfrac{\partial^{2}\tau}{\partial\bm{x}_{a}\partial\bm{x}_{b}}=1-\dfrac{\partial^{2}\psi}{\partial\bm{x}_{a}\partial\bm{x}_{b}}.$
(32)

When the phase of the integrand varies rapidly (large values of $\Omega A$),
the individual stationary points are separable into discrete images on the sky
(and into discrete events in the time domain, for transient events). The first
Fresnel zone of each image is observationally separable from the others, and
they interfere with each other to form an oscillatory fringe pattern. This is
called the eikonal limit. When the observing frequency is so high (extremely
large values of $\Omega A$) that the oscillations are washed out by summing
over the finite extent of the source in the source plane, ordinary geometric
optics applies. This is the regime in which most gravitational lensing
observations have been conducted.

In this work, we make the distinction between the Eikonal limit (i.e. the
stationary-phase approximation of Eq. 28) and geometric optics (which we take
to be the regime in which the Eikonal images are added together incoherently
at the observer). Some authors do not make this distinction and simply refer
to the Eikonal limit as geometric optics as well. Grillo & Cordes (2018)
distinguish the incoherent and coherent regimes as “zeroth" and “first-order"
geometric optics. We also note that the Eikonal limit of optics that we refer
to in this section is distinct from the Eikonal approximation used in Sec. 2
in order to go from Eq. 5 to Eq.9.

For small values of $\Omega A$, the first Fresnel zones for each image blend
together. In this regime the images are inseparably blurred together on the
sky, and are separated by less than a period in the time domain. There are
still multiple stationary points of the time-delay potential, so the lensing
is still “strong lensing”, but it is impossible to associate a particular
photon to a particular image. For this reason, it is helpful to think of this
regime as a perturbation of the unlensed ($F=1$) case. We will see that the
lensing is indeed suppressed in this limit.

To organize these limits, we consider possible values of $\Omega$ and
$\kappa$. With modern observational capabilities ranging from ground-based
(and soon, space-based) gravitational wave detectors to gamma-ray telescopes,
$\Omega$ can vary by 30 orders of magnitude for different lenses and
wavelengths, and $A$ can vary from order unity to $\sim 10^{3}$ in rare
extreme-magnification scenarios. Roughly speaking, we can focus on $\Omega$,
keeping in mind that $A$ may play an important role under circumstances where
the stars align. Fig. 2 shows $\Omega$ as a function of the mass scale of the
lens and the wavelength of the source, while Fig. 3 schematically depicts the
spectrum of a lensed source as a function of $\Omega$.

Figure 3: A schematic depicting different regimes in gravitational lensing as
a function of $\Omega$. Left red: When $\Omega\sim 1$, we are in the
diffractive or wave optics regime where the lensing is present but suppressed
($|F|^{2}\to 1$) for $\Omega<<1$. Center blue: In the eikonal regime, the
interference fringe between the two images, here represented as two amplitudes
$F_{+}$ and $F_{-}$, is clearly visible. The total magnification has a perfect
sinusoidal dependence (see Eq. 38) with a characteristic spectral oscillation
scale of $1/\tau_{lens}$. Right black: Finally, at high frequencies, we enter
the geometric optics (traditional lensing) regime when the lens’s resolving
power is sufficiently great to resolve the finite extent of the source. In
this limit, the interference effects become washed out, and the lensing
magnification becomes achromatic. In this limit the total flux magnification
is the sum of the image flux magnifications: $|F_{+}|^{2}+|F_{-}|^{2}$.

Within Fig. 2, we have also delineated a second boundary, which denotes
whether the radiation is detected incoherently (without phase information,
e.g. with a flux detector) or coherently (with phase information, e.g. with an
antenna). At long wavelengths (radio and gravitational waves), it is possible
with current technology to detect the phase of the radiation field using
electromagnetic/gravitational antennas. These situations where we have access
to the phase of the electromagnetic/gravitational field are extremely powerful
and enable qualitively different observations to be carried out. For example,
recording the field amplitude instead of the intensity (voltages instead of
light-curves) allows antennas to reach Nyquist-limited time resolution for
resolving shorter lensing delays, regardless of the flux variability timescale
of the source. This is similar to radio interferometry, where coherent voltage
recording allows precise time-delay measurements between different receivers
even when the sources do not exhibit time variability in their light curves.
Improving time resolution through coherent detectors and recording
instrumentation therefore opens up new phase space towards smaller time delays
(corresponding to smaller mass scales) in gravitational lensing.

###  3.1 Beyond Scalar Wave Optics

These aforementioned regimes in scalar wave optics cover much of the relevant
phenomenology of interference effects in gravitational lensing. From an
observational standpoint, the eikonal and wave optics regimes are still
uncharted territory. Virtually all observations of gravitational lensing thus
far fall squarely in the geometric optics regime.

However, in an era of gravitational-wave observations from the ground and from
space, there are compelling reasons to push further theoretically. Notably, we
have worked with scalar waves in the weak gravity limit, ignoring cases where
$U\sim 1$, for which Eq. 20 is not valid. Numerically simulating wave
propagation in strongly-curved spacetimes becomes necessary and is an active
area of research in general relativity (Zenginoǧlu & Galley, 2012; Yang et
al., 2014). These numerical solutions of scalar wave propagation reveal rich
phenomenology, particularly near caustics. To give an idea of the
phenomenology, Zenginoǧlu & Galley (2012) computes the Green’s function of Eq.
1 near the Schwarzchild radius of the lens, and observe multipath propagation
for integer numbers of cycles around the photon ring. The observed “caustic
echo” has a period in good agreement with the light-crossing time of the
photon ring (Stewart, 1989).

Other aspects of lensing are not captured by the scalar theory, such as
polarization-dependent effects. Two examples from general relativity are the
gravitational Faraday rotation induced by rotating gravitational potentials
(Ishihara et al., 1988; Nouri-Zonoz, 1999; Asada & Kasai, 2000; Li et al.,
2022) and the gravitational spin-Hall effect, in which different polarizations
of light traverse different trajectories through a gravitational potential
(Gosselin et al., 2007; Oancea et al., 2020).

Note also that while the scalar wave formalism for both gravitational and
plasma lensing are effectively identical (the only difference being that
gravitational potentials are frequency independent, but plasma lensing
potentials have a $\sim f^{-2}$ dependence), magnetic fields in ionized plasma
cause different polarizations to experience different refractive indices in a
phenomenon known as birefringence (Li et al., 2019). Thus, while the scalar
wave theory of gravitational lensing and plasma lensing share significant
overlap, notable differences arise in the vector propagation.

A final limitation is that we have so far limited our discussion to the
effects of scalar perturbations to the metric on the propagation of waves
through space-time. In general, metric perturbations may themselves have
vector and tensor components (e.g. gravitational waves). Gravitational waves
will induce additional time delays in the propagation of light from background
sources which may be observable (indeed, this is the basis for pulsar timing
array experiments). Gravitational waves themselves may therefore lens coherent
sources of radiation (Boyle & Pen, 2012; Rahvar, 2018) inducing characteristic
interference patterns; this is an under-explored and active area of research.
Polarization effects, such as the gravitational Faraday rotation, may also be
more pronounced in lensing by tensor and vector perturbations, and may have
observable impacts on CMB lensing (Dai, 2014). However, in light of the
complexity and magnitude of these effects, for the remainder of this review
paper, we will focus on scalar wave optics in the aforementioned regimes:
eikonal scalar optics, diffractive scalar optics, and the transition from
coherent to incoherent (traditional) gravitational lensing.

##  4 Eikonal Optics

Lensing in eikonal optics is similar to treating the source as a plane wave
hitting an opaque screen with a set of well-separaged pinholes of negligible
size. Each pinhole’s size and phase relative to the other pinholes is set by
the traditionally-calculated flux magifications and time delays. To see this
mathematically, we simplify Eq. 28 using the saddle-point approximation.
Conceptually, this collapses the integral over the lens plane to a coherent
superposition of a small number of paths defined by the stationary points of
the lens potential. Taking the limit as $\Omega\to\infty$ results in rapid
oscillations of the phase of the integrand as a function of $\bm{x}$. In this
limit, only the stationary points of $\Omega\tau(x,y)$ as a function of $x$,
which we index by $j$, contribute constructively to the integrand. These
positions $\bm{x}_{j}$ correspond to the image positions as calculated with
geometric optics. Near the $j^{th}$ stationary point, the time delay potential
can be approximated by a Taylor expansion:

$\displaystyle\tau_{j}(\bm{x},\bm{y})\approx\tau(\bm{x}_{j})$
$\displaystyle+\dfrac{1}{2}A^{j}_{ab}(\bm{x}-\bm{x}_{j})_{a}(\bm{x}-\bm{x}_{j})_{b}+\ldots$
(33) where $A^{j}_{ab}=\partial_{a}\partial_{b}\tau(\bm{x}_{j})$ is the
inverse magnification matrix evaluated at the $j$th image. In Fig. 1, we plot
in the top row the phase $\Omega\tau(\bm{x},\bm{y})$ modulo $2\pi$, and in the
bottom two rows the saddle-point approximation of the phase potential about
each of the two stationary points. We approximate the Fresnel-Kirchhoff
integral as the coherent superposition of both of these expansions. This is
called the stationary-phase (or saddle-point) approximation, which relies on
the fact that contributions far away from the stationary points cancel each
other out. $\displaystyle\int d^{2}\bm{x}\exp(i\Omega\tau(\bm{x}))$
$\displaystyle\approx\sum_{j}\int d^{2}\bm{x}\exp(i\Omega\tau_{j}(\bm{x}))$
(34) $\displaystyle=\sum_{j}\exp(i\Omega\tau(\bm{x}_{j}))$ $\displaystyle\int
d^{2}\bm{x}\exp(i\Omega
A^{j}_{ab}(\bm{x}-\bm{x}_{j})_{a}(\bm{x}-\bm{x}_{j})_{b}/2)$ (35) We evaluate
the above expression analytically, being a sum of complex Gaussian integrals.
Since $A^{j}_{ab}$ is a symmetric matrix, it is diagonalizable via a rotation
of $\bm{x}$; we call its eigenvalues $\lambda^{j}_{1}$ and $\lambda^{j}_{2}$.
The Gaussian integral associated with each image contributes a factor of
$\dfrac{2\pi i}{\Omega\sqrt{\lambda_{1}^{j}}\sqrt{\lambda^{j}_{2}}}$. The
eigenvalues $\lambda^{j}_{1},\lambda^{j}_{2}$ may be both positive if the
image is a minimum of the time-delay potential, both negative if the image is
a maximum of the time-delay potential, or they may have different signs if the
image is at a saddle point of the potential. If we define $F_{j}$ by
$\displaystyle F_{j}^{-1}$
$\displaystyle=\sqrt{\lambda_{1}^{j}}\sqrt{\lambda^{j}_{2}}$ (36) the
amplification becomes $\displaystyle F(\Omega)$
$\displaystyle=\sum_{j}F_{j}\exp(i\Omega\tau(\bm{x}_{j})-i\pi n_{j}).$ (37)

where $n_{j}=0,1/2,1$ keeps track of the extra factor of $-i$ arising from the
sign of the eigenvalues. This phase shift arises from the image parity and is
called the Morse phase; in the eikonal limit it is discretely-
valued666However, this is not necessarily true in the diffractive regime when
the saddle-point approximation breaks down – see Appendix A of (Jow et al.,
2022). It has been known for decades in (Blandford & Narayan, 1986; Nakamura &
Deguchi, 1999) but only recently have the observational ramifications become
observationally relevant, in particular in the parameter space of
gravitational-wave lensing (Dai & Venumadhav, 2017; Dai et al., 2020). The
overall Morse phase is not predicted by waveform models. In the specific
setting of gravitational waves generated by the nearly-circular inspirals of
compact binaries, an overall Morse phase shift of $180^{\circ}$ is degenerate
with a $180^{\circ}$ orbital phase shift of the binary as long as we restrict
our analysis to the dominant quadrupole mode of the radiation. Orbital
eccentricity or the inclusion of higher-order modes in model waveforms can
also break this degeneracy (Janquart et al., 2021; Ezquiaga et al., 2021).

The relative Morse phase between two images of the same waveform can be used
to reduce backgrounds in searches for lensed GWs (Dai et al., 2020). In the
time domain, a relative Morse phase of $180^{\circ}$ induces a sign flip of
the electromagnetic/gravitational waveform (Type II images; $n_{j}=1$)
relative to a Type I image. In the case of saddle points (Type III images;
$n_{j}=1/2$), the Morse phase applies a Hilbert transform to the unlensed
waveform. This latter distortion is non-trivial and leads to biases in
parameter estimation unless this is properly taken into account in search
templates (Vijaykumar et al., 2022).

Interference effects are imprinted in not only the phase, but also the
amplitude information from a lensing event. Interference induces ripple-like
modulations of the intrinsic spectrum of a source. In this scenario, knowing
the intrinsic spectrum allows us to measure $|F|^{2}$ as a function of
$\Omega$. Expanding Eq. 37, we compute

$|F|^{2}=\sum_{j}|F_{j}|^{2}+2\sum_{j<k}F_{j}F_{k}\cos(\Omega(\tau(\bm{x}_{j})-\tau(\bm{x}_{k}))-\pi(n_{j}-n_{k}))$
(38)

where we recognize $|F_{j}|^{2}=1/\det{A^{j}_{ab}}$ as the flux magnification
of each image present in geometric optics. Eikonal optics induces interference
terms which induce a ripple-like modulation, similar to pulsar scintillation,
in the spectrum of the observed source over a bandwidth
$\Delta\Omega=2\pi/(\tau(\bm{x}_{j})-\tau(\bm{x}_{k}))$, as shown in Fig. 3.
These modulations are potentially quite detectable due to the favorable
scaling of Eikonal optics over geometric optics (Jow et al., 2020; Takahashi &
Nakamura, 2003). To see this, consider a detector which detects the brightest
image, which has magnification $|F_{j}|^{2}$, and whose noise floor is
expressed as some fraction $\varepsilon^{2}$ of this:
$\varepsilon^{2}|F_{j}|^{2}$. In geometric optics, the detectability of a
dimmer image $k$ if it is distinguishable from the brighter one (e.g. has a
large angular separation on the sky) depends whether the ratio
$|F_{k}|^{2}/|F_{j}|^{2}\sim\varepsilon^{2}$. In Eikonal optics, however, an
image $k$ can be detected and spectrally separated from image $j$ using just
the combined spectrum of both images. The detection is made by searching for
the interference term in the combined spectrum, which induces sinusoidal
modulations in the spectrum whose size is
$|F_{j}F_{k}|/|F_{j}|^{2}\sim\varepsilon$.

The Eikonal limit, or stationary-phase approximation, is a particularly useful
approximation because the behaviour of the stationary-phase points, a.k.a. the
images, is described by the rich mathematics of catastrophe theory.
Catastrophe theory characterizes the topology of caustics, which are the
boundaries in parameter space at which the number of stationary-phase points
changes. The essence of what makes catastrophe theory so powerful is the
result that the topology of caustics can be classified by a small number of
“elementary catastrophes" (Thom, 1967). The result is that lenses will behave
qualitatively similarly based on which elementary catastrophes characterize
the lenses, even though they may differ significantly in functional form (Nye,
1999). While, technically, the Eikonal limit breaks down at caustics (at
caustics, two or more of the inverse magnification matrices are degenerate,
leading to formally infinite magnifications), the area of the region around
the caustic for which the Eikonal approximation fails goes to zero as
$\Omega\to\infty$. Moreover, the full wave optics diffraction patterns of some
of the elementary catastrophes are known (Berry & Upstill, 1980). For example,
as we will see in Sec. 9, the diffraction pattern near a fold catastrophe is
given exactly by the Airy function (see also (Nakamura & Deguchi, 1999)).
Therefore, the known diffraction patterns near caustics coupled with the
Eikonal limit away from the caustics can provide a full description of a
generic lens. This approach is taken in Grillo & Cordes (2018).

The effects mentioned so far – Morse phase shifts, spectral modulation of the
amplification factor, and favorable scaling of lensing cross-sections – arise
as generic interference phenomena present in Eikonal optics. Going into the
fully diffractive regime adds additional features which we will discuss next.

##  5 Diffractive optics ($\Omega\gtrsim 1$)

In the case where $\Omega$ is not much greater than $1$, Eq. 28 cannot be
simplified using the saddle-point approximation. For an axisymmetric lens, the
potential $\psi(\bm{x})\to\psi(x)$, and the amplification can be written as

$\displaystyle F(\Omega,y)=$ $\displaystyle-i\Omega\exp(i\Omega
y^{2}/2)~{}\times$ (39) $\displaystyle\int_{0}^{\infty}~{}dxxJ_{0}(\Omega
xy)\exp\left[i\Omega(x^{2}/2-\psi(x))\right]$ (40)

where $J_{0}$ is the zeroth-order Bessel function. In certain cases this may
be computed analytically. For a point mass where $\psi(x)=\log(x)$ this
simplifies to

$\displaystyle F(\Omega,y)$
$\displaystyle=\exp\left[\dfrac{\pi\Omega}{4}+i\dfrac{\Omega}{2}\left(\ln\left(\dfrac{\Omega}{2}\right)\right)\right]\times$
(41)
$\displaystyle\Gamma\left(1-\dfrac{i\Omega}{2}\right)~{}_{1}F_{1}\left(\dfrac{i\Omega}{2},1;\dfrac{i\Omega
y^{2}}{2}\right).$ (42)

In this equation (see also Eq. 21 of (Deguchi & Watson, 1986a), Eq. 2.16 (Katz
et al., 2018) for slightly different forms), $\Gamma$ is the gamma function
and ${}_{1}F_{1}$ is the hypergeometric function777Wolfram MathWorld:
Hypergeometric ${}_{1}F_{1}$.

The point mass is the simplest analytical solution for the amplification
factor in wave optics in gravitational lensing. However, other analytical
solutions exist which we briefly mention: the solution for a Kerr (rotating)
black hole is very similar to that of the point mass; it includes a shift in
the point-mass lensing potential
$\psi(\bm{x})=\log(x)-\dfrac{\bm{x}\cdot\bm{\alpha}}{x^{2}}$, where the shift
$\vec{\alpha}=\hat{n}\times\vec{J}/(MR_{E})$ depends on the physical Einstein
radius $R_{E}$ and the lens’s angular momentum vector per unit mass
$\vec{J}/M$ relative to the line-of-sight unit vector $\hat{n}$. The
interference pattern is the same as that of Eq. 42, being circularly-
symmetric, albeit shifted by an angle $\alpha$ (Baraldo et al., 1999). This is
similar to the Sagnac effect for a rotating optical interferometer. Analytic
solutions for the singular isothermal sphere (Matsunaga & Yamamoto, 2006), an
infinite cosmic string (Suyama et al., 2006), and a one-dimensional sinusoidal
phase screen (Beach & Lovelace, 1997) are also known.

Solutions to the diffraction integral beyond the analytical ones presented
here are often computed numerically. However, this can be difficult because
the phase in integrals such as Eq. 28 oscillates faster for larger values of
the argument $\bm{x}$. Some of the first numerical efforts based on contour
integration were developed for the calculation of femtolensing spectra (Ulmer
& Goodman, 1995); this method has also been applied successfully to
gravitational waves (Mishra et al., 2021). Approximate methods using Eikonal
optics (Grillo & Cordes, 2018) are effective for large lensing phase shifts
($\Omega T\gg 1$) and near caustics (see also discussion of Sec. 1 of (Jow et
al., 2021)). Grillo & Cordes (2018) also point out that the Fresnel integral
can be written as a convolution and evaluated via FFT. This performs well near
caustics, but requires the lens plane to be densely sampled with a finite and
regular two-dimensional grid fine enough to capture the details of the
oscillatory integrand.

More recently, Picard-Lefschetz theory has been used to assist in the analysis
and evaluation of Fresnel integrals (Feldbrugge et al., 2019; Feldbrugge &
Turok, 2020; Jow et al., 2021; Shi & Xu, 2021; Jow et al., 2022). These
methods essentially involve deforming the two lens plane integrals into a set
of carefully-chosen contours in the complex plane along which contour
integrals are non-oscillatory and absolutely convergent. One feature of
Picard-Lefschetz integration is that while it converges slowly for $\Omega\sim
1$, performance improves as the integrands become more oscillatory ($\Omega\gg
1$), so it is complementary to many other techniques. A limitation is that the
lens phase function $\phi(\bm{x})$ (e.g.
$\phi(\bm{x})=\Omega\tau(\bm{x},\bm{y})$) must be able to be evaluated when
the domain is complexified from $\mathbb{R}^{2}\to\mathbb{C}^{2}$. This is
straightforward for the broad class of rational functions (quotients of
polynomials), which have no branch cuts (also a desirable property) and a
small number of isolated poles. The point mass lens has $\psi(x)=\log(|x|)$
which does have a branch cut; in a follow-up work to (Feldbrugge et al.,
2019), a technique to evaluate the isolated and the binary point lens
potentials by reparametrizing the radial integral was reported (Feldbrugge &
Turok, 2020). However, well-known lens models like the singular isothermal
sphere ($\psi(x)=|x|$) have not yet been evaluated using this method. Despite
this limitation, the utility of Picard-Lefschetz theory is not limited to its
numerical application, but also includes its conceptual power. By extending
the analysis of optics into the complex plane, Picard-Lefschetz theory gives a
well-defined prescription for separating contributions to the total observed
flux from individual geometric images, even at low frequencies. Picard-
Lefschetz theory can be seen as a way to extend the “discrete images”
description from geometric optics lensing deep into the wave regime. In this
way, Picard-Lefschetz theory can conceptually bridge the gap between these two
regimes(Jow et al., 2022).

In Fig. 4 we show several magnification maps, evaluated using the Picard-
Lefschetz method for several values of $\Omega$. The map coordinates
$\bm{\mu}=(\mu_{x},\mu_{y})$ are related to the source plane and observer
plane coordinates ($\bm{y}$ and $\bm{y}_{\textrm{obs}}$ respectively) by
$\bm{\mu}=\bm{y}\dfrac{D_{L}}{D_{S}}+\bm{y}_{\textrm{obs}}\dfrac{D_{LS}}{D_{S}}$.
This parametrization reflects the symmetry of the source and the observer in
optics: in gravitational lensing we may think of fixing
$\bm{y}_{\mathrm{obs}}=0$, and treat $\bm{\mu}$ as a source plane coordinate,
while in a laboratory diffraction experiment we may instead fix $\bm{y}=0$ and
observe the diffraction pattern on the observer’s screen, which spans many
values of $\bm{y}_{\mathrm{obs}}$. The Fresnel integral evaluated is that of
the Chang-Refsdal lens (Chang & Refsdal, 1979), for which

$\tau(\bm{x},\bm{y})=\dfrac{1}{2}|\bm{x}-\bm{\mu}|^{2}-\log(|\bm{x}|)-\dfrac{\gamma}{2}(x_{1}^{2}-x_{2}^{2})$
(43)

whose single free parameter $\gamma>0$ describes the external shear.

Figure 4: A map of $|F|^{2}(\Omega,\bm{\mu})$ for a point mass+external shear
model (Eq. 43), where the map coordinates $\bm{\mu}$ are the scaled source
plane coordinates $\bm{y}$. Columns: From left to right, each column
corresponds to $\Omega=100,75,50,25$. Rows: From top to bottom, each row
corresponds to $\gamma=0.1,0.3,0.5$. The caustic, at which the geometric
optics magnification diverges, is overlaid as a white diamond on each map.
Reproduced with permission from source code used for (Feldbrugge & Turok,
2020).

###  5.1 A useful analytical approximation

The point mass model is extremely useful, but in practice, hypergeometric
functions with complex arguments can be cumbersome to numerically evaluate
(Pearson et al., 2014). There is also a useful “semiclassical” approximation
to Eq. 42 for situations when $\Omega\gg 1$ but when the sinusoidal form of
Eq. 37 is too simple. Solving for the stationary points of

$\displaystyle\tau(x)=$ $\displaystyle\dfrac{x^{2}}{2}-\psi(x)$ (44) gives the
geometric-optics image positions $x_{E}$ such that
$x_{E}=\psi^{\prime}(x_{E})$. To second order in $\Delta x=x-x_{E}$,
$\displaystyle\tau(x)$
$\displaystyle=\tau(x_{E})+\dfrac{1}{2}(1-\psi^{\prime\prime}(x_{E}))\Delta
x^{2}=\tau(x_{E})+(1-\kappa)\Delta x^{2}.$ (45) The region of the lens plane
where the rays interfere constructively is defined by $\displaystyle
2\pi\gtrsim$ $\displaystyle(1-\psi^{\prime\prime}(x_{E}))\Delta x^{2}\Omega,$
(46) and this domain has a width of $\Delta
x\sim(\Omega|1-\psi^{\prime\prime}(x_{E})|)^{-1/2}$. Within this domain, the
Bessel function argument of Eq. 40 changes by $\Omega\Delta xy$, and outside
this domain, it oscillates quickly. If $y\lesssim\Omega^{-1/2}$, then
$1\gtrsim\Omega\Delta xy$, and the Bessel function does not change over the
important range of the integration. Pulling $xJ_{0}(\Omega xy)$ out of the
integral sign leaves a Gaussian integral from $0$ to $\infty$:
$\displaystyle|F_{\textrm{wave}}|^{2}$ $\displaystyle\approx\pi\Omega
x_{E}^{2}|1-\psi^{\prime\prime}(x_{E})|^{-1}J_{0}^{2}(\Omega xy)$ (47)

$\leq\pi\Omega x_{E}^{2}|1-\psi^{\prime\prime}(x_{E})|^{-1}=2\pi\Omega
x_{E}^{2}|1-\kappa(x_{E})|^{-1}.$

where we have used $2(1-\kappa(x_{E}))=1-\psi^{\prime\prime}(x_{E})$.

##  6 Angular broadening/finite-size effects

So far, we have assumed a point-like source of monochromatic radiation in the
source plane. In reality, we observe radiation from a finite patch of
$\bm{y}$-values in the source plane. In general, integrating over a finite
patch of the source plane washes out the chromatic interference patterns seen
in both the eikonal and wave optics regimes, e.g. Eqs. 38 and 42. This is the
main obstacle to observation of interference effects in gravitational lensing.

Before we dive into a deeper discussion of decoherence, it is worth defining
what type of coherence we care about. Under one definition, “phase-coherence”
is strictly a property of the emitted radiation, which can be considered as a
time series. If one can predict the field amplitude as a function of time (or
frequency, if one prefers the Fourier domain) with a model, that source would
be considered phase-coherent.

Gravitational waves are an excellent example of a phase-coherent source by
this definition. Given the time of the merger, it is possible to write down a
model for the strain amplitude as a function of time. The reason this is
possible is that the emitter is a simple binary system, and that the emission
region is much smaller than a GW wavelength. It is therefore possible to
coordinate the entire emission region at once.

Most electromagnetic waves emitted by astrophysical sources are not phase-
coherent. Quantitatively, this is because the transverse size of virtually any
astrophysical emission region is usually much larger than a wavelength of
electromagnetic radiation. This means that different patches of the emission
region cannot be causally connected, i.e. synchronized to better than one
wavelength. While individual emitters may be phase-coherent (e.g. a single
electron orbiting a field line emitting synchrotron radiation is phase
coherent), the fact that there exists an ensemble of emitters, all radiating
with some random phase with respect to each other, means that the phase of the
total time series is not predictable.

Source emission regions do not need their transverse sizes to smaller than one
wavelength (i.e. they do not need to be phase coherent) in order for wave-
mechanical features to be observed. The actual requirement on the source size
is much looser; it is a constraint on the angular coherence of the source and
involves the observer’s geometry in addition to the emission physics. In the
terminology of interferometry and interstellar optics, this is similar to a
source appearing “unresolved” when viewed with a particular receiver setup or
through an inhomogeneous medium with a certain geometry.

Whether a source is angularly-coherent, or unresolved, depends the apparent
angular extent of the source. In some cases this corresponds to the true size
of the emission region. However, in the case of radio waves, interstellar
scattering can enlarge the apparent size of the source via angular broadening.
In addition, in areas of the source plane with high magnifications (caustics),
small displacements in the source plane correspond to large changes in the
image positions on the lens plane (and therefore the lensing time delays).
These effects all work to increase the apparent angular size of the source.
Astrophysical masers, however, may have physical sizes which are larger than
their apparent sizes (Goldreich & Keeley, 1972), which have been observed to
be inversely correlated with their isotropic luminosities (Johnston et al.,
1997).

We begin with a generic analysis of angular decoherence in the eikonal and the
wave optics regimes, and study specific cases later (point caustics in Sec. 8
and briefly comment on the case of fold caustics in Sec. 9). Our quantitative
criteria for angular decoherence is whether the chromatic oscillations (e.g.
those depicted in Fig. 3) are washed out by integrating over a finite patch of
the source plane. We require the total phase across the extent of the source
to change by $\leq 1$ radian, when viewed at a distance through the
gravitational lens. We find nontrivial differences by considering this
criteria in the Eikonal and the diffractive regimes.

Figure 5: In the eikonal limit, a lensed source remains coherent if the
spectral oscillations do not significantly change over the source’s finite
angular size, here shown as a blue patch. When the images are well-separated,
“rays” from the source traverse the lens plane through the image positions
$\bm{x}_{\pm}$ calculated via geometric optics. For the point-mass lens, one
image passes outside the Einstein ring at a transverse physical scale $R_{E}$,
and the other passes inside the Einstein ring on the opposite side of the
lens. The interference fringes from lensing are present if the source looks
point-like as observed by an interferometer with a baseline of Eq. 51 (i.e. a
pair of antennas placed in the lens plane at the apparent locations of the
gravitational lens images in the lens plane).

In the Eikonal limit that the stationary points are well-separated into
discrete images $j$ and $k$, there exists a well-defined time delay
$T_{j}-T_{k}=T_{jk}$ which can be calculated for a point-like source at any
location in the source plane. We require that when integrating over a patch of
radius $\Delta\bm{y}$, changes in $T_{jk}$ must be small. In other words,
different patches of the source must not experience different time delays;
otherwise, we are averaging over many fringe spacings in the spectrum as
depicted in Fig. 3. Denoting $\overline{\bm{y}}$ as the displacement from the
center of the source, and $T_{j}$ and $T_{k}$ denote the time delays of two
distinct images generated by the gravitational lens. The phase shift is

$\displaystyle 1\geq\Delta\varphi$ $\displaystyle=\Omega\dfrac{\partial
T_{jk}}{\partial y_{a}}\cdot\overline{\bm{y}}_{a}$ (48)
$\displaystyle=\Omega(\bm{x}_{j}-\bm{x}_{k})_{a}\overline{\bm{y}}_{a}.$ (49)
Using Eq. 26 gives the width of the Fresnel zone: $\displaystyle\Delta y$
$\displaystyle\lesssim\dfrac{1}{\theta_{E}}\dfrac{\lambda}{2\pi(1+z_{L})D_{\textrm{img}}}.$
(50) $D_{\mathrm{img}}$ resembles the Rayleigh resolution limit for a radio
interferometer observing at a wavelength of $\lambda/(1+z_{L})$ with baseline
corresponding to the physical scale of the image separation $\displaystyle
D_{\textrm{img}}$
$\displaystyle=\dfrac{D_{L}D_{S}}{D_{LS}}|\bm{\theta}_{j}-\bm{\theta}_{k}|,.$
(51)

This can be interpreted as the effective separation of two telescopes located
in the lens plane (Schneider, 1983), at the geometric-optics positions of the
images, observing the source at a frequency of $\Omega(1+z_{L})$, as shown in
Fig. 5. The exact form of the fall-off is studied in (Matsunaga & Yamamoto,
2006), and applied to the particular case of femtolensing in (Stanek et al.,
1993). In the diffractive limit, where the images in the time domain are
separated by less than $\approx 1/\Omega$, the two stationary points blend
together coherently, i.e. they belong to the same Fresnel zone (see e.g. the
image positions denoted in the top left panel of Fig. 1). Therefore, the
criterion we derived in Eq. 50 treating the paths as well-separated does not
apply. In this regime, in the absence of a template for the unlensed signal,
it does not make sense to think of the two images as well-separated,
distinguishable paths through the lens plane, because all parts of the lens
plane within the first few Fresnel zones contribute significantly to
$F(\Omega,\bm{y})$.

Figure 6: Regimes of coherence in gravitational lensing. We plot the apparent
radius of the source against the observing wavelength (ignoring redshift
effects). Solid lines: If the lens mass is in the eikonal limit, then the
source must be below the solid lines corresponding to Eq. 50. Dashed lines: In
the diffractive/wave optics regime, the source size must be below the dashed
lines corresponding to Eq. 56. We have roughly labeled the apparent sizes of
various sources and wavelengths at which they have been observed in the
various boxes (Oguri, 2019; Leung et al., 2022; Katz et al., 2018).

If slightly different patches of the source have very different image
configurations, this could lead to very different time delays. One might
expect this effect to be small, but there are scenarios where this might
become significant. Near the Einstein ring caustic in a point-mass lens, for
example, small displacements in the source plane can lead to large differences
in the image positions.

Near the $j$th stationary point of $T$, where $dT/d\bm{x}=0$, we analyze the
curvature of $T$ with respect to $\bm{x}$:

$\Delta\varphi=\Omega\Delta T=\frac{1}{2}\Omega\frac{\partial^{2}T}{\partial
x_{a}\partial
x_{b}}\overline{\bm{x}}_{a}\overline{\bm{x}}_{b}={\textstyle\frac{1}{2}}\Omega(A^{j}_{ab})^{-1}\overline{\bm{x}}_{a}\overline{\bm{x}}_{b}$
(52)

where $A^{j}_{ab}$ is the magnification matrix evaluated at the stationary
point $j$. Hence the allowable displacement on the image plane that causes a
one-radian phase change is defined by

$1\geq\Delta\varphi=\Omega
A^{-1}_{ab}\overline{\bm{x}}_{a}\overline{\bm{x}}_{b}$ (53)

Now we equate $\overline{\bm{x}}_{a}$ to the shift in the geometric-optics
image for the source whose extent is $\overline{\bm{y}}_{c}$. From the
definition of the magnification matrix, we have

$\overline{\bm{x}}_{a}=A^{j}_{ac}\overline{\bm{y}}_{c}$ (54)

Eq. 53 then becomes

$1\geq\Delta\varphi=\dfrac{1}{2}\Omega
A_{cd}\overline{\bm{y}}_{c}\overline{\bm{y}}_{d}.$ (55)

Substituting Eq. 26 gives our estimate of the Fresnel zone’s shortest axis, if
$\lambda_{\mathrm{max}}$ is the largest eigenvalue of $A_{cd}$:

$\Delta
y\lesssim\dfrac{1}{\theta_{E}|\lambda_{\mathrm{max}}|}\sqrt{\dfrac{\lambda}{2\pi(1+z_{L})D_{\mathrm{eff}}}}$
(56)

where $D_{\mathrm{eff}}=\dfrac{D_{L}D_{S}}{D_{LS}}$. Eq. 56 tells us that if
the source size is on the order of the Fresnel angular scale
($\sqrt{\lambda/D_{\mathrm{eff}}}$), the pattern in Fig. 6 will get washed out
because the image positions shift significantly over the source’s finite
extent, yielding different delays relative to the lensing delay evaluated
towards the center of the source. Magnification plays a role in that higher
flux magnifications shrinks the Fresnel scale in the source plane. We
illustrate Eqs. 50 and 56 in Fig. 6, using $R_{src}\sim\Delta
y\theta_{E}D_{\mathrm{eff}}$ and ignoring factors of $1+z_{L}$. We see that
Eq. 50 is the more stringent condition that applies in most astrophysical
situations, though for short wavelengths and extreme magnification scenarios
(e.g. proposals to use the sun as a gravitational lens (Maccone, 2010;
Turyshev, 2017; Engeli & Saha, 2022)), the latter may become relevant.

##  7 Exploiting Chromaticity in Lensing

As illustrated in Fig. 3, one powerful feature of interference effects in
lensing is the introduction of an additional length scale (the wavelength)
into the problem. In the eikonal and wave optics regimes, the magnification
factor picks up an oscillatory dependence on the wavelength (through $\Omega$)
and geometry (through $\bm{y}$). By observing the source over a wide enough
range of wavelengths to measure the spectral oscillations, more astrophysical
information can be recovered than in the geometric-optics case alone. One
particularly interesting setup is a multi-messenger observation of a lensed
gravitational-wave + electromagnetic source. This would allow for exquisite
measurements of wave optical effects, since due to their low $\Omega$,
gravitational waves are potent probes of diffractive effects, while
electromagnetic waves are typically in the eikonal or the geometric optics
regimes. Observing simultaneously-emitted low-frequency (gravitational) images
and high-frequency (electromagnetic) images can expose wave optics effects
such as modifications to the phase velocity of diffractively-lensed images
(Takahashi, 2017), though analyses of the causal “wavefront” (Suyama, 2020;
Ezquiaga et al., 2020) reveals that no violation of causality occurs, though
an apparent violation of causality may arise from incorrect separation of the
observed image superposition into individual images (Tanaka & Suyama, 2023).

Another example is using the chromaticity of the interference pattern to break
the mass-sheet degeneracy in the eikonal and the wave optics regimes.
Cremonese et al. (2021) shows that the point mass can be disentangled from a
mass sheet in both regimes using the spectral dependence of the amplification
factor (the oscillatory terms in Eq. 57 and Eq. 47, which encode information
about the mass scale $M$). In the diffractive limit, even though the spectral
dependence becomes more complex than in the simple eikonal limit, the lensed
images become more difficult to separate from one another since their temporal
separation is less than one wavelength.

In a similar vein, the chromatic oscillations can be used to directly measure
time delays in lensing in both the eikonal and diffractive limits. This idea
is applicable to various physical systems, e.g. in radio pulsar/fast radio
bursts lensed by sub-solar mass compact objects and planets (Jow et al., 2020)
and sub-Hertz gravitational waves generated by chirping binaries lensed by
$\sim 10^{8}M_{\odot}$ objects (Takahashi & Nakamura, 2003). For example,
consider a microlensing scenario where a transiting foreground lens with
proper motion $\mu$ achieves a minimum impact parameter $y_{0}$ at a time
$t_{0}$. In diffractive optics, the magnification varies as a function of
$(t,\Omega,t_{E},y_{0},t_{0})$, where
$t_{E}\approx\dfrac{1}{\mu}\sqrt{R_{E}/D_{\mathrm{eff}}}\propto\sqrt{M/D_{\mathrm{eff}}}/\mu$
is the Einstein crossing time of the lens. The lens mass enters through both
$\Omega$ and $t_{E}$. In geometric optics, the magnification varies as a
function of ($t,t_{E},y_{0},t_{0}$) only (since there is no frequency
dependence). Since the mass dependence enters only through $t_{E}$ it is
perfectly degenerate with the lens distance $D$. Without an independent
distance measurement, it is difficult to measure the lensing mass in this
situation.

These examples show that there is clearly additional information encoded by
the frequency dependence, it is helpful to conceptualize, using a thought
experiment, where the information really comes from. Consider a delta-function
electric-field pulse emitted by a source which gets lensed and subsequently
recorded by a detector which records the raw electric field as a function of
time. In this situation it would be possible to simultaneously measure the
magnification ratio from the light curves, as well as the Shapiro delay ($\sim
20\mu\mathrm{s}\times(M/M_{\odot})$), which would give a handle on the mass.
By Nyquist’s theorem, the time resolution of the Shapiro delay measurement is
limited by the inverse bandwidth ($1/\Delta\nu$) of the detector. Larger
bandwidths allow measuring $F(\Omega,y)$ for a wide range of $\Omega$ at a
single value of $y$; this multi-frequency coverage allows the mass to be
constrained. As was shown earlier in Eq. 37, the oscillations in $F(\Omega,y)$
vary on scales proportional to the inverse Shapiro delay. In a Fourier-
transform sense, measuring light curves over a wide range of wavelengths
provides a new way to directly constrain the Shapiro time delay, and therefore
the mass scale. This is completely distinct from the method of measuring the
characteristic timescale of a lensing transit, which indirectly probes the
mass via long-term monitoring and measurement of the Einstein crossing
timescale and the relative velocities of the source and observer.

These Shapiro delay measurements are not limited to hypothetical delta-
function pulses, or even astrophysical transients. Shapiro delay measurements
can be made on non-variable sources, if they are compact enough, since the
measurement uses fluctuations in the electric field rather than the total
power. However, the stringent size constraints on preserving the interference
fringes excludes most, if not all, currently known persistent sources from
this type of measurement. Transient sources from compact objects are precisely
those which are known to be sufficiently point-like to realize the
interferometric measurement of Shapiro delays using the chromatic modulations
of the dynamic spectrum.

In summary, using chromaticity to extract astrophysical information from lens
systems offers an interesting trade-off – the stringent requirements for
source compactness demand transients, which often do not repeat. For these
one-off transients, we have given up our ability to monitor the Einstein
crossing timescale for a one-time, interferometric measurement of the Shapiro
delay timescale free of uncertainties from estimating the relative velocities
of the source and lens. Using a lens model, the Shapiro delay and the flux
ratio completely determine the (redshifted) lens mass and impact parameter
(see also (Itoh et al., 2009)). In the case of persistent emitters (e.g.
optical microlensing) or recurring transient events such as pulsar pulses,
measurement of a third quantity (the Einstein crossing timescale) can be used
to extract additional information (such as the relative motion of the lens and
the source).

Figure 7: Top left: Frequency-dependent amplification
$|F_{\textrm{wave}}(\Omega,y)|^{2}$ spectrum of a compact, broadband radio
sources as it transits behind a point-mass lens of $5M_{\otimes}$, $D=1$ kpc
with a minimum impact parameter of $y_{0}=0.5$. The amplification oscillates
like $\cos^{2}(\Omega(\tau_{+}-\tau_{-}))$ as a function of frequency for
large impact parameters where the Eikonal optics (see Eq. 38) approximation
does not break down. Top right: Fractional discrepancy between
$|F_{\textrm{wave}}(\Omega,y)|^{2}$ and $|F_{\textrm{geo}}(\Omega,y)|^{2}$ for
the point-mass lens. Diffractive effects distort the regular fringe pattern
predicted from geometric optics. Bottom left: The Fourier transform of the top
left panel along the frequency axis, also referred to as the “secondary
spectrum”. The arc tracks the instantaneous lensing delay
$(\tau_{+}(t)-\tau_{-}(t))$ as it changes throughout the transit. Bottom
right: The Fourier transform of the top right panel along the frequency axis.
Adapted with permission from (Jow et al., 2020).

Another way of understanding this tradeoff is by analogy to radio
scintillation (see e.g. (Narayan, 1992; Goodman et al., 1987)). The
interference pattern in Fig. 3 is similar to diffractive scintillation in
pulsar astronomy, where inhomogeneous plasma in front of a compact background
source causes an intrinsically featureless spectrum to exhibit stochastic
ripples over some characteristic bandwidth. Like in radio scintillation, the
only sources which “twinkle”, or scintillate, are not necessarily those which
show intrinsic variability, but those which appear compact and unresolved by
the cosmic telescope created by the plasma inhomogeneities: a star need not
exhibit intrinsic variations to twinkle in our atmosphere.

To push the analogy further, another manifestation of scintillation is
refractive scintillation, where the clumpy spectrum caused by diffractive
scintillation changes slowly over longer timescales as the inhomogeneities
drift with respect to the observer with some characteristic crossing time
(hours to days for scintillation of Galactic pulsars). Monitoring the dynamic
spectrum of a source over long timescales probes the motion of plasma screens
via scintillation (see e.g. (Walker et al., 2004)), similar to how one might
measure Einstein crossing timescales in lensing. Scintillation shows up as
parabolic, arc-like features in the time-lag correlation function, estimated
from the Fourier transform of the source’s spectrum, when the time-lag
correlation function is concatenated over many days such that it has both a
delay axis (measuring timescales on the scale of the inverse bandwidth) and a
time axis (measuring timescales of days or weeks).

A drifting gravitational lens, or a collection of many gravitational lenses
drifting together, can similarly produce a parabolic arc in this space (see
lower panel of Fig. 7) as it crosses the lens. The effect of a collection of
moving lenses causing gravitational scintillation has been considered;
(Macquart, 2004) derives analytical expressions for the phase perturbations
sourced by clumpy matter distributions; (Congedo et al., 2006) considers the
particular case of gravitational scintillation by stellar clusters in the
Galactic center. (Takahashi et al., 2005) considers the validity of the Born
approximation and the impact of chromatic (diffractive) effects in
gravitational scintillation. Diffractive effects as a probe of a clumpy matter
distributions on small scales are considered in (Takahashi, 2006; Oguri &
Takahashi, 2020; Inamori & Suyama, 2021; Oguri & Takahashi, 2022).

##  8 Finite magnifications near the Einstein ring

Having generically analyzed eikonal and wave optics and the transition to
incoherent (traditional) gravitational lensing, we can now turn to studying
more specific phenomena which arise in wave optics. One failure of geometric
optics is that for certain regions of the source plane, magnification ratios
become infinite. These regions are typically isolated points or one-
dimensional curves in the source plane where geometric-optics rays “cross” and
as such are regions where interference effects are most crucial. These are
caustics (Ohanian, 1983): the most famous example of a caustic is the Einstein
ring for a point-mass lens. Rays emanating from the point caustic at
$\bm{y}_{0}=0$ pierce the lens plane at a ring of points called the “critical
curve” (here, the critical curve is defined by $|\bm{x}_{0}|=1$) on their way
to the observer.

The infinite magnification of geometric optics lensing is removed by two
distinct effects: finite source size (only an infinitely small patch is
magnified infinitely), and diffraction (Herlt & Stephani, 1976; Benson &
Cooke, 1979). In other words, to truly attain infinite magnification, one
would need an infinitely-small point source which radiates at an infinitely-
short wavelength. In Fig. 8 we calculate magnification maps for both a point
source close to the caustic, and a finite-sized source centered on the
caustic, in the case of infinite and finite frequency (geometric optics and
wave optics respectively).

The magnification profiles in geometric optics and diffractive optics can be
computed as follows. In geometric optics the total magnification of a point
source at some impact parameter $y$ is

$\displaystyle|F_{\textrm{geo}}|^{2}$
$\displaystyle=\dfrac{y^{2}+2}{y\sqrt{y^{2}+4}}~{}$ (57) To obtain the finite-
size result, we average Eq. 57 over a finite-sized top-hat of radius $\Delta
y>0$ in the source plane centered at $y=0$. This gives a finite total
brightness. The amplification factor for a circular, on-axis source of radius
$\Delta y$ is $\displaystyle\langle|F_{\textrm{geo}}|^{2}\rangle_{\Delta y}$
$\displaystyle=\sqrt{1+\dfrac{4}{\Delta y^{2}}}~{}$ (58) In wave optics, we
repeat the top-hat averaging over the point-mass magnification map derived
from wave optics (Eq. 47). This yields
$\displaystyle\langle|F_{\textrm{wave}}(\Omega)|^{2}\rangle_{\Delta y}$
$\displaystyle=\dfrac{\pi\Omega x_{E}^{2}[J_{0}(\Omega x_{E}\Delta
y)^{2}+J_{1}(\Omega x_{E}\Delta y)^{2}]}{|1-\psi^{\prime\prime}(x_{E})|}~{}$
(59)

In Fig. 8, we plot Eqs. 57 and 47 in the top panel, and Eqs. 59 and 58 in the
bottom panel. While both diffraction and finite source size effects modify the
behavior of $|F|^{2}$ near $y=0$, they have very different qualitatitive
features – even a tiny departure from a point source washes out the large
oscillations characteristic of diffractive optics. It is also evident that
when the source size is on the caustic, the finite source size solution
deviates considerably from the point source solution (Witt & Mao, 1994;
Sugiyama et al., 2020). The discrepancy between using a point-source formula
and the extended source matter most when the source dimension $\Delta
y\lesssim\Omega^{-1}$; this is most relevant for gravitational waves, which
have low values of $\Omega$ and have source sizes
$\lesssim$100\text{\,}\mathrm{km}$$ (Oguri, 2019).

Figure 8: The behavior of the magnification $|F|^{2}$ taking into account
finite source size and wave optics effects near the Einstein ring ($y\to 0$).
Top: The magnification map as a function of $\bm{y}$ near the origin saturates
at $\bm{y}\sim 1/\Omega$ at a magnification of $\pi\Omega$. Bottom: the
magnification map for an on-axis source, now as a function of source radius
$\Delta y$.

##  9 Finite magnifications near fold caustics

The well-known Einstein ring is an example of a caustic: loci where the
geometric optics magnification diverges. Caustics may be isolated points or
one-dimensional curves in the source plane. Similar to the Einstein ring, the
infinite magnifications near “fold caustics” get regularized by finite source
size as well as diffraction. Here, following (Nakamura & Deguchi, 1999), we
consider the Fresnel integral (Eq. 28) in the scenario where a source at
$\bm{y}_{0}$ lies directly on the caustic. We expand the time-delay potential
about a point $\bm{x}_{0}$ on the critical curve. We can choose coordinate
systems in both the source and lens planes with unit vectors
($\hat{\bm{y}}_{1}$, $\hat{\bm{y}}_{2}$), ($\hat{\bm{x}}_{1}$,
$\hat{\bm{x}}_{2}$). What does the time delay potential look like near the
critical curve in the neighborhood of $\bm{x}_{0}$? We can begin with a third-
order Taylor expansion about $\bm{x}_{0}$ and reason through each term. Here,
$\overline{\bm{x}}=\bm{x}-\bm{x}_{0}$, and
$\overline{\bm{y}}=\bm{y}-\bm{y}_{0}$, and the indices $a,b,c$ run over the
two dimensions of the lens plane.

$\displaystyle\tau(\bm{x},\bm{y})-\tau(\bm{x_{0}},\bm{y})$
$\displaystyle=\dfrac{\partial\tau}{\partial
x_{a}}\bm{\overline{x}_{a}}+\dfrac{1}{2}\dfrac{\partial^{2}\tau}{\partial
x_{a}\partial x_{b}}\bm{\overline{x}_{a}}\bm{\overline{x}_{b}}$ (60)
$\displaystyle+\dfrac{1}{6}\dfrac{\partial^{3}\tau}{\partial x_{a}\partial
x_{b}\partial
x_{c}}\bm{\overline{x}_{a}}\bm{\overline{x}_{b}}\bm{\overline{x}_{c}}$ (61)
The gradient of the potential is straightforward to evaluate at $\bm{x}_{0}$
and is $-\overline{\bm{y}}$. The magnification at $\bm{x}_{0}$
($\det(A^{0}_{ab})$) diverges if $\bm{y}=\bm{y}_{0}$; one of the eigenvalues
of $A^{0}$ (say, $\lambda_{1}$) is zero. The potential is therefore free of
curvature along the $\hat{\bm{x}}_{1}$ direction (thought it has a slope of
$-\overline{y}_{1}$). Along the $\hat{\bm{x}}_{2}$ direction, the curvature is
$1-\psi_{22}$, where the subscripts in $\psi_{22}$ refer to the two
derivatives taken (e.g. $\psi_{22}=\partial^{2}_{x_{2}}\psi$).
$\displaystyle\tau(\bm{x},\bm{y})-\tau(\bm{x_{0}},\bm{y})$
$\displaystyle\approx-\overline{\bm{x}}_{0}\cdot\overline{\bm{y}}_{0}+\dfrac{1-\psi_{22}(\bm{x_{0}})}{2}x_{2}^{2}$
(62) As we turn $\bm{y}\to\bm{y}_{0}$, the first derivative term vanishes.
This forms an image at $\bm{x}_{0}$ whose magnification is roughly the size of
the Fresnel zone (the patch in the lens plane over which the phase changes by
$2\pi$; see Fig. 1 for visualization). Traversing the $\bm{\hat{x}}_{2}$
direction (fixing $x_{1}=0$), the second derivative of the potential sets the
width of the Fresnel zone as usual. Along the $\bm{\hat{x}}_{1}$ axis (fixing
$x_{2}=0$), since the potential has no second derivative, the third derivative
sets the extent of the Fresnel zone along the $\bm{\hat{x}}x_{1}$ direction:
$\displaystyle\partial^{2}\tau/\partial x_{1}^{3}$
$\displaystyle=-\psi_{111}(\bm{x}_{0}).$ (63) If we define
$p=(1-\psi_{22}(\bm{x}_{0}))^{-1}$ and $q=-2/\psi_{111}(\bm{x}_{0})>0$, the
integrand is a quadratic along the $x_{2}$ direction and cubic along the
$x_{1}$ direction: $\displaystyle F=\dfrac{\Omega}{2\pi i}$ $\displaystyle\int
d\overline{x}_{1}\exp\left[i\Omega\left(-\overline{x}_{1}\overline{y}_{1}+\dfrac{\overline{x}_{1}^{3}}{3q}\right)\right]\times$
(64) $\displaystyle\int
d\overline{x}_{2}\exp\left[i\Omega\left(-\overline{x}_{2}\overline{y}_{2}+\dfrac{\overline{x}^{2}_{2}}{2p}\right)\right]$
(65) Taking $\mu_{*}=|p|q^{2/3}\Omega^{1/3}$, $Y=q^{-1/3}\Omega^{-2/3}$, and
letting Ai denote the Airy function999see Wolfram MathWorld, the amplification
factor is $\displaystyle|F|^{2}=2\pi\mu_{*}\textrm{Ai}(-\overline{y}_{1}/Y)$
(66)

The characteristic fringe spacing near the fold caustic is $\sim
Y\propto\Omega^{-2/3}$, which has a weaker frequency dependence than the point
caustic $(\sim 1/\Omega)$. Jaroszynski & Paczynski (1995) discusses
diffractive effects in caustic crossings at optical wavelengths, and points
out that it is possible to use the caustic crossing timescale to measure the
finite linear size of a source.

##  10 Modified lensing probabilities in wave optics

Modifying the amplification factor’s dependence on the impact parameter near
the Einstein ring modifies any optical depth calculations for which the
observations are flux-limited. For forecasts and constraints, it is important
to characterize the conditions under which diffraction modifies the optical
depth relative to the standard case of geometric optics. One effect, discussed
previously in Sec. 4, is that interference terms, if they are not washed out,
are easier to detect than autocorrelation terms. With a search method that
utilizes spectral information, this leads to an increase in the optical depth
(Katz et al., 2020; Jow et al., 2020; Kader et al., 2022) which becomes
significant when the source is much brighter than the limiting flux threshold.

In Sec. 8, we discussed how the infinite magnifications in the point-mass
model which appear in geometric optics get regularized by wave optics. One
might initially suspect that this decreases the optical depth in a wave optics
calculation relative to the geometric-optics calculation, but the reality is
that it depends on the magnification scale (Sugiyama et al., 2020), and the
angular motion of the source over the observation period. Diffraction
eliminates the possibility of chance alignments creating extremely high
magnification events ($|F|^{2}\gtrsim\pi\Omega$) (Bontz & Haugan, 1981;
Nakamura & Deguchi, 1999; De Paolis et al., 2002), as can be seen by the
flattening of the wave optics magnification curves at small $y$ in the bottom
panel of Fig. 8 (for compact sources, magnifications of up to
$|F|^{2}\approx\pi\Omega$ can be achieved). However, for more modest
amplification factors, the sidelobes of $F(\Omega,y)$ contribute non-
negligibly to the cross-section. Let us take a graphical approach to
understand this. Note the rapidly-varying oscillations in the top panel of
Fig. 8. In a typical optical depth calculation with a certain sensitivity
threshold $|F|^{2}_{\textrm{min}}$, we integrate over all $y$ in the source
plane for which $|F(y)|^{2}>|F|^{2}_{\textrm{min}}$ to obtain the optical
depth $\sigma(>|F_{\textrm{min}}|^{2})$. In the geometric optics case,
$|F_{\textrm{geo}}|^{2}(y)$ is a decreasing function. Therefore, the region of
the source plane that counts towards the optical depth is a circle of radius
$y_{*}$ such that $|F_{\textrm{geo}}(y_{*})|^{2}=|F|^{2}_{\textrm{min}}$. Its
area, the cross-section to lensing, is
$\sigma(>|F_{\textrm{min}}|^{2})=\sigma_{\textrm{geo}}=\pi y_{*}^{2}$. In the
wave optics case, it is not as straightforward to calculate a threshold value
of $y$, because the amplification factor is oscillatory, as seen in the top
panel of Fig. 8. We numerically evaluate Eq. 47 on a dense grid of points
$\\{y_{i}\\}$, checking for each $y_{i}$ whether
$|F(y_{i})|^{2}>|F_{\textrm{min}}|^{2}$ and accumulating $2\pi y_{i}\Delta
y_{i}$ to $\sigma(>|F_{\textrm{min}}|^{2})$ if the magnification is above the
threshold. Repeating this as a function of $|F_{\textrm{min}}|^{2}$, for
$|F_{\textrm{geo}}|^{2}$ (Eq. 57), as well as $|F_{\textrm{wave}}|^{2}$ (we
use the approximation in Eq. 47) allows us to calculate the optical depth in
wave optics and in geometric optics as a function of the sensitivity cutoff
$|F_{\textrm{min}}|^{2}$; the ratio is plotted in Fig. 9 (see also Appendix B
of (Sammons et al., 2022), which studies the effect of high magnification
lensing events on fast transient luminosity functions).

Figure 9: The lensing cross-section $\sigma$ for a minimum detectable
amplification factor $|F|^{2}$. Top: we plot several values of $\Omega$ as
well as the cross-section calculated in geometric optics. Bottom: the ratio of
$\sigma_{\textrm{wave}}/\sigma_{\textrm{geo}}$ is well-approximated by a
boxcar between $|F|^{2}=[1,2\pi\Omega]$ with a height of $y=3/2$. The
conclusion is that interference effects enhances lensing cross-sections for
high magnifications by a factor of $\approx 3/2$ while cutting off huge
amplifications greater than $\approx\pi\Omega$. This result has the caveat
that Eq. 47 is used to approximate Eq. 42, which is an excellent approximation
for $\Omega\gtrsim 1$.

Finally, for $\Omega\lesssim 1$, the lensing is somewhat suppressed. This was
considered in the context of kilohertz gravitational waves, which do not
“feel” mass distributions smaller than $\sim 10^{2}M_{\odot}$ (Nakamura,
1998). The bad news is that this limits the mass reach of searches for
gravitational wave lensing. The good news is that existing search pipelines
are reasonably sensitive to diffractively-lensed gravitational waves, and that
the subtle distortions due to diffractive lensing will be detectable for
louder events (Dai et al., 2018). Moreover, the chromaticity induces a mass
scale dependence of the lensing effect. This enables lensing magnifications to
probe intervening density fluctuations as a function of mass scale (i.e. the
matter power spectrum) (Takahashi, 2006; Oguri & Takahashi, 2020) and dark
matter on dwarf galaxy scales (Gil Choi et al., 2021; Guo & Lu, 2022). In
addition, diffraction suppresses the effect of microlensing on the measured
magnification ratios of gravitational waves macro-lensed by galaxies or
clusters, unless the macro-magnifications $\gtrsim 15$ (Mishra et al., 2021).
This can offer cleaner time-delay cosmography measurements and allow
gravitational waves to be localized (Cheung et al., 2021).

Since $F(\Omega,\bm{y})$ is oscillatory as a function of $\Omega$, the cross-
section of microlensing is higher at low frequencies. However, as the
frequency approaches $\Omega\leq 1$, no microlensing occurs (see next
section). According to Fig. 2, the reach of optical microlensing is also cut
off at low masses due to diffractive effects. These two effects lead to
overestimating cross-sections calculated in the geometric optics limit. Due to
finite source size effects (e.g. (Witt & Mao, 1994)) and wave effects, the
maximum magnification is less than infinity but drops off more slowly. See
e.g. Fig. 5 of (Sugiyama et al., 2020) for a detailed study of the effect of
finite source size effect on the light curve of a transiting optical source.

##  11 Observational prospects and discussion

We have presented a summary of the various aspects of interference in
gravitational lensing. In some ways, interference in gravitational lensing
makes lensing easier to detect. For example, in the Eikonal limit, the
presence of interference terms between images makes searching for the
interference terms much more sensitive than searches for individual images:
one simply searches for a regular sinusoidal pattern in the spectrum whose
spacing allows for a direct measurement of the lensing time delay. For
diffractive optics with a point-mass lens, the sinusoidal pattern becomes
distorted at low frequencies as $\Omega$ tends to 1; the unique form of this
spectral distortion can break degeneraces in microlensing, allowing for a
measurement of the impact parameter and (redshifted) mass of the lens
simultaneously.

In addition, diffraction limits the maximum attainable magnification near
caustics, where interference effects become the most pronounced. We find in
diffractive lensing that extremely large magnifications near cusps and folds
are no longer possible, but that “sidelobes” in the magnification pattern can
boost the cross-section to lensing. We calculate this boost, and find (in
agreement with (Sammons et al., 2022), which appeared while this manuscript
was being prepared) that to a very good approximation, the lensing cross-
section is higher by a factor of $3/2$ as compared to geometric optics for
large ($1<<|F|^{2}<\pi\Omega$) magnifications.

Finally, we characterized the primary challenge to observing interference in
gravitational lensing: the paucity of point-like, or angularly-coherent,
sources of radiation. A true point source, even one that does not exhibit any
significant light curve variations, could be used to measure lensing delays to
Nyquist-limited time resolution (nanoseconds at radio frequencies). Such
compact steady sources are not yet known, but compact transient sources at
cosmological distances such as supernovae and fast radio bursts could be used.
Gravitational waves from compact binary mergers observed by current and
upcoming ground- and space-based detectors are particularly promising
candidates for observing diffractive lensing. We comment briefly on the
different observational prospects here, but refer the reader to (Oguri, 2019)
for a comprehensive discussion of transient lensing.

In the field of lensed transients, supernovae have a head start since a
handful of lensed supernovae have already been discovered (e.g. (Kelly et al.,
2015)). Since $\Omega$ is large for optical and near-infrared supernovae, they
can potentially be extremely highly magnified (Diego, 2019). However,
according to Fig. 8, supernovae photospheres are not sufficiently compact,
except at cosmological distances.

Gamma-ray bursts (GRBs) are among the most luminous and distant transients
known, leading to favorable lensing rates and femtolensing searches over the
last few decades (Gould, 1992), but because of the extremely short wavelengths
involved, their angular sizes are not negligible (Stanek et al., 1993; Witt &
Mao, 1994; Matsunaga & Yamamoto, 2006). While the apparent size of the
emission is dependent on the beaming factor (Katz et al., 2018; Khangulyan et
al., 2022), and can vary from burst to burst, most GRBs are not point-like
enough to observe the spectral features of femtolensing; furthermore the range
of lens sizes required to observe interference effects in GRBs may not be
astrophysical (see Fig. 2).

FRBs are extremely compact, but like active galactic nuclei, they may be
angularly broadened (Rickett, 1990; Lazio et al., 2008) into a “scattering
disk” at the low (sub-GHz) radio frequencies at which they are typically
observed(Katz et al., 2020). Nevertheless, angular coherence in the sense of
Eq. 50 is preserved for sub-solar mass lenses (Kader et al., 2022; Leung et
al., 2022). The bulk of FRBs detected to date have been surveyed by the
Canadian Hydrogen Intensity Mapping Experiment at low (sub-gigahertz)
frequencies (CHIME/FRB Collaboration et al., 2021), but because of the strong
frequency dependence of the size of the scattering disk
($R\propto\lambda^{2}$), these FRBs are sensitive to angular broadening by
inhomogeneous plasma, which potentially washes out interference fringes. This
problem can be mitigated with FRBs surveyed at even slightly higher
frequencies by upcoming experiments like CHORD and the DSA-2000 (Vanderlinde
et al., 2019; Hallinan et al., 2019); these may enable observation of coherent
interference patterns in FRBs lensed by stars (Connor & Ravi, 2022).

Perhaps the most promising way to detect diffractive effects is via
gravitational wave detectors like LIGO (LIGO Scientific Collaboration et al.,
2015), VIRGO (Acernese et al., 2015), and KAGRA (Kagra Collaboration et al.,
2019), whose pristine sources are not contaminated by propagation effects,
except perhaps by microlenses near the images of a macro-lens (Mishra et al.,
2021). Gravitational-wave sources are extremely compact
($\sim$100\text{\,}\mathrm{km}$$; see (Oguri, 2019)), and searches have
already been designed and conducted in both the Eikonal and diffractive
regimes (Haris et al., 2018; Hannuksela et al., 2019; Dai et al., 2020; Abbott
et al., 2021b), though no candidates have been conclusively confirmed. The
observational challenge of GW lensing is that with the current generation of
gravitational wave detectors (LIGO, VIRGO, KAGRA, and GEO600), gravitational
waves are difficult to localize to a host galaxy without an electromagnetic
counterpart. However, rapid progress, particularly in FRBs and gravitational-
wave instruments, will soon lead to detections of interference and diffractive
effects in gravitational lensing, and the creative application of lensing to
unique measurements for astronomy and fundamental physics.

###### Acknowledgements.

We thank the International Space Science Institute in Bern (ISSI) for their
hospitality and the conveners for organizing the workshop on Strong
Gravitational Lensing where this paper was initially written. C.L. was
supported by the U.S. Department of Defense (DoD) through the National Defense
Science & Engineering Graduate Fellowship (NDSEG) Program. L.D. acknowledges
research grant support from the Alfred P. Sloan Foundation (Award Number
FG-2021-16495). M.O. was supported by JSPS KAKENHI Grant Number JP22K18720. We
thank Paul Schechter for helpful discussions throughout the preparation of
this manuscript.

## References

  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102, doi: 10.1103/PhysRevLett.116.061102
  * Abbott et al. (2017) —. 2017, Phys. Rev. Lett., 119, 161101, doi: 10.1103/PhysRevLett.119.161101
  * Abbott et al. (2021a) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021a, ApJ, 915, L5, doi: 10.3847/2041-8213/ac082e
  * Abbott et al. (2021b) —. 2021b, ApJ, 923, 14, doi: 10.3847/1538-4357/ac23db
  * Acernese et al. (2015) Acernese, F., Agathos, M., Agatsuma, K., et al. 2015, Classical and Quantum Gravity, 32, 024001, doi: 10.1088/0264-9381/32/2/024001
  * Asada & Kasai (2000) Asada, H., & Kasai, M. 2000, Progress of Theoretical Physics, 104, 95, doi: 10.1143/PTP.104.95
  * Baraldo et al. (1999) Baraldo, C., Hosoya, A., & Nakamura, T. T. 1999, Phys. Rev. D, 59, 083001, doi: 10.1103/PhysRevD.59.083001
  * Beach & Lovelace (1997) Beach, T. L., & Lovelace, R. V. E. 1997, Radio Science, 32, 913, doi: 10.1029/97RS00063
  * Benson & Cooke (1979) Benson, J. R., & Cooke, J. H. 1979, ApJ, 227, 360, doi: 10.1086/156739
  * Berry (2021) Berry, M. V. 2021, Journal of Optics, 23, 065604, doi: 10.1088/2040-8986/abfee1
  * Berry & Upstill (1980) Berry, M. V., & Upstill, C. 1980, Progress in Optics, 18, 257
  * Blandford & Narayan (1986) Blandford, R., & Narayan, R. 1986, ApJ, 310, 568, doi: 10.1086/164709
  * Bontz & Haugan (1981) Bontz, R. J., & Haugan, M. P. 1981, Ap&SS, 78, 199, doi: 10.1007/BF00654034
  * Boyle & Pen (2012) Boyle, L., & Pen, U.-L. 2012, Phys. Rev. D, 86, 124028, doi: 10.1103/PhysRevD.86.124028
  * Chang & Refsdal (1979) Chang, K., & Refsdal, S. 1979, Nature, 282, 561, doi: 10.1038/282561a0
  * Cheung et al. (2021) Cheung, M. H. Y., Gais, J., Hannuksela, O. A., & Li, T. G. F. 2021, MNRAS, 503, 3326, doi: 10.1093/mnras/stab579
  * CHIME/FRB Collaboration et al. (2021) CHIME/FRB Collaboration, Amiri, M., Andersen, B. C., et al. 2021, ApJS, 257, 59, doi: 10.3847/1538-4365/ac33ab
  * Chwolson (1924) Chwolson, O. 1924, Astronomische Nachrichten, 221, 329, doi: 10.1002/asna.19242212003
  * Congedo et al. (2006) Congedo, G., de Paolis, F., Longo, P., et al. 2006, International Journal of Modern Physics D, 15, 1937, doi: 10.1142/S0218271806009248
  * Connor & Ravi (2022) Connor, L., & Ravi, V. 2022, arXiv e-prints, arXiv:2206.14310.  https://arxiv.org/abs/2206.14310
  * Cremonese et al. (2021) Cremonese, P., Ezquiaga, J. M., & Salzano, V. 2021, Phys. Rev. D, 104, 023503, doi: 10.1103/PhysRevD.104.023503
  * Dai (2014) Dai, L. 2014, Phys. Rev. Lett., 112, 041303, doi: 10.1103/PhysRevLett.112.041303
  * Dai et al. (2018) Dai, L., Li, S.-S., Zackay, B., Mao, S., & Lu, Y. 2018, Phys. Rev. D, 98, 104029, doi: 10.1103/PhysRevD.98.104029
  * Dai & Miralda-Escudé (2020) Dai, L., & Miralda-Escudé, J. 2020, AJ, 159, 49, doi: 10.3847/1538-3881/ab5e83
  * Dai & Venumadhav (2017) Dai, L., & Venumadhav, T. 2017, arXiv e-prints, arXiv:1702.04724.  https://arxiv.org/abs/1702.04724
  * Dai et al. (2020) Dai, L., Zackay, B., Venumadhav, T., Roulet, J., & Zaldarriaga, M. 2020, arXiv e-prints, arXiv:2007.12709.  https://arxiv.org/abs/2007.12709
  * Davydov (2013) Davydov, A. S. 2013, Quantum mechanics: international series in natural philosophy, Vol. 1 (Elsevier)
  * De Paolis et al. (2002) De Paolis, F., Ingrosso, G., Nucita, A. A., & Qadir, A. 2002, A&A, 394, 749, doi: 10.1051/0004-6361:20021258
  * Deguchi & Watson (1986a) Deguchi, S., & Watson, W. D. 1986a, Phys. Rev. D, 34, 1708, doi: 10.1103/PhysRevD.34.1708
  * Deguchi & Watson (1986b) —. 1986b, ApJ, 307, 30, doi: 10.1086/164389
  * Diego (2019) Diego, J. M. 2019, A&A, 625, A84, doi: 10.1051/0004-6361/201833670
  * Doyle & Carico (2009) Doyle, L. R., & Carico, D. P. 2009, The Open Astronomy Journal, 2, 63, doi: 10.2174/1874381100902010063
  * Engeli & Saha (2022) Engeli, S., & Saha, P. 2022, MNRAS, 516, 4679, doi: 10.1093/mnras/stac2522
  * Ezquiaga et al. (2021) Ezquiaga, J. M., Holz, D. E., Hu, W., Lagos, M., & Wald, R. M. 2021, Phys. Rev. D, 103, 064047, doi: 10.1103/PhysRevD.103.064047
  * Ezquiaga et al. (2020) Ezquiaga, J. M., Hu, W., & Lagos, M. 2020, Phys. Rev. D, 102, 023531, doi: 10.1103/PhysRevD.102.023531
  * Falco et al. (1985) Falco, E. E., Gorenstein, M. V., & Shapiro, I. I. 1985, ApJ, 289, L1, doi: 10.1086/184422
  * Feldbrugge et al. (2019) Feldbrugge, J., Pen, U.-L., & Turok, N. 2019, arXiv e-prints, arXiv:1909.04632.  https://arxiv.org/abs/1909.04632
  * Feldbrugge & Turok (2020) Feldbrugge, J., & Turok, N. 2020, arXiv e-prints, arXiv:2008.01154.  https://arxiv.org/abs/2008.01154
  * Gil Choi et al. (2021) Gil Choi, H., Park, C., & Jung, S. 2021, arXiv e-prints, arXiv:2103.08618.  https://arxiv.org/abs/2103.08618
  * Goldreich & Keeley (1972) Goldreich, P., & Keeley, D. A. 1972, ApJ, 174, 517, doi: 10.1086/151514
  * Goodman et al. (1987) Goodman, J. J., Romani, R. W., Blandford, R. D., & Narayan, R. 1987, MNRAS, 229, 73, doi: 10.1093/mnras/229.1.73
  * Gosselin et al. (2007) Gosselin, P., Bérard, A., & Mohrbach, H. 2007, Phys. Rev. D, 75, 084035, doi: 10.1103/PhysRevD.75.084035
  * Gould (1992) Gould, A. 1992, ApJ, 386, L5, doi: 10.1086/186279
  * Grillo & Cordes (2018) Grillo, G., & Cordes, J. 2018, arXiv e-prints, arXiv:1810.09058.  https://arxiv.org/abs/1810.09058
  * Guo & Lu (2022) Guo, X., & Lu, Y. 2022, Phys. Rev. D, 106, 023018, doi: 10.1103/PhysRevD.106.023018
  * Hallinan et al. (2019) Hallinan, G., Ravi, V., Weinreb, S., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 255, doi: 10.48550/arXiv.1907.07648
  * Hannuksela et al. (2019) Hannuksela, O. A., Haris, K., Ng, K. K. Y., et al. 2019, ApJ, 874, L2, doi: 10.3847/2041-8213/ab0c0f
  * Haris et al. (2018) Haris, K., Mehta, A. K., Kumar, S., Venumadhav, T., & Ajith, P. 2018, arXiv e-prints, arXiv:1807.07062.  https://arxiv.org/abs/1807.07062
  * Herlt & Stephani (1976) Herlt, E., & Stephani, H. 1976, International Journal of Theoretical Physics, 15, 45, doi: 10.1007/BF01807086
  * Hewitt et al. (1988) Hewitt, J. N., Turner, E. L., Schneider, D. P., Burke, B. F., & Langston, G. I. 1988, Nature, 333, 537, doi: 10.1038/333537a0
  * Inamori & Suyama (2021) Inamori, M., & Suyama, T. 2021, ApJ, 918, L30, doi: 10.3847/2041-8213/ac2142
  * Ishihara et al. (1988) Ishihara, H., Takahashi, M., & Tomimatsu, A. 1988, Phys. Rev. D, 38, 472, doi: 10.1103/PhysRevD.38.472
  * Itoh et al. (2009) Itoh, Y., Futamase, T., & Hattori, M. 2009, Phys. Rev. D, 80, 044009, doi: 10.1103/PhysRevD.80.044009
  * Jacques et al. (2007) Jacques, V., Wu, E., Grosshans, F., et al. 2007, Science, 315, 966, doi: 10.1126/science.1136303
  * Janquart et al. (2021) Janquart, J., Seo, E., Hannuksela, O. A., Li, T. G. F., & Van Den Broeck, C. 2021, ApJ, 923, L1, doi: 10.3847/2041-8213/ac3bcf
  * Jaroszynski & Paczynski (1995) Jaroszynski, M., & Paczynski, B. 1995, ApJ, 455, 443, doi: 10.1086/176593
  * Johnston et al. (1997) Johnston, K. J., Gaume, R. A., Wilson, T. L., Nguyen, H. A., & Nedoluha, G. E. 1997, ApJ, 490, 758, doi: 10.1086/304891
  * Jow et al. (2020) Jow, D. L., Foreman, S., Pen, U.-L., & Zhu, W. 2020, MNRAS, 497, 4956, doi: 10.1093/mnras/staa2230
  * Jow et al. (2021) Jow, D. L., Lin, F. X., Tyhurst, E., & Pen, U.-L. 2021, MNRAS, 507, 5390, doi: 10.1093/mnras/stab2337
  * Jow et al. (2022) Jow, D. L., Pen, U.-L., & Feldbrugge, J. 2022, arXiv e-prints, arXiv:2204.12004.  https://arxiv.org/abs/2204.12004
  * Kader et al. (2022) Kader, Z., Leung, C., Dobbs, M., et al. 2022, arXiv e-prints, arXiv:2204.06014.  https://arxiv.org/abs/2204.06014
  * Kagra Collaboration et al. (2019) Kagra Collaboration, Akutsu, T., Ando, M., et al. 2019, Nature Astronomy, 3, 35, doi: 10.1038/s41550-018-0658-y
  * Katz et al. (2018) Katz, A., Kopp, J., Sibiryakov, S., & Xue, W. 2018, J. Cosmology Astropart. Phys, 2018, 005, doi: 10.1088/1475-7516/2018/12/005
  * Katz et al. (2020) —. 2020, MNRAS, 496, 564, doi: 10.1093/mnras/staa1497
  * Kelly et al. (2015) Kelly, P. L., Rodney, S. A., Treu, T., et al. 2015, Science, 347, 1123, doi: 10.1126/science.aaa3350
  * Khangulyan et al. (2022) Khangulyan, D., Barkov, M. V., & Popov, S. B. 2022, ApJ, 927, 2, doi: 10.3847/1538-4357/ac4bdf
  * Lazio et al. (2008) Lazio, T. J. W., Ojha, R., Fey, A. L., et al. 2008, ApJ, 672, 115, doi: 10.1086/520572
  * Leung et al. (2018) Leung, C., Brown, A., Nguyen, H., et al. 2018, Phys. Rev. A, 97, 042120, doi: 10.1103/PhysRevA.97.042120
  * Leung et al. (2022) Leung, C., Kader, Z., Masui, K. W., et al. 2022, arXiv e-prints, arXiv:2204.06001.  https://arxiv.org/abs/2204.06001
  * Li et al. (2019) Li, D., Lin, F. X., Main, R., et al. 2019, MNRAS, 484, 5723, doi: 10.1093/mnras/stz374
  * Li et al. (2022) Li, Z., Qiao, J., Zhao, W., & Er, X. 2022, J. Cosmology Astropart. Phys, 2022, 095, doi: 10.1088/1475-7516/2022/10/095
  * Lightman et al. (1975) Lightman, A. P., Press, W. H., Price, R. H., & Teukolsky, S. A. 1975, Problem Book in Relativity and Gravitation
  * LIGO Scientific Collaboration et al. (2015) LIGO Scientific Collaboration, Aasi, J., Abbott, B. P., et al. 2015, Classical and Quantum Gravity, 32, 074001, doi: 10.1088/0264-9381/32/7/074001
  * Lin et al. (2022) Lin, R., van Kerkwijk, M. H., Main, R., et al. 2022, arXiv e-prints, arXiv:2211.05209, doi: 10.48550/arXiv.2211.05209
  * Maccone (2010) Maccone, C. 2010, Acta Astronautica, 67, 521, doi: 10.1016/j.actaastro.2010.03.013
  * Macquart (2004) Macquart, J. P. 2004, A&A, 422, 761, doi: 10.1051/0004-6361:20034512
  * Matsunaga & Yamamoto (2006) Matsunaga, N., & Yamamoto, K. 2006, J. Cosmology Astropart. Phys, 2006, 023, doi: 10.1088/1475-7516/2006/01/023
  * Miller & Wheeler (1997) Miller, W. A., & Wheeler, J. A. 1997, Foundations Of Quantum Mechanics In The Light Of New Technology. Series: Advanced Series in Applied Physics, 4, 72, doi: 10.1142/9789812819895_0008
  * Mishra et al. (2021) Mishra, A., Meena, A. K., More, A., Bose, S., & Bagla, J. S. 2021, MNRAS, 508, 4869, doi: 10.1093/mnras/stab2875
  * Nakamura (1998) Nakamura, T. T. 1998, Phys. Rev. Lett., 80, 1138, doi: 10.1103/PhysRevLett.80.1138
  * Nakamura & Deguchi (1999) Nakamura, T. T., & Deguchi, S. 1999, Progress of Theoretical Physics Supplement, 133, 137, doi: 10.1143/PTPS.133.137
  * Narayan (1992) Narayan, R. 1992, Philosophical Transactions of the Royal Society of London Series A, 341, 151, doi: 10.1098/rsta.1992.0090
  * Nemiroff et al. (1993) Nemiroff, R. J., Norris, J. P., Wickramasinghe, W. A. D. T., et al. 1993, ApJ, 414, 36, doi: 10.1086/173054
  * Nouri-Zonoz (1999) Nouri-Zonoz, M. 1999, Phys. Rev. D, 60, 024013, doi: 10.1103/PhysRevD.60.024013
  * Nye (1999) Nye, J. F. 1999, Natural focusing and fine structure of light: caustics and wave dislocations
  * Oancea et al. (2020) Oancea, M. A., Joudioux, J., Dodin, I. Y., et al. 2020, Phys. Rev. D, 102, 024075, doi: 10.1103/PhysRevD.102.024075
  * Oguri (2019) Oguri, M. 2019, Reports on Progress in Physics, 82, 126901, doi: 10.1088/1361-6633/ab4fc5
  * Oguri & Takahashi (2020) Oguri, M., & Takahashi, R. 2020, ApJ, 901, 58, doi: 10.3847/1538-4357/abafab
  * Oguri & Takahashi (2022) —. 2022, Phys. Rev. D, 106, 043532, doi: 10.1103/PhysRevD.106.043532
  * Ohanian (1973) Ohanian, H. C. 1973, Phys. Rev. D, 8, 2734, doi: 10.1103/PhysRevD.8.2734
  * Ohanian (1974) —. 1974, International Journal of Theoretical Physics, 9, 425, doi: 10.1007/BF01810927
  * Ohanian (1983) —. 1983, ApJ, 271, 551, doi: 10.1086/161221
  * Pearson et al. (2014) Pearson, J. W., Olver, S., & Porter, M. A. 2014, arXiv e-prints, arXiv:1407.7786, doi: 10.48550/arXiv.1407.7786
  * Peters (1974) Peters, P. C. 1974, Phys. Rev. D, 9, 2207, doi: 10.1103/PhysRevD.9.2207
  * Rahvar (2018) Rahvar, S. 2018, MNRAS, 479, 406, doi: 10.1093/mnras/sty1369
  * Rickett (1990) Rickett, B. J. 1990, ARA&A, 28, 561, doi: 10.1146/annurev.aa.28.090190.003021
  * Saha (2000) Saha, P. 2000, AJ, 120, 1654, doi: 10.1086/301581
  * Sammons et al. (2022) Sammons, M. W., James, C. W., Trott, C. M., & Walker, M. 2022, MNRAS, doi: 10.1093/mnras/stac3013
  * Schneider (1983) Schneider, P. 1983, in Liege International Astrophysical Colloquia, Vol. 24, Liege International Astrophysical Colloquia, ed. J.-P. Swings, 131–133
  * Shi & Xu (2021) Shi, X., & Xu, Z. 2021, MNRAS, 506, 6039, doi: 10.1093/mnras/stab2108
  * Stanek et al. (1993) Stanek, K. Z., Paczynski, B., & Goodman, J. 1993, ApJ, 413, L7, doi: 10.1086/186946
  * Stewart (1989) Stewart, J. M. 1989, Proceedings of the Royal Society of London Series A, 424, 239, doi: 10.1098/rspa.1989.0078
  * Sugiyama et al. (2020) Sugiyama, S., Kurita, T., & Takada, M. 2020, MNRAS, 493, 3632, doi: 10.1093/mnras/staa407
  * Suyama (2020) Suyama, T. 2020, ApJ, 896, 46, doi: 10.3847/1538-4357/ab8d3f
  * Suyama et al. (2006) Suyama, T., Tanaka, T., & Takahashi, R. 2006, Phys. Rev. D, 73, 024026, doi: 10.1103/PhysRevD.73.024026
  * Takahashi (2006) Takahashi, R. 2006, ApJ, 644, 80, doi: 10.1086/503323
  * Takahashi (2017) —. 2017, ApJ, 835, 103, doi: 10.3847/1538-4357/835/1/103
  * Takahashi & Nakamura (2003) Takahashi, R., & Nakamura, T. 2003, ApJ, 595, 1039, doi: 10.1086/377430
  * Takahashi et al. (2005) Takahashi, R., Suyama, T., & Michikoshi, S. 2005, A&A, 438, L5, doi: 10.1051/0004-6361:200500140
  * Tanaka & Suyama (2023) Tanaka, S., & Suyama, T. 2023, arXiv e-prints, arXiv:2303.05650.  https://arxiv.org/abs/2303.05650
  * Thom (1967) Thom, R. 1967, Structural Stability And Morphogenesis (CRC Press).  https://books.google.ca/books?id=YVNPDwAAQBAJ
  * Turyshev (2017) Turyshev, S. G. 2017, Phys. Rev. D, 95, 084041, doi: 10.1103/PhysRevD.95.084041
  * Ulmer & Goodman (1995) Ulmer, A., & Goodman, J. 1995, ApJ, 442, 67, doi: 10.1086/175422
  * Vanderlinde et al. (2019) Vanderlinde, K., Liu, A., Gaensler, B., et al. 2019, in Canadian Long Range Plan for Astronomy and Astrophysics White Papers, Vol. 2020, 28, doi: 10.5281/zenodo.3765414
  * Venumadhav et al. (2017) Venumadhav, T., Dai, L., & Miralda-Escudé, J. 2017, ApJ, 850, 49, doi: 10.3847/1538-4357/aa9575
  * Vijaykumar et al. (2022) Vijaykumar, A., Mehta, A. K., & Ganguly, A. 2022, arXiv e-prints, arXiv:2202.06334.  https://arxiv.org/abs/2202.06334
  * Walker et al. (2004) Walker, M. A., Melrose, D. B., Stinebring, D. R., & Zhang, C. M. 2004, MNRAS, 354, 43, doi: 10.1111/j.1365-2966.2004.08159.x
  * Witt & Mao (1994) Witt, H. J., & Mao, S. 1994, ApJ, 430, 505, doi: 10.1086/174426
  * Wucknitz et al. (2021) Wucknitz, O., Spitler, L. G., & Pen, U. L. 2021, A&A, 645, A44, doi: 10.1051/0004-6361/202038248
  * Xiao et al. (2022) Xiao, H., Dai, L., & McQuinn, M. 2022, arXiv e-prints, arXiv:2206.13534.  https://arxiv.org/abs/2206.13534
  * Yamamoto (2003) Yamamoto, K. 2003, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/0309696
  * Yang et al. (2014) Yang, H., Zhang, F., Zimmerman, A., & Chen, Y. 2014, Phys. Rev. D, 89, 064014, doi: 10.1103/PhysRevD.89.064014
  * Zenginoǧlu & Galley (2012) Zenginoǧlu, A., & Galley, C. R. 2012, Phys. Rev. D, 86, 064030, doi: 10.1103/PhysRevD.86.064030