The pixel corresponding to the spectral channel centered on $\overline{k}$
subtends an angle $\overline{\theta}_{\mathrm{F}}/\mathcal{R}$ to capture
wavenumbers in the range
$k\in\left[\bar{k}-\sigma_{k}/2,\bar{k}+\sigma_{k}/2\right]$, but it also
needs to be sufficiently large to contain the entire seeing disk of the
source(s). Atmospheric aberrations cause a spatial spread (which furthermore
fluctuates on millisecond time scales) on the order of the magnified seeing
angle $\theta_{\mathrm{seeing}}f/f^{\prime}$. The illuminated area of the
grating has a linear size $Nd_{\mathrm{g}}$ which is related to the aperture
diameter $D$ of the telescope by the same (inverse) magnification factor (the
focal length ratio of primary and secondary mirrors):
$Nd_{\mathrm{g}}=Df^{\prime}/f$. These relations imply that containment of the
seeing disk is equivalent to:

$\displaystyle\frac{\overline{\theta}_{\mathrm{F}}}{\mathcal{R}}\gtrsim\frac{f}{f^{\prime}}\theta_{\mathrm{seeing}}\quad\Leftrightarrow\quad\frac{N}{\mathcal{R}}\gtrsim\frac{\theta_{\mathrm{seeing}}}{\theta_{\mathrm{diff}}}\approx
40\left(\frac{D}{4\,\mathrm{m}}\right)\left(\frac{500\,\mathrm{nm}}{\overline{\lambda}}\right)\left(\frac{\theta_{\text{seeing}}}{1\,\mathrm{arcsec}}\right),$
(3.11)

which is nothing more than the usual requirement in spectroscopy that gratings
have to be oversized in practice to account for seeing, since atmospheric
fluctuations reduce the grating’s nominal resolution $N$ to the much smaller
effective resolution $\mathcal{R}$. When this inequality is
saturated—$N/\mathcal{R}\simeq\theta_{\mathrm{seeing}}/\theta_{\mathrm{diff}}$—then
the suppression found in eq. 3.10 is equal to:

$\displaystyle
C\leq\frac{\pi\sigma_{d}}{\sigma_{k}}\simeq\frac{\theta_{\mathrm{diff}}}{\theta_{\mathrm{seeing}}},$
(3.12)

i.e. the parametric upper bound stated in eq. 3.5 in the beginning of the
section. In figure 12, we show that Phase I of EPIC would already be limited
by (and saturate) eq. 3.12, and thus not achieve its theoretical SNR without
modification; the suppression would be even worse for Phases II and III.

#### Potential solutions

There are at least two (and possibly more) straightforward solutions to
mitigate this spectroscopic suppression, albeit at the expense of additional
complexity and cost.

The first and simplest strategy is to use multiple $N_{\mathrm{subap}}$
subapertures in each EPIC telescope, since the correlation ceiling from eq.
3.12 scales as $C\propto\theta_{\mathrm{diff}}$ and thus with the inverse
linear size of photocollection. One could take the configuration from figure
5(b), and after the path-extension stage, geometrically split the light with
$N_{\mathrm{subap}}$ separate gratings, focusing mirrors, and photodetector
arrays. Each photodetector array would thus have a smaller effective
photocollection size (in the spectroscopic direction), and operate near its
nominal fringe contrast of $C\sim 1/c\sigma_{k}\sigma_{t}$ with the
appropriate subaperture number $N_{\mathrm{subap}}\sim 4,50,100$ for Phases I,
II, III, respectively (see figure 12).333A related technique involves
correlation of sub-pupil intensities to perform aperture synthesis [112]. The
advantages of this approach are that tolerances on the primary and secondary
mirrors are not severe, and that the photodetector reset time constraints are
relaxed by a factor $N_{\mathrm{subap}}^{-1}$. The disadvantage (besides cost
and complexity) is that the photodetector DCR requirements are more stringent
by the same factor $N_{\mathrm{subap}}^{-1}$.

A second potential solution is the employment of adaptive optics (AO), a well-
established technology to recover the diffraction limit of the telescope by
using wavefront sensors and thin deformable mirrors to cancel atmospheric
aberrations. AO is already in use on many spectrographs (see ref. [113] for a
review). In an EPIC context, this would allow for pixel sizes
$\overline{\theta}_{F}/\mathcal{R}\simeq\theta_{\mathrm{diff}}$, for which the
spectroscopic suppression from eq. 3.12 does not occur. The advantage of this
approach is that it can allow for even higher spectral resolutions
($\mathcal{R}\gtrsim 100{,}000$) and correspondingly larger fringe contrast
$C$. This elegant approach would however require tight tolerances on the
polishing of _all_ optically active surfaces, and come at considerable cost
especially for large apertures.

Besides those brute-force approaches, there may exist a method to correct for
these intra-pixel dispersive phases with clever optical configurations,
possibly consisting of arrays of microlenses and/or microprisms, which have
demonstrated lenslet pitch on the order of tens of microns to fit SNSPD or
SPAD pixel sizes [114]. Any such optical elements would have to cancel the
relative propagation phases of a fractional wavenumber range of order
$1/\mathcal{R}$, i.e. “undo” the spectral splitting within the pixel itself so
that $\mathcal{N}\simeq\mathcal{D}$ in the limit $\sigma_{t}\to 0$ in eqs. 3.8
and 3.9.444A minimal configuration would be a diverging lens to recollimate
the beam, a corrective dispersive element such as a microprism, and
potentially a final refocusing lens, all directly above each pixel. A simple
calculation of the requirements of the microprism dictates resolving power
$b\frac{\text{d}n}{\text{d}\lambda}=N$, where $b$ is the base of the prism and
$n$ its index of refraction, which should be accurate to about one part in
$N/N_{\text{eff}}$, where $\sigma_{t}\equiv\bar{k}/{\pi N_{\text{eff}}}$, i.e.
$N_{\text{eff}}\sim 15{,}000$ for a $10\,\mathrm{ps}$ timing resolution. The
feasibility of achieving this on the scale of individual pixels is left to
future study. We have checked that simple modifications—defocusing,
recollimating the light onto the pixel, irregular gratings, or a non-co-phased
segmented mirror—by themselves do not resolve the spectroscopic suppression.
We leave investigations of possible clever optical configurations to future
work.

###  3.5 Tolerances

#### Optical tolerances

One of the attractive features of intensity interferometry is the mildness of
the required tolerances on the telescope optics, in stark contrast to the
tight tolerances for imaging telescopes and (especially optical) amplitude
interferometers. In the simplest realization of intensity interferometry,
without spectral splitting nor path extension, the primary requirement is that
the wavefronts in either telescope do not acquire relative path differences
(or equivalently, time delays) larger than the relative timing resolution
$\sigma_{t}$ of the photodetectors. The phototubes used at the NSII in the
1970s had electronic bandwidths on the order of $60\,\mathrm{MHz}$ [37] or
timing resolutions of $\sigma_{t}\sim 10\,\mathrm{ns}$. This poor timing
translated to a tolerance on the order of a meter on the location of the
receivers, and was therefore clearly not the limiting factor of the
observatory. Even for our Phase III parameter of $\sigma_{t}=3\,\mathrm{ps}$,
significant loss of contrast can avoided as long as there are no relative path
differences greater than a millimeter, about 4–5 orders of magnitude above the
sub-wavelength rms error needed for conventional diffraction-limited
astronomical observations. For the same reason, intensity interferometry is
also quite immune to atmospheric distortions and delays (section 2.5).

The elements with the most stringent tolerances of (multichannel) EPIC fall
into two categories: those that may cause differential phases to the light
from both sources, and those with common effects. The path extension belongs
to the former, as it is designed to (spatially) split the light from both
sources: the two-photon amplitude depicted in the left panel of figure 2
corresponds to the fringe useable for wide-angle differential astrometry (cfr.
eqs. 2.25 and 2.26). Therefore, the relative path length $\ell_{p}$ needs to
be _known_ in each telescope $p=1,2$ to a fraction of a wavelength. In
principle, it needs to be _adjustable_ to within a range of order
$\sigma_{k}^{-1}$ from the optimal values of eqs. 2.25 and 2.26 to maintain
fringe contrast, though it is advisable in practice to keep it well within
that range to avoid fringe confusion. These considerations apply to the two
beamsplitters and two reflecting flat mirrors of the path extension module. We
provide a simple analysis of the beamsplitters in appendix C.2, where we also
elaborate on additional phases from the splitters’ finite width and index of
refraction, which also need to be characterized to achieve optimal operation
of an EPIC interferometer.

Optical elements that affect the light of both sources simultaneously—the
primary and secondary mirrors, diffraction grating, and camera—should be
positioned and polished well enough for the resulting image to fit within the
pixel corresponding to each spectral channel, so that the desired spectral
resolution $\mathcal{R}$ is achieved. This tends to be a significantly
stricter requirement than the aforementioned polishing to within $\sigma_{t}$;
numerically, it amounts to upper bounds on primary aberrations of
$\mathcal{O}(\text{few}\,\mathrm{\mu m})$ for the parameters under
consideration. This is a much looser requirement than on conventional imaging
telescopes aiming to achieve diffraction-limited resolution: tolerances on an
EPIC telescope are looser by a factor of about
$f^{\prime}D/f\overline{\lambda}\mathcal{R}\sim 200$ for $D=10\,\mathrm{m}$,
$\overline{\lambda}=500\,\mathrm{nm}$, $f/f^{\prime}=10$ and
$\mathcal{R}=10{,}000$, the fiducial parameters for EPIC Phase II. We provide
a more detailed account on effects from aberrations in appendix C.3. These
less restrictive tolerances imply that some of the expensive and time-
consuming aspects of conventional large-mirror polishing could be skipped
[115], helping to make the cost-effective construction of large EPIC arrays a
realistic goal.

As we explained in section 3.4, our assumed spectral resolution in combination
with exquisite timing resolution may lead to a suppression in fringe contrast
that has to be mitigated. The tolerance requirements mentioned above are for
optical configurations than can correct for this intra-pixel dispersion. If
one decides to solve the issue raised in section 3.4 via multiple sub-
apertures, then the requirement is that the aberrations yield an image for a
point source that is no larger than the seeing disk. If one employs adaptive
optics instead, then the aberration requirements tighten to those of a
conventional imaging telescope operated at the diffraction limit.

#### Geodety, synchronization, and calibration

Extracting astrometrically useful information (such as the angle
$\boldsymbol{\mathbf{\theta}}_{ba}$ between two sources) from the intensity
correlation $C$ requires excellent geodety and synchronization—a precise
knowledge of the relative positions (i.e. baseline
$\boldsymbol{\mathbf{d}}=\boldsymbol{\mathbf{r}}_{2}-\boldsymbol{\mathbf{r}}_{1}$)
and times (i.e. timing offset $\tau=t_{2}-t_{1}$) at each telescope. These
techniques have been developed in support of VLBI programs over decades and
continue to improve [116, 117].

Geological shifts, tides, and Earth’s precession and nutation can cause drifts
in the relative positions of the telescopes to each other as well as to the
celestial reference frame, typically on the order of centimeters per year
[117]. Accurate geodetic data are thus key for accurate measurements, since
the baseline $\boldsymbol{\mathbf{d}}$ is to be known to a precision of
$c\sigma_{t}$ (about $1\,\mathrm{cm}$ for Phase I, and $1\,\mathrm{mm}$ for
Phase III, cfr. table 1). On the $1\,\mathrm{km}$–$100\,\mathrm{km}$ scale of
the baselines considered here, the requirements are somewhat less stringent
compared to VLBI. In particular, for sources observed near zenith, the
vertical motion of Earth’s surface, which is less precisely known, contributes
subdominantly to the differential astrometric phase of interest
($\overline{k}\boldsymbol{\mathbf{\theta}}_{ba}\cdot\boldsymbol{\mathbf{d}}$).
The VLBI global geodetic measurements are planned to improve by an order of
magnitude compared to current precision as part of the next-generation VLBI
Global Observing System (VGOS) [118, 119], which is currently being
implemented. A geodetic analysis of the data from a recent demonstration of a
subset of the VGOS system yielded a root-mean-square deviation of the baseline
length residuals about the weighted mean of $1.6\,\mathrm{mm}$ [120]. The
expected improvements will result in residual path length errors of order a
millimeter, or timing precision of a few picoseconds across the globe [121],
sufficient for the anticipated tolerances for EPIC.

The fast photon detection and short intensity fluctuation times set stringent
requirements on the timing stability of the telescopes. In particular, the
timing offset between the two telescopes needs to be known to a precision
better than the quadrature sum of the photodetectors’ timing resolutions (cfr.
eq. 2.8) in order to be a subleading cause for fringe contrast loss. A
relative synchronization of $\sigma_{t}\lesssim 10\,\mathrm{ps}$ over an
integration time of $t_{\text{obs}}=10^{4}\,\mathrm{s}$ corresponds to a
required Allan deviation of $10^{-14}$. For shorter baselines, it may be
possible to achieve this by sharing a common frequency reference, but for
longer baselines, each telescope (site) needs to be equipped with its own
high-precision atomic clock. Allan deviations at this level are achievable
with commercially available hydrogen masers, both active [122] and passive
[123], and are employed and tested by e.g. the Event Horizon Telescope [28].

If $\mathrm{SNR}>1$ can be achieved on the intensity correlation fringe on
shorter integration times with large apertures and/or bright sources, the
requirements on the timing stability can be relaxed, making possible the use
of cheaper rubidium [124, 125] or cesium [126] frequency standards, which in
commercial packages can have Allan deviations down to $\mathcal{O}(10^{-13})$.
We note that in addition to the time-delayed (ghost) intensity correlation
fringe at e.g. $\ell=\ell_{1}^{\mathrm{opt}}$ and
$\ell^{\prime}=-\ell_{2}^{\mathrm{opt}}$ most useful for relative astrometry,
one can also employ the “main” intensity correlation fringe in the first line
of eq. 2.24, which has about _four_ times the contrast of the ghost fringe and
can thus be detected _sixteen_ times faster (at the same SNR). This main
fringe can be used to self-calibrate and tease out any drifts in baseline and
time difference, and can therefore ensure the telescopes’ clocks remain
synchronized over the course of the observation.

Finally, each telescope’s optical path and path extension needs to be
calibrated and regularly monitored. This can be done experimentally with a
known light source (e.g. laser) at the other “input ports” on the right of the
beamsplitters in figure 5(b). A calibration should measure (frequency-
dependent) phase shifts from propagation _within_ the beamsplitter, the
distance $w_{p}$, and angles $\gamma$ and $\gamma^{\prime}$. Likewise, the
mapping of pixels to wavelengths needs to be established, but this is a
standard procedure in spectroscopy, and can be performed by observation of a
bright source with a known spectrum.

##  4 Applications

Historically, stellar intensity interferometry has mainly been used to
investigate the morphology of bright emission regions. In particular, one can
determine stellar radii and shapes (projected on the celestial sphere) to high
precision, and study star spots, stellar winds, and limb darkening, which are
all quantified by the form factors $\mathcal{F}(\overline{k},\mathbf{d})$ in
eq. 2.34. For more information, we refer the reader to the classic intensity
interferometry papers on this subject cited in section 1. Emission region
morphology studies of _single sources_ are also possible in a extended path
intensity correlator, _at the same time as differential astrometric
measurements between two sources_. For example, one can collect data for
differential astrometry on two sources $a$ and $b$, and afterwards select the
appropriate fringe by adjusting $\tau$ for the desired application. For
differential astrometry, one would use $\tau$ equal to that of eq. 2.27, while
the emission morphology of source $j=a,b$ could be investigated by taking e.g.
$\tau=-\hat{\boldsymbol{\mathbf{\theta}}}_{j}\cdot\boldsymbol{\mathbf{d}}/c$
(cfr. eq. 2.24), as well as potentially other fringe combinations.

In this section, we present a (partial) science case for the differential
astrometry capabilities of an extended-path intensity correlator (see ref.
[127] for a general overview). The path-extension modification of an intensity
interferometer increases its effective field of view up to a few arcseconds,
limited by atmospheric fluctuations, cfr. figure 9. This (literally and
figuratively) increases the scope of stellar intensity interferometry,
enabling _extreme-precision differential astrometry between two widely
separated sources_ , with the hierarchy:
$\sigma_{\delta\theta}\ll\sigma_{\text{res}}\ll\theta_{ba}\lesssim\theta_{0}$.
In other words, the light-centroiding precision $\sigma_{\delta\theta}$ from
eq. 2.51 is much better than the angular resolution
$\sigma_{\theta_{\text{res}}}$ from eq. 2.9, which can be orders of magnitudes
smaller than the source separation $\theta_{ba}$. The fringe is still detected
at maximum contrast as long as the sources are separated by less than the
isoplanatic angle $\theta_{0}$ from eq. 2.61.

We identify five use cases where differential astrometry between two bright
sources with EPIC has the potential to provide unprecedented accuracy and
precision, opening up new prospects for astrophysical and cosmological
measurements. The five uses cases are: binary-orbit characterization (section
4.1), exoplanet detection (section 4.2), stellar acceleration measurements
(section 4.3), stellar microlensing (section 4.4), and quasar microlensing by
low-mass dark matter halos (section 4.5). This last signature is entirely new,
and will be explored in more detail in a follow-up paper [128]. This list of
applications is almost certainly not exhaustive, but hopefully sufficient to
motivate research and development in support for EPIC Phase I and subsequent
Phases with larger apertures and arrays and/or better spectroscopy and timing.

###  4.1 Binaries

Differential astrometry over a narrow field of view is ideally suited to
characterization of binary orbits, since the two sources are automatically
nearby, and orbital measurements only demand good precision on the _relative_
separation of the sources. Close binary characterization and distance
measurement was already achieved in the early days of intensity interferometry
[42]. Our EPIC setup allows for observation of binary pairs that are much more
widely separated and fainter.

With data only from an intensity interferometer, the period $P$ and all angles
(including inclination) of the binary orbit can be determined. With additional
spectroscopic radial-velocity measurements of each of the binary constituents,
the full Keplerian orbit can be fixed without degeneracies, including the
masses $M_{a}$, $M_{b}$, and the line-of-sight distance $D_{\text{s}}$ to the
binary. We first review some basic notation of binary orbits, and then present
mock analyses on potential binary targets.

Let
$\boldsymbol{\mathbf{r}}\equiv\boldsymbol{\mathbf{r}}_{b}-\boldsymbol{\mathbf{r}}_{a}\equiv\\{x,y,z\\}$
be the relative separation between the two sources in the observer’s frame,
with $z$ denoting the line-of-sight direction. In this frame, we can describe
any Keplerian orbit of period $P$ with the six orbital elements: semi-major
axis $a$, eccentricity $e$, inclination $i$, longitude of ascending node
$\Omega$, argument of periapsis $\omega$, and true anomaly $\nu$:

$\displaystyle x$
$\displaystyle=r\left[\cos\Omega\cos(\omega+\nu)-\sin\Omega\sin(\omega+\nu)\cos
i\right];$ (4.1) $\displaystyle y$
$\displaystyle=r\left[\sin\Omega\cos(\omega+\nu)+\cos\Omega\sin(\omega+\nu)\cos
i\right];$ (4.2) $\displaystyle z$ $\displaystyle=r\left[\sin(\omega+\nu)\sin
i\right].$ (4.3)

The separation $r$ evolves as $a(1-e^{2})/(1+e\cos\nu)$. With that
substitution, all time dependence is encapsulated though the variation of the
true anomaly $\nu$, which is related to the eccentric anomaly $E$, and to the
mean anomaly $\mathcal{M}$ evolving linearly in time $t$ as:

$\displaystyle\nu=2\arctan\left\\{\sqrt{\frac{1+e}{1-e}}\tan\frac{E}{2}\right\\};\quad
E-e\sin E=\mathcal{M}=\frac{2\pi}{P}(t-t_{0}),$ (4.4)

where $t_{0}$ is the reference epoch.

With 7 or more purely astrometric measurements (ideally spaced out roughly
evenly over the orbit) obtained by an intensity interferometer, which can only
measure the projected angles on the sky, namely
$\\{x/D_{\text{s}},y/D_{\text{s}}\\}$, one can fit for the six “dimensionless”
orbital elements $\\{a/D_{\text{s}},e,i,\Omega,\omega,t_{0}/P\\}$ and one
dimensionful quantity, the orbital period $P=2\pi a^{3/2}(GM)^{-1/2}$ with
$M=M_{a}+M_{b}$ the total mass. However, there would always remain one size-
distance-mass degeneracy along the “direction” $a\propto
D_{\text{s}}\propto(GM)^{1/3}$. This final degeneracy can be broken with
radial velocity measurements of the constituents: $\dot{z}_{a}$ and
$\dot{z}_{b}$. Their difference $\dot{z}_{b}-\dot{z}_{a}=\dot{z}$ constrains
the ratio $\sqrt{a/GM}$ independent of $D_{\text{s}}$, while their ratio
$\dot{z}_{b}/\dot{z}_{a}$ (after subtracting out the time-averaged mean radial
velocity of the binary’s center of mass) constrains the mass ratio
$M_{a}/M_{b}$ in direct proportionality.

As a simplified example, we perform an Markov Chain Monte Carlo (MCMC)
analysis on a mock data set reminiscent of the spectroscopic binary W0135-052
composed of an unresolved but detached pair of DA white dwarfs [129, 130]. The
parameters of this system are: binary period $P=1.556\,\mathrm{day}$,
effective temperatures $T_{a}=7450\pm 500\,\mathrm{K}$, $T_{b}=6920\pm
500\,\mathrm{K}$, component masses $M_{a}=0.47\pm 0.05\,M_{\odot}$ and
$M_{b}=0.52\pm 0.05\,M_{\odot}$, velocity semi-amplitudes
$K_{a}=76.5\,\mathrm{km/s}$ and $K_{b}=66.3\,\mathrm{km/s}$, and a parallax-
based distance $D_{\text{s}}=12.35\pm 0.43\,\mathrm{pc}$. The masses are only
estimated after measurements of the surface gravity and assumptions about the
stellar composition and mass-radius relation [130].

Figure 13: Mock data (red/orange points with error bars) for astrometric (left
panel) and spectroscopic (right panel) measurements of the doubly-degenerate
binary WD0135-052A/B. The EPIC Phase I astrometric uncertainties
$\sigma_{\delta\theta}=5.6\times 10^{-11}\,\mathrm{rad}$ are magnified by a
factor of 10 for clarity.

For simplicity, we assumed “truth parameters” for the mock data set of
$P=1.556\,\mathrm{day}$, $D_{\text{s}}=12.35\,\mathrm{pc}$, velocity semi-
amplitudes of $K_{a}=76.5\,\mathrm{km/s}$ and $K_{b}=66.3\,\mathrm{km/s}$,
primary mass $M_{a}=0.47\,M_{\odot}$, secondary mass
$M_{b}=M_{a}K_{a}/K_{b}=0.5423\,M_{\odot}$, ignored the systemic line-of-sight
radial velocity, and assumed a circular orbit ($e=0$, $\omega$ undefined) with
inclination $i=0.886\,\mathrm{rad}=50.8^{\circ}$ and longitude of ascending
node $\Omega=0$. We set $t_{0}=0$, and simulated 40 measurements taken over 40
days (1 day between each epoch), with astrometric precision of
$\sigma_{\delta\theta}=5.6\times 10^{-11}\,\mathrm{rad}$ along a random axis
in $(x/D_{\text{s}},y/D_{\text{s}})$ space, and radial velocity precision of
$\sigma_{\dot{z}_{i}}=1\,\mathrm{km/s}$, both per epoch. This astrometric
precision is in line with the expectation for this binary for EPIC Phase I,
and the spectroscopic precision to that of refs. [131, 132]. We show mock data
for these truth parameters in figure 13.

Figure 14: MCMC posteriors on the parameters
$\\{a,i,\Omega,M_{a},M_{b},D_{\text{s}}\\}$ for a (restricted) binary orbit of
WD0135-052A/B, based on the mock data of figure 13. The truth parameters are
indicated by green lines (and squares). The contours are $n=1,2,3$ sigma
contours (levels $1-e^{-n^{2}/2}$) for the 2D slices, and the vertical black
dashed lines indicate the 16%, 50%, and 84% quantiles (also quoted numerically
in the titles) in the 1D histograms.

In figure 14, we show the results of a MCMC parameter estimation of
$\\{a,i,\Omega,M_{a},M_{b},D_{\text{s}}\\}$ based on the mock data from figure
13. For this system, the “astrometric angles” such as $a/D_{\text{s}}$, $i$,
and $\Omega$ can be determined at $10^{-3}\,\mathrm{rad}$ level of precision,
whereas the dimensionful quantities $M_{i}$, $a$, and $D_{\text{s}}$
separately are slightly less tightly constrained (but still at a precision of
about $0.5\%$ or better fractionally) due to the limited spectroscopic
precision. Nevertheless, it is clear that even with these conservative
assumptions, binary mass and distance measurements should be possible at
unprecedented precision. For example, observations of WD0135-052A/B (or other
similar nearby systems [133]) with results as depicted in figure 14 would
yield valuable constraints on the mass-radius relation of white dwarfs [134,
135, 136], as the same intensity interferometry observations would also yield
both of the white dwarfs’ radii. With EPIC Phase II, differential astrometric
precision of about $\sigma_{\delta\theta}\approx 5.6\times
10^{-12}\,\mathrm{rad}$ should be possible (cfr. figure 8), while the
spectroscopic precision of $\sigma_{\dot{z}_{i}}\sim 10\,\mathrm{m/s}$ should
be possible with dedicated observations by modern instruments such as ESPRESSO
at the VLT [137]. A similar observation run with such capabilities should be
able to characterize all orbital parameters at a fractional precision of
$10^{-4}$.

This type of analysis can give even more striking results on very bright
stars. For example, the spectroscopic binary $\gamma^{2}$ Velorum, comprising
a O-type blue supergiant and a Wolf-Rayet (WR) star at $D_{\text{s}}=379\pm
4\,\mathrm{pc}$. An astrometric precision per epoch of
$\sigma_{\delta\theta}\sim 10^{-13}\,\mathrm{rad}$ should already be
achievable with EPIC Phase I, and radial velocity precision of
$\sigma_{\dot{z}_{i}}\sim 20\,\mathrm{m/s}$ per epoch should also be possible
on such a bright binary, yielding component masses with a fractional precision
at or below the $10^{-4}$ level. At that level (or better), one could
conceivably start to measure the mass loss rate of the WR star at
$10^{-4}\,M_{\odot}/\mathrm{yr}$, and obtain even more precise line-of-sight
distance and semi-major axis. These “spectro-geometric” distance
determinations on binaries involving red giants and/or Cepheid variables may
also help calibrate the cosmic distance ladder independently, though not to
very large line-of-sight distances.

###  4.2 Exoplanets

In the previous section, we have seen that EPIC can be a powerful tool for
characterizing binary orbits, or orbits of more complicated triple/quadruple
systems, comprising _visible_ constituents. However, as long as there are two
visible sources—either a visible binary or a double star from an accidental
chance alignment—one can search for anomalous astrometric deviations due to
invisible dark companions. This includes exoplanet searches, as outlined in
[1]. Depending on the eccentricity and inclination of the exoplanet orbit, the
gravitational attraction from the exoplanet will introduce an additional
astrometric wobble on its host star, which will trace out an ellipse projected
onto the celestial sphere (on top its unperturbed motion). For the special
(but most motivated) case of a circular orbit with mass $M_{\text{p}}$, semi-
major axis $a_{\text{p}}$, inclination $i_{\text{p}}$, period $P_{\text{p}}$,
and reference epoch $t_{0}$, the magnitude of the host star’s astrometric
wobble is:

$\displaystyle\Delta\theta_{a}\simeq\frac{M_{\text{p}}}{M_{a}}\frac{a_{\text{p}}}{D_{a}}\sqrt{1-\sin^{2}i_{\text{p}}\,\sin^{2}\left(2\pi\frac{t-t_{0}}{P_{\text{p}}}\right)}.$
(4.5)

The host star’s mass is $M_{a}$ and its line-of-sight distance $D_{a}$.

The method of looking for anomalous astrometric deviations of stellar light
centroids has not been the most prolific exoplanet discovery technique
historically, but the Gaia satellite has already found dozens of candidates in
early data releases [5], and is expected to find thousands [138] of high-mass
planets, preferentially at large semimajor axis $a_{\text{p}}$ around nearby
stars, due to the scalings of eq. 4.5. In figure 15, we plot the expected
sensitivity of Gaia at its fifth data release (DR5) containing ten years of
observations, for a Sun-like star at $D_{a}=20\,\text{pc}$ in dashed blue. The
sensitivity contour corresponds to a final-mission sensitivity of
$\sigma_{\Delta\theta_{a}}\approx 7\,\mathrm{\mu as}$, and a 3-$\sigma$
detection threshold, for orbits shorter than the mission duration, i.e.
$a_{\text{p}}\lesssim 4.6\,\mathrm{AU}$.

Figure 15: Sensitivity of exoplanet detection techniques as a function of
semimajor axis $a_{\text{p}}$ and mass $M_{\text{p}}$. The astrometric
3-$\sigma$ detection sensitivity of EPIC Phases I-III is shown by the blue
regions for a Sun-like star at $D=20\,\mathrm{pc}$. The markers depict
confirmed exoplanets in multiple-star systems (those where EPIC is definitely
applicable), classified by discovery method from the NASA Exoplanet Archive
[139, acquired May 2023]. The astrometric sensitivity of Gaia DR5 is outlined
by the dashed blue line, that of radial velocity (RV) surveys by the orange
lines, transit methods by the purple lines, direct imaging by the red line,
and microlensing by the green lines (adapted from ref. [140]). The eight Solar
System planets and planets discovered eclipse timing variations (TV) are also
overlaid (image credit: NASA). Figure from companion paper [1].

An EPIC setup can surpass this sensitivity by a couple of orders of magnitude
for nearby stars, or alternatively, extend the volume over which the technique
is applicable. For a pair of Sun-like stars at $D_{a}=20\,\text{pc}$, the per-
epoch light centroiding precision is
$\sigma_{\delta\theta}\approx\\{4.5,0.29,0.011\\}\,\mathrm{\mu as}$ for EPIC
Phases $\\{\mathrm{I,II,III}\\}$. With $N_{\text{obs}}=300$ observations (e.g.
a cadence of about one per month over 30 years), the final wobble sensitivity
is of order
$\sigma_{\Delta\theta_{a}}\approx\\{0.26,0.017,0.00065\\}\,\mathrm{\mu as}$,
indeed better than Gaia by a large factor for this particular application (for
nearby stars, with a companion in close proximity). We show contours of
3-$\sigma$ sensitivity as blue regions in figure 15 for
$D_{a}=20\,\mathrm{pc}$. Earth-Sun-like systems ($M_{a}=M_{\odot}$,
$M_{\text{p}}=M_{\oplus}$, and $a_{\text{p}}=1\,\mathrm{AU}$) should be
observable by EPIC Phase II at distances of $D_{a}\lesssim 20\,\mathrm{pc}$
($D_{a}\lesssim 400\,\mathrm{pc}$ by Phase III), as long as a sufficiently
bright reference star is nearby.

The exoplanet parameter space accessible to an EPIC is complementary to that
of other methods. Exoplanet transits causing repeating, temporary dips in the
host star’s brightness have been detected from the ground [141, 142] and
especially with the Kepler satellite [143], which has produced the largest
exoplanet catalog to date; the approximate parameter space is depicted by
purple lines in figure 15. These surveys have a moderate bias toward small
semimajor axes due to the higher transit frequency. Likewise, many exoplanets
are detected via the radial velocity (RV) they impart on their host star
[144], but this method also favors small semimajor axes, as the radial
velocity scales as
$\dot{z}_{a}\sim\sqrt{GM_{\text{p}}^{2}/M_{a}a_{\text{p}}}$. We plot the RV
detection sensitivity for precisions of $3\,\mathrm{m\,s^{-1}}$ in solid
orange, and $0.3\,\mathrm{m\,s^{-1}}$ in dashed orange. Direct imaging also
favors large planets far away from their host star (red) [145]. Finally,
chance alignments of exoplanetary systems with background stars can lead to
detection of even very low-mass systems; the approximate sensitivity with
ground-based observations is depicted in solid green, and with space-based
observations (e.g. Nancy Grace Roman Space Telescope) in dashed green,
following ref. [146].

###  4.3 Accelerations

The exquisite astrometric light-centroiding of EPIC enables precise
measurements of the relative transverse angular accelerations between two
stars, and thus of the difference in their gravitational potential gradient
(or at least two components thereof). This is a _direct, local_ probe of the
gravitational potential and the underlying matter distribution, unlike
techniques based on the MW’s rotation curve [147] (not local, since they are
based on the enclosed mass within the orbit) or on velocity dispersions [148,
149, 150, 151, 152] (not direct, since they rely on assumptions of equilibrium
and symmetry, which are furthermore violated [153, 154, 155]).

The line-of-sight component of the MW acceleration (relative to that of the
Solar barycenter) can be measured by next-generation spectrographs through the
resulting changes in radial velocity [156, 157, 158] and by pulsar timing
arrays through changes in the pulsar’s spin period [159] or in the orbital
period of pulsar binaries [159, 160, 161]. The latter orbital period drift is
already constraining the vertical acceleration towards the Galactic disk at
moderate signal-to-noise ratio, with further improvements expected over time.
Ref. [15] has proposed using Gaia and future astrometric observatories
(potentially in tandem) to measure the transverse angular acceleration of MW
stars, forecasting that Gaia alone should have significant evidence for the
MW’s transverse angular acceleration by the end of its mission, if large-angle
systematics can be kept under control.

In this section, we revisit the idea from ref. [15] and show that EPIC can
significantly surpass the sensitivity of Gaia by a large factor, albeit on a
much smaller number of stars due to its limitations on apparent brightness.
For a _single_ pair of stars, it may approach a (relative) angular
acceleration sensitivity of $\sigma_{\delta\alpha}\sim 10^{-3}\,\mathrm{\mu
as\,yr^{-2}}$ over a 30-year observing campaign, close to the expected size of
Galactic accelerations.

Figure 16: Apparent angular acceleration vector field in the Milky Way in
galactic coordinates at three line-of-sight distances
$D_{\text{s}}=1\,\mathrm{kpc}$ (top), $D_{\text{s}}=200\,\mathrm{pc}$ (bottom
left), and $D_{\text{s}}=5\,\mathrm{kpc}$ (bottom right). The angular
accelerations towards the Galactic disk are dominant (but no larger than
$10^{-3}\,\mathrm{\mu as}\,\mathrm{yr}^{-2}$, as indicated by the color bars)
for small distances, while for $D_{\text{s}}\gtrsim 1\,\mathrm{kpc}$ the bulge
and DM halo also give important contributions. The nodes at
$(l,b)\approx\\{(0^{\circ},0^{\circ}),(90^{\circ},0^{\circ}),(180^{\circ},0^{\circ}),(270^{\circ},0^{\circ}),(0^{\circ},90^{\circ}),(0^{\circ},-90^{\circ})\\}$
for small $D_{\text{s}}$ are dictated by symmetry.

To estimate the size of the Galactic accelerations, we follow ref. [15] and
use a simple model of the MW matter components, which is a combination of a
bulge, disk, and spherical DM halo. We use the the Hernquist bulge model [162]
with total mass $M_{\text{bulge}}=5\times 10^{9}\,M_{\odot}$ and scale radius
$r_{\text{bulge}}=0.6\,\mathrm{kpc}$, a Miyamoto-Nagai disk model [163] with
mass $M_{\text{disk}}=6.8\times 10^{10}\,M_{\odot}$, scale length
$a_{\text{disk}}=3\,\mathrm{kpc}$, and scale height
$b_{\mathrm{disk}}=280\,\mathrm{pc}$, and the Navarro-Frenk-White (NFW)
profile [164] for the DM halo with scale radius $r_{s}=16\,\mathrm{kpc}$ and
local density (at $r_{\odot}=8\,\mathrm{kpc}$ from the Galactic Center) of
$\rho_{\odot}=0.3\,\mathrm{GeV/cm^{3}}=0.008\,M_{\odot}/\mathrm{pc}^{3}$.
Within this model, it is then straightforward to calculate the gravitational
potential $\Phi$ and the resulting acceleration
$\boldsymbol{\mathbf{a}}(\boldsymbol{\mathbf{r}})=-\boldsymbol{\mathbf{\nabla}}\Phi$
at any position $\boldsymbol{\mathbf{r}}$ in the MW. The apparent angular
acceleration is then given by

$\displaystyle\boldsymbol{\mathbf{\alpha}}(l,b,D_{\text{s}})=\frac{\delta
a}{D_{\text{s}}}\left[\hat{\boldsymbol{\mathbf{\delta
a}}}-\hat{\boldsymbol{\mathbf{D}}}_{s}(\hat{\boldsymbol{\mathbf{D}}}_{s}\cdot\hat{\boldsymbol{\mathbf{\delta
a}}})\right],$ (4.6)

in terms of the acceleration relative to that of the solar barycenter
$\boldsymbol{\mathbf{\delta
a}}\equiv\boldsymbol{\mathbf{a}}(\boldsymbol{\mathbf{r}})-\boldsymbol{\mathbf{a}}(\boldsymbol{\mathbf{r}}_{\odot})$,
the distance vector
$\boldsymbol{\mathbf{D}}_{s}\equiv\boldsymbol{\mathbf{r}}-\boldsymbol{\mathbf{r}}_{\odot}$
to the star, and their unit vectors $\hat{\boldsymbol{\mathbf{\delta a}}}$ and
$\hat{\boldsymbol{\mathbf{D}}}_{s}$, respectively. We plot the fiducial
angular acceleration of eq. 4.6 as a function of of galactic longitude $l$ and
latitude $b$ in figure 16, for three distances
$D_{\text{s}}=\\{200\,\mathrm{pc},1\,\mathrm{kpc},5\,\mathrm{kpc}\\}$. Typical
Galactic acceleration magnitudes are at the level of $0.3\times
10^{-3}\,\mathrm{\mu as}\,\mathrm{yr}^{-2}$.

An EPIC telescope can track the relative angular separation
$\boldsymbol{\mathbf{\theta}}_{ba}$ between two sources $a$ and $b$, and thus
the relative angular acceleration between two sources
$\boldsymbol{\mathbf{\alpha}}_{ba}\equiv\boldsymbol{\mathbf{\alpha}}_{b}-\boldsymbol{\mathbf{\alpha}}_{a}$
according to eq. 4.6. In practice, both of the components of
$\boldsymbol{\mathbf{\theta}}_{ba}$ would be fit to the model:

$\displaystyle\theta_{ba}^{j}(t)\simeq\varpi_{0}^{j}\cos\left(2\pi\frac{t}{\mathrm{yr}}+\gamma^{j}\right)+\theta_{0}^{j}+\mu_{0}^{j}(t-t_{0})+\alpha_{0}^{j}\frac{(t-t_{0})^{2}}{2}.$
(4.7)

where $j=1,2$ is a component superscript (e.g. $(\text{ra},\text{dec})$), and
$\varpi_{0}^{j}$, $\theta_{0}^{j}$, $\mu_{0}^{j}$, and $\alpha_{0}^{j}$ are
the fit parameters for the relative parallax, angular position, velocity, and
acceleration, respectively. The phases $\gamma^{j}$ are known based on the
position of the source. Assuming measurements are taken at a fixed repetition
frequency $f_{\text{rep}}$ over a survey time $\mathcal{T}$ (so a total number
of measurement epochs $f_{\text{rep}}\mathcal{T}$), the covariance matrix on
the estimated fit parameters for one component $p$ from a least-squares
analysis would then be:

$\displaystyle\mathrm{cov}\left[\widehat{\theta}_{i},\widehat{\theta}_{j}\right]=\frac{\sigma_{\delta\theta}^{2}}{f_{\mathrm{rep}}}\begin{pmatrix}\frac{2}{\mathcal{T}}&\frac{15}{\pi^{2}\mathcal{T}^{3}}&0&\frac{-360}{\pi^{2}\mathcal{T}^{5}}\\\
\frac{15}{\pi^{2}\mathcal{T}^{3}}&\frac{9}{4\mathcal{T}}&0&\frac{-30}{\mathcal{T}^{3}}\\\
0&0&\frac{12}{\mathcal{T}^{3}}&0\\\
\frac{-360}{\pi^{2}\mathcal{T}^{5}}&\frac{-30}{\mathcal{T}^{3}}&0&\frac{720}{\mathcal{T}^{5}}\end{pmatrix};\qquad\widehat{\theta}_{i}=\left\\{\widehat{\varpi}_{0}^{j},\widehat{\theta}_{0}^{j},\widehat{\mu}_{0}^{j},\widehat{\alpha}_{0}^{j}\right\\},$
(4.8)

where $\sigma_{\delta\theta}$ is the single-epoch light-centroiding precision
of eq. 2.50 as before.

The Galactic acceleration is prohibitively small to be measured by EPIC Phase
I, and difficult to detect on single star pairs even for EPIC Phase II, though
perhaps positive evidence for the effect can be obtained on aggregate. EPIC
Phase III, however, could map out the Galactic acceleration field on single
star pairs. Its relative light-centroiding precision for two Sun-like stars at
a distance of $200\,\mathrm{pc}$ would be $\sigma_{\delta\theta}\approx
0.11\,\mathrm{\mu as}$, and the resulting error on the relative acceleration
would be
$\sigma_{\delta\alpha}\simeq\sqrt{720/f_{\text{rep}}\mathcal{T}^{5}}\sigma_{\delta\theta}\approx
0.18\times 10^{-3}\,\mathrm{\mu as}\,\mathrm{yr}^{-2}$ for
$f_{\text{rep}}=1/\text{month}$ and $\mathcal{T}=30\,\mathrm{yr}$, smaller or
comparable to the expected accelerations of stars in the Milky Way shown in
figure 16. Brighter stars would yield higher accuracy and sensitivity to
accelerations at larger distances.

We leave a detailed projection on how a EPIC could be used to infer the
structure of the Milky Way (e.g. as done in ref. [15]) to future work. Such an
analysis would entail finding suitable bright source pairs, where the two
constituents of the pair would have significantly different line-of-sight
distances, otherwise their predicted relative accelerations would be too
similar, with little signal remaining in their relative angular accelerations.
An ideal configuration would consist of bright nearby stars close to distant
sources, such as red/blue giants or quasars.

###  4.4 Stellar microlensing

Gravitational microlensing is the phenomenon whereby a compact foreground
object, either luminous (e.g. a star) or dark (e.g. a black hole, neutron
star, or compact dark matter object), temporarily passes nearly in front of a
background source, bending its light and splitting its image into two (often
unresolved) images [165]. Such an event produces transient signals in both
photometric and astrometric observables [165, 166, 167, 168]. Due to
experimental realities, the photometric signal—a characteristic non-repeating
magnification event in the light curve of a background source—has historically
been the primary driver in the search for compact lenses [169, 170, 171], even
excluding primordial black holes as the totality of dark matter over a wide
mass range. Ground-based microlensing surveys such as Optical Gravitational
Lensing Experiment (OGLE) [172], Microlensing Observations in Astrophysics
(MOA) [173], and the Korea Microlensing Telescope Network (KMTNet) [174]
continue to detect thousands of photometric microlensing events every year in
dense stellar fields, and will soon be joined by a similar survey [146] with
the Nancy Grace Roman Space Telescope [175]. Astrometric microlensing has been
observationally challenging [176, 177, 178] but has recently become a useful
tool, yielding positive stellar microlensing detections [179, 180], confirming
the discovery of the first isolated black hole [3, 181, 182, 4], and spurring
forecasts for astrometry-only analyses of luminous lenses [183] and “blind”
searches for dark compact lenses [14, 184].

Since the optical depth for Galactic microlensing is tiny and building a vast
catalog with EPIC is infeasible, blind searches for dark compact lenses with
EPIC are not practical. However, EPIC’s excellent astrometric capabilities
enable precise parameter estimation of two classes of microlensing events.
Firstly, weak microlensing events between two luminous sources as predicted by
e.g. Gaia may pinpoint the lens mass along the lines of ref. [183] but with
higher precision. Secondly, for photometrically-flagged strong microlensing
events of dark compact lenses, the relative separation of the two source
images can be measured (as proposed in ref. [185]) by EPIC, again with great
potential for parameter estimation. In this section, we will explore EPIC’s
potential for these two use cases.

The characteristic angular scale of a microlensing event by a point-like lens
of mass $M_{\mathrm{L}}$ at a line-of-sight distance $D_{\mathrm{L}}$ of a
source at distance $D_{\mathrm{S}}$ is set by the (angular) Einstein radius:

$\displaystyle\theta_{\mathrm{E}}=\sqrt{\frac{4GM_{\mathrm{L}}}{c^{2}D_{\mathrm{L}}}\frac{D_{\mathrm{S}}-D_{\mathrm{L}}}{D_{\mathrm{S}}}}\approx
2.85\,\mathrm{mas}\sqrt{\frac{M_{\mathrm{L}}}{M_{\odot}}\frac{1\,\mathrm{kpc}}{D_{\mathrm{L}}}\frac{D_{\mathrm{S}}-D_{\mathrm{L}}}{D_{\mathrm{S}}}}.$
(4.9)

At an angular impact parameter $\boldsymbol{\mathbf{\beta}}$ from the lens to
the (unlensed) source location, the source image is split into two (denoted by
$+$ and $-$ subscripts), with astrometric deflection angles
$\boldsymbol{\mathbf{\Delta\theta}}_{\pm}$ and photometric magnifications
$\mu_{\pm}$:

$\displaystyle\boldsymbol{\mathbf{\Delta\theta}}_{\pm}$
$\displaystyle=\frac{\pm\sqrt{\beta^{2}+4\theta_{\mathrm{E}}^{2}}-\beta}{2}\hat{\boldsymbol{\mathbf{\beta}}};$
(4.10) $\displaystyle\mu_{\pm}$
$\displaystyle=\frac{\beta^{2}+2\theta_{\mathrm{E}}^{2}}{2\beta\sqrt{\beta^{2}+4\theta_{\mathrm{E}}^{2}}}\pm\frac{1}{2};$
(4.11) $\displaystyle\beta^{j}$
$\displaystyle\equiv\left(\hat{\boldsymbol{\mathbf{\theta}}}_{\mathrm{S}}-\hat{\boldsymbol{\mathbf{\theta}}}_{\mathrm{L}}\right)^{j}_{\mathrm{unlensed}}\equiv\varpi_{0}^{j}\cos\left(2\pi\frac{t}{\mathrm{yr}}+\gamma^{j}\right)+\theta_{0}^{j}+\mu_{0}^{j}(t-t_{0}).$
(4.12)

On the RHS of the last line, we have parametrized the unlensed trajectory by 5
parameters: the relative parallax $\varpi_{0}$ (whose direction and phases
$\gamma_{j}$ are known), and the relative separation
$\boldsymbol{\mathbf{\theta}}_{0}$ and proper motion
$\boldsymbol{\mathbf{\mu}}_{0}$ (two components each, not to be confused with
the magnifications $\mu_{\pm}$) at some reference time $t_{0}$ (taken to be
the midpoint of the observation period).

In this notation, the two aforementioned signatures for EPIC are as follows.
Analogous to the Gaia forecasts of ref. [183], EPIC can measure the relative
separation of a luminous lens and (the primary “$+$” image of) a source

$\displaystyle\boldsymbol{\mathbf{\theta}}_{\mathrm{SL}}\equiv\hat{\boldsymbol{\mathbf{\theta}}}_{\mathrm{S}+}-\hat{\boldsymbol{\mathbf{\theta}}}_{\mathrm{L}}=\boldsymbol{\mathbf{\beta}}+\boldsymbol{\mathbf{\Delta\theta}}_{+}\quad\text{(weak
lensing, luminous lens)}.$ (4.13)

Likewise, EPIC can also be employed to execute the proposal from ref. [185] to
measure the relative separation $\boldsymbol{\mathbf{\theta}}_{+-}$ between
the two images of the source in a microlensing event by a non-luminous lens:

$\displaystyle\boldsymbol{\mathbf{\theta}}_{+-}\equiv\hat{\boldsymbol{\mathbf{\theta}}}_{\mathrm{S}+}-\hat{\boldsymbol{\mathbf{\theta}}}_{\mathrm{S}-}=\sqrt{\beta^{2}+4\theta_{\mathrm{E}}^{2}}\,\hat{\boldsymbol{\mathbf{\beta}}}\quad\text{(strong
lensing, dark lens)}.$ (4.14)

The observables in eqs. 4.13 and 4.14 can both be fit by 6 parameters: the 5
parameters in eq. 4.12, and the lens mass $M_{\mathrm{L}}$, using the well-
known fact that one can write
$\theta_{\mathrm{E}}^{2}=4GM_{\mathrm{L}}\varpi_{0}/\mathrm{AU}$ in terms of
the lens mass $M_{\mathrm{L}}$ and the relative parallax $\varpi_{0}$. In both
cases, the lens mass can be extracted directly, unlike in photometry-only
microlensing measurements with only the total magnification $\mu_{+}+\mu_{-}$
as an observable.

To assess EPIC’s astrometric microlensing capabilities, we generate two mock
data sets. The first is for the weak microlensing event
4410342766254199168-35742 between two luminous stars predicted from Gaia data
in ref. [186]. We assume “truth values” for the system of
$\\{-46.1\,\text{mas}$, $45.8\,\text{mas}$, $108.04\,\text{mas/yr}$,
$-118.78\,\text{mas/yr}$, $11.4\,\text{mas}$, $1.3\,M_{\odot}\\}$ for
$\\{\theta_{0}^{\alpha}$, $\theta_{0}^{\delta}$, $\mu_{0}^{\alpha*}$,
$\mu_{0}^{\delta}$, $\varpi_{0}$, $M_{\text{L}}\\}$, the two components each
(right ascension $\alpha$ and declination $\delta$, with
$\mu_{0}^{\alpha*}\equiv\mu_{0}^{\alpha}\cos\delta$) of relative separation
$\boldsymbol{\mathbf{\theta}}_{0}$ and proper motion
$\boldsymbol{\mathbf{\mu}}_{0}$, and for the relative parallax $\varpi_{0}$
and lens mass $M_{\mathrm{L}}$. We set the parallax phases $\gamma_{j}$ to
zero and $3\pi/2$ for the right ascension and declination components,
respectively. We generate a mock data set according to eq. 4.13 comprised of
six relative angular separation measurements, spanning a period of a year
around the closest approach. We assume an astrometric light-centroiding
precision of $\sigma_{\delta\theta}=\,11\mathrm{\mu as}$ per epoch, which is
the expectation per observing night given the quoted apparent magnitude of the
two stars in ref. [186] and Phase II instrumental parameters of an EPIC
interferometer. We perform an MCMC parameter estimation of
$\\{\theta_{0}^{\alpha}$, $\theta_{0}^{\delta}$, $\mu_{0}^{\alpha*}$,
${\mu_{0}^{\delta}}$, $\varpi_{0}$, $M_{\text{L}}\\}$, whose posteriors we
show in figure 17.

The second mock data set is for the black-hole microlensing event
MOA-2011-BLG-191/OGLE-2011-BLG-0462 from refs. [3, 181, 182, 4]. We create a
(counterfactual) simulation of the dark lens’ mass estimation had EPIC taken
data around the time that OGLE detected the photometric amplification. For our
mock data generation, we assume a relative separation of
$\\{\theta_{0}^{\alpha},\theta_{0}^{\delta}\\}=\\{-0.23,-0.17\\}\,\mathrm{mas}$
at the time of closest (unlensed) approach (taken to be the reference epoch
$t_{0}$), commensurate with the reported closest approach of
$\mathrm{min}\,\beta\approx 0.05\,\theta_{\mathrm{E}}$ and Einstein angle
$\theta_{\mathrm{E}}\approx 5.81\,\mathrm{mas}$ from ref. [182]. We take the
mean values of their inferred velocities as our “truth parameters”, namely
$\\{\mu^{\alpha*}_{\text{L}},\mu^{\delta}_{\text{L}}\\}=\\{4.36,3.06\\}\,\mathrm{mas/yr}$
and
$\\{\mu^{\alpha*}_{\text{S}},\mu^{\delta}_{\text{S}}\\}=\\{2.263,-3.597\\}\,\mathrm{mas/yr}$,
with the relative proper motion
$\boldsymbol{\mathbf{\mu}}_{0}\equiv\boldsymbol{\mathbf{\mu}}_{\mathrm{S}}-\boldsymbol{\mathbf{\mu}}_{\mathrm{L}}$.
Likewise, we take the mass of the dark lens to be
$M_{\mathrm{L}}=7.1\,M_{\odot}$ and a relative parallax magnitude of
$\varpi_{0}=0.6\,\mathrm{mas}$. We also do not study the effects of the
unequal brightness of the two images, nor their lensing-induced elongation
(primarily along the axis perpendicular to
$\hat{\boldsymbol{\mathbf{\theta}}}_{+-}$), which could slightly affect the
form factor.

Given the faintness of this particular source star, whose I-band magnitude
varies in the range $13\lesssim m_{\mathrm{I}}\lesssim 16$ during the
microlensing event, Phase II of EPIC would yield a low SNR, so Phase III
parameters might be needed for high significance. On the other hand, the
source is far enough that larger baselines can be employed: a Sun at the
$1.5\,\mathrm{kpc}$ distance to the lens, would have an optimal baseline of
$4\,\mathrm{km}$ baseline, yielding an angular resolution of
$\sigma_{\theta_{\mathrm{res}}}\approx 6\,\mathrm{\mu as}$ in the I band, and
even better light centroiding. Assuming a conservative light-centroiding
precision of $\sigma_{\delta\theta}\simeq\sigma_{\theta_{\mathrm{res}}}$, we
generate a mock data set which comprises six relative-angle measurements
between the two images of the source, according to eq. 4.14, interspersed
within two years around the closest approach. We perform an MCMC parameter
estimation of
$\\{\theta_{0}^{\alpha},\theta_{0}^{\delta},\mu_{0}^{\alpha*},\mu_{0}^{\delta},\varpi_{0},M_{\text{L}}\\}$,
whose posteriors we show in figure 18.

The capability of Gaia in determining stellar masses from star-star
microlensing is projected to be at best a few percent: according to ref.
[183], the mission can achieve stellar mass determination of about 5 to 10%
for about seven events and below 15% for thirteen events total. As seen from
the marginalized posterior for the lens mass in figure 17, even with the
conservative assumptions of a Phase II interferometer and just six observing
epochs, EPIC can greatly surpass Gaia’s capability per microlensing event by
at least an order of magnitude, thereby achieving a 0.46% determination of the
lens mass.

The microlensing event MOA-2011-BLG-191/OGLE-2011-BLG-0462 from an isolated
black hole yielded a determination of the dark lens’ mass to about 18% using
HST astrometry [3, 181, 182, 4]. The observable in this case was the
luminosity-weighted light centroid separation of _both_ images of the source
star, since _HST_ cannot resolve them from each other. In contrast, an EPIC
interferometer, much like the proposal of ref. [185] to use the VLT, can
directly resolve the angular separation of the two images, thereby allowing
for an order of magnitude better measurement of the dark lens’ mass, down to a
precision level of 1.7% with Phase II specifications, as seen from figure 18.
Phase III could improve the significance and light centroiding by more than an
order of magnitude, far surpassing the capabilities of optical
interferometers, whose baselines can be at most a few hundred meters and
light-centroiding precision of about $10\,\mathrm{\mu as}$ [185]. It would be
interesting to pursue a more detailed analysis of dark lens microlensing
events, including the modification of the stellar sizes due to elongation,
which could give a form factor suppression akin to that of eq. 2.34 if
resolved, but could also provide information about the lens structure.

Figure 17: MCMC posteriors of the full relative astrometric solution
$\\{\theta_{0}^{\alpha},\,\theta_{0}^{\delta},\,\mu_{0}^{\alpha*},\,\mu^{\delta}_{0},\,\varpi_{0},\,M_{\text{L}}\\}$
from a mock data set of six relative separation measurements of the
microlensing event 4410342766254199168-35742 predicted in ref. [186]. The
measurements span a period of one year around the time of closest approach of
the two stars, each having a precision $\sigma_{\theta}=11\,\mathrm{\mu as}$
Remarkably, EPIC Phase II can determine the mass of the stellar lens at the
$0.46\%$ level, even with just enough measurements to break the parameter
degeneracy. Figure 18: MCMC posteriors of the full relative astrometric
solution
$\\{\theta_{0}^{\alpha},\,\theta_{0}^{\delta},\,\mu_{0}^{\alpha*},\,\mu_{0}^{\delta},\,\varpi_{0},\,M_{\text{L}}\\}$
from a mock data set of six relative separation measurements of the two images
of MOA-2011-BLG-191/OGLE-2011-BLG-0462 [3, 181, 182, 4]. The measurements span
a period of two years around the time of closest approach of the star to the
dark lens, each having a precision $\sigma_{\theta}=6\,\mathrm{\mu as}$
assuming Phase II instrumental parameters. The achieved determination of the
mass of the dark lens—in this case an isolated black hole—is at the $1.7\%$
level.

###  4.5 Strongly lensed quasars

Intensity interferometry, including the extended-path variation detailed in
this work, is only feasible on compact sources with a high surface brightness.
The most obvious sources _at cosmological distances_ are quasars, extremely
luminous emission regions powered by the accretion of matter onto a
supermassive black hole. (Other potential targets at large distances include
dense stellar clusters and supernovae, but intensity interferometry
applications for these sources are left to future work.) Their overall
apparent brightness is not too small for advanced Phases of EPIC. The
morphology of quasar emission regions is not well understood, and only an
intensity interferometer may have the requisite effective angular resolution
to directly resolve the nano-arcsecond structure of quasar central engines.
Complementary indirect methods have also been proposed [187].

About 1 in 500 quasars are known to be strongly lensed into two or more images
[188]. These images are often magnified by about an order of magnitude in
apparent brightness, and are typically separated by about an arcsecond.
(Smaller image separations are also possible but are not easily resolved, so
some “single quasars” may actually consist of two or more duplicate images.)
The first such system was $\mathrm{Q0957+561}$ [189], while the most striking
and best-studied system is $\mathrm{Q2237+0305}$ [190, 191]. A few hundred
multiply-imaged quasars are known at present [192, 193], with an explosion of
this sample (approaching $10^{4}$) expected from upcoming surveys [194], most
notably LSST but also the Dark Energy Survey, DESI, and Euclid.

EPIC is ideally suited to perform relative astrometry on any pair of strongly-
lensed quasar images. The relative locations of these images are not so
interesting per se, but with its exquisite light-centroiding precision, EPIC
can track _changes_ in their relative separation. Such changes can occur due
to the (relative) peculiar velocities of the quasar-lens system, which lead to
differential proper motion between the various quasar images because of the
relative distortions caused by the lensing potential. The apparent image
proper motions also sweep through the fine-grained structure of the
gravitational potential, most notably the micro-caustic network of stars in
the lens galaxy. Near any micro-caustic crossing, one expects (correlated)
photometric [191, 195, 196, 197, 198] and astrometric [188, 199, 200]
variability. A companion paper [128] to this work fleshes out an entirely new
signature—stochastic astrometric weak lensing fluctuations from low-mass dark
matter halos.

In this work, we will content ourselves with making rough projections on the
relative astrometric light-centroiding precision between two images, and how
this information can be used to measure the effective proper motion and macro-
lensing Jacobians. These measurements would allow, for the first time, highly
significant determinations of transverse proper motions of objects at
cosmological distances and strong-lensing geometries for multiply-imaged
quasars in their respective rest frames. The lensing map from an image-plane
location $\boldsymbol{\mathbf{\theta}}$ to a source-plane location
$\boldsymbol{\mathbf{\beta}}$ is given by the lensing map
$\boldsymbol{\mathbf{\varphi}}$,

$\displaystyle\boldsymbol{\mathbf{\beta}}=\boldsymbol{\mathbf{\varphi}}(\boldsymbol{\mathbf{\theta}})=\boldsymbol{\mathbf{\theta}}-\boldsymbol{\mathbf{\alpha}}(\boldsymbol{\mathbf{\theta}}),$
(4.15)

with $\boldsymbol{\mathbf{\alpha}}(\boldsymbol{\mathbf{\theta}})$ the well-
known (reduced) deflection angle, which can be written as:

$\displaystyle\boldsymbol{\mathbf{\alpha}}(\boldsymbol{\mathbf{\theta}})=2\int_{0}^{D_{\mathrm{q}}}\mathrm{d}D^{\prime}\frac{D_{\mathrm{q}}-D^{\prime}}{D_{\mathrm{q}}}\boldsymbol{\mathbf{\nabla}}_{\perp}\Phi(D^{\prime}\boldsymbol{\mathbf{\theta}},D^{\prime}).$
(4.16)

In the above, $D_{\mathrm{q}}$ is the angular diameter distance to the quasar,
and $\Phi$ is the 3D gravitational potential, which in almost all cases of
relevance is dominated by the lens galaxy at some distance $D_{\mathrm{L}}$.
The macro-lensing Jacobian is the determinant of the lensing map, and is given
by:

$\displaystyle A_{ij}\equiv\frac{\partial\beta_{i}}{\partial\theta_{j}}.$
(4.17)

Let $B$ be the inverse of the $2\times 2$ matrix $A$, i.e.
$B_{ij}=(A^{-1})_{ij}$, which is a function of $\boldsymbol{\mathbf{\theta}}$.
The total image magnification is given by $\det B$. The inverse map
$\boldsymbol{\mathbf{\varphi}}^{-1}:\boldsymbol{\mathbf{\beta}}\to\boldsymbol{\mathbf{\theta}}$
is, in general, multi-valued (leading to multiple images). Let us denote the
multiple quasar image locations as
$\boldsymbol{\mathbf{\theta}}^{\mathrm{I}}=\boldsymbol{\mathbf{\varphi}}^{-1}(\boldsymbol{\mathbf{\beta}}_{\mathrm{q}})|_{\mathrm{I}}$
with $\mathrm{I}=\mathrm{A,B,\ldots}$ a discrete label for the images, and
$\boldsymbol{\mathbf{\beta}}_{\mathrm{q}}$ the quasar’s location on the source
plane.

The effective apparent angular motion between the quasar and the lens galaxy
(in the source plane) is

$\displaystyle\boldsymbol{\mathbf{\mu}}\equiv\frac{\mathrm{d}\boldsymbol{\mathbf{\beta}}}{\mathrm{d}t}=\frac{\boldsymbol{\mathbf{v}}_{\mathrm{q}}}{1+z_{\mathrm{q}}}\frac{1}{D_{\mathrm{q}}}-\frac{\boldsymbol{\mathbf{v}}_{\mathrm{L}}}{1+z_{\mathrm{L}}}\frac{1}{D_{\mathrm{L}}}+\frac{\boldsymbol{\mathbf{v}}_{\mathrm{o}}}{1+z_{\mathrm{L}}}\frac{D_{\mathrm{q}}-D_{\mathrm{L}}}{D_{\mathrm{L}}D_{\mathrm{q}}}$
(4.18)

where $\boldsymbol{\mathbf{v}}_{\mathrm{q}}$,
$\boldsymbol{\mathbf{v}}_{\mathrm{L}}$, and
$\boldsymbol{\mathbf{v}}_{\mathrm{o}}$ are the proper velocities of the
quasar, lens galaxy, and observer, respectively [201]. With these definitions,
the trajectory for each of the images is:

$\displaystyle\theta^{\mathrm{I}}_{i}(t)=\theta^{\mathrm{I}}_{i}(t_{0})+\sum_{j}B^{\mathrm{I}}_{ij}\left[\mu_{j}(t-t_{0})+\varpi_{j}\cos\left(2\pi\frac{t}{\mathrm{yr}}+\beta_{j}\right)\right]$
(4.19)

with $B^{\mathrm{I}}_{ij}\equiv
B_{ij}[\boldsymbol{\mathbf{\theta}}^{\mathrm{I}}(t_{0})]$. For example, for
Huchra’s lens at $z_{\mathrm{L}}=0.0394$ which quadruply lenses the quasar
$\mathrm{Q2237{+}0305}$, we have

$\displaystyle\mu\simeq\frac{v_{\mathrm{L}}}{(1+z_{\mathrm{L}})D_{\mathrm{L}}}\approx
0.252\,\mathrm{\mu
as/yr}\left(\frac{v_{\mathrm{L}}}{200\,\mathrm{km/s}}\right);\qquad\varpi\simeq\frac{\mathrm{AU}}{(1+z_{\mathrm{L}})D_{\mathrm{L}}}\approx
0.0060\,\mathrm{\mu as}.$ (4.20)

These astrometric deviations are then further magnified by the strong lens,
parametrized by the (inverse) Jacobian $B^{\mathrm{I}}_{ij}$, which can
further amplify astrometric perturbations by up to a factor of 10 (or even
higher [202]) along some direction(s), depending on the image and strong lens
under consideration.

Figure 19: Light-centroiding precision $\sigma_{\delta\theta}$ from eq. 2.51
(at the optimal baseline distance $d^{\text{opt}}$ of eq. 2.52) between two
identical quasar images with apparent angular optical radius $\theta_{500}$.
For EPIC Phase III, we project that
$\sigma_{\delta\theta,\mathrm{q}}^{\mathrm{opt}}\approx 4.64\times
10^{-13}\,\mathrm{rad}\,(10^{-11}\,\mathrm{rad}/\theta_{500})$ in our simple
quasar emission model.

An error estimate analogous to the one in eq. 4.8 for this simpler model
without acceleration gives uncertainties on the estimators of the parallax,
relative separation, and proper motion (in each direction $j=1,2$) of:

$\displaystyle\begin{pmatrix}\sigma_{\widehat{\varpi}_{0}^{j}}\\\
\sigma_{\widehat{\theta}_{0}^{j}}\\\
\sigma_{\widehat{\mu}_{0}^{j}}\end{pmatrix}=\frac{\sigma_{\delta\theta}}{\sqrt{f_{\mathrm{rep}}\mathcal{T}}}\begin{pmatrix}\sqrt{2}\\\
1\\\
\frac{2\sqrt{3}}{\mathcal{T}}\end{pmatrix}\approx\begin{pmatrix}0.168\,\mathrm{\mu
as}\\\ 0.119\,\mathrm{\mu as}\\\ 0.0138\,\mathrm{\mu
as}/\mathrm{yr}\end{pmatrix}\left(\frac{\sigma_{\delta\theta}}{10^{-11}\,\mathrm{rad}}\right)\left(\frac{10/\mathrm{yr}}{f_{\mathrm{rep}}}\right)^{1/2}$
(4.21)

with the latter estimate at $\mathcal{T}=30\,\mathrm{yr}$. The quoted
uncertainties are for the apparent parallax and proper motion in the image
plane; the magnification by the $B^{\mathrm{I}}_{ij}$ factors for each image
in the pair will further reduce the uncertainties for the quantities in the
source plane. We plot the expected light-centroiding precision
$\sigma_{\delta\theta}$ from eq. 2.51 (at the optimal baseline distance
$d^{\text{opt}}$ of eq. 2.52) between two identical quasars in figure 19 for
the three Phases of EPIC. Comparison of eqs. 4.20 and 4.21 shows that the
relative motions between any image pair should be measurable at high
fractional precision over a long time period with EPIC Phase II and III, while
the parallax is extremely challenging to measure, even for this especially
close lens.

In the simple model of eq. 4.19, we have ignored the effects of fine-grained
structures in the lens galaxy, both from the stellar population and from dark
matter substructures. Astrometric deflections from stellar substructure have
been proposed as potentially observed signatures in ref. [188, 200] and from
medium-scale dark matter subhalos in ref. [199]; these are all interesting
effects that can be explored with an intensity interferometer. A companion
paper [128] will explore the astrometric noise caused by dark matter
substructures at extremely small scales (substellar masses).

##  5 Conclusions

We have proposed the _Extended-Path Intensity Correlator_ (EPIC) [1], an
optical modification of intensity interferometers that introduces an
adjustable, differential geometric extension in the light paths from widely-
separated sources. EPIC expands the field of view of traditional intensity
interferometry—by up to six orders of magnitude beyond the angular
resolution—while preserving the light-centroiding precision. In detail, the
EPIC setup consists of a large primary mirror and a collimating secondary
mirror that guides the light into the path-extension system (figure 5). The
latter has the form of a Mach-Zehnder interferometer of unequal path lengths,
and serves a dual role: it introduces the geometric path extension and
corrects for the relative tilt of the sources’ wavefronts. The wavefront tilt
correction avoids the fringe suppression present when a constituent telescope
resolves the source separation. The path extension allows for angular
distances between sources at arcsecond separations to be determined with
microarcsecond precision.

Our proposal arrives at a time when technological advances can be employed
together to increase the signal-to-noise ratio and extend the intensity
interferometry technique to a broader range of targets. We develop a
multichannel observation scheme over a broad frequency range, combined with
large telescopes, fast photodetectors, and synchronization over large baseline
distances. We outline the current capabilities in spectroscopy and picosecond-
level timing, and identify the specific advances needed from the recently-
demonstrated kilopixel single-photon detector arrays for the EPIC proposal.
Modern development of atomic clocks and infrastructure for large baseline
radio interferometry make possible picosecond-level synchronization and
millimeter-level baseline characterization over the hours-long observation
times and kilometer-plus distances needed for intensity interferometry.

We find that EPIC will be able to target stars as faint as magnitude 15, as
well as bright quasars, with telescope apertures no larger than
$10\,\mathrm{m}$ in diameter. The maximum angular source separation for EPIC
is limited by the atmosphere’s isoplanatic angle—a few arcseconds in the
optical—within which atmospheric fluctuations are correlated.

This work improves on differential astrometry strategies with intensity
interferometry in several ways. The feasibility of measuring the relative
separation of sources with traditional interferometry was first demonstrated
on the close spectroscopic binary $\alpha$ Vir at the original NSII [42],
which determined all six Keplerian orbital parameters. Recently, ref. [60] has
proposed to target widely separated sources (by more than the isoplanatic
angle) with a _fixed_ path extension using a single beam splitter and a GHz
photon bandwidth per spectral channel. Such a setup imposes exacting (and
potentially infeasible) requirements on the adaptive optics and spectroscopic
systems. Using the _variable_ path extension proposed here, we can track the
separation of sources at or near the main path-extended fringe (center of
figure 4(b)) and employ the differential wavefront tilt correction to avoid
fringe suppression for sources separated by more than a few milliarcseconds.
Our approach yields a per-epoch light-centroiding precision that scales as
$t_{\mathrm{obs}}^{-1/2}$ (eq. 2.51), compared to the fringe-scanning method
of ref. [60] at extremely high fringe orders, which could lead to fringe
confusion and scales as $t_{\mathrm{obs}}^{-3/2}$.

In addition, while multi-band observation schemes have been proposed in recent
years [50, 60, 61], this work is the first to study the optical tolerances
induced by such a scheme. The few picosecond time-resolution assumed in this
work allows for errors in optical paths as long as a millimeter. The desired
spectral resolution and the small pixel sizes of ultra-fast photon counters
are the limiting factors, and can tolerate primary-mirror imperfections as
large as tens of wavelengths, in contrast to imaging telescopes and amplitude
interferometers that require sub-wavelength polishing [115].

We also provide the first quantitative study of a fringe contrast suppression
for very fast photodetectors and large pixels—unique to ultra-fast and multi-
channel intensity interferometry and an entirely new effect unveiled in this
work. This perhaps surprising effect originates for setups with photodetector
arrays and spectroscopic splitting: atmospheric and instrumental aberrations
require the angle subtended by each pixel to be larger than the diffraction
limit of each telescope, resulting in uncorrelated fluctuations within a
pixel. We identify several possible solutions to this problem, and focus on
the conceptually simplest: splitting the beam exiting the EPIC optical delay
system into several sub-beams, each with its own dispersion element and
photodetector array.

Finally, we calculate effects from atmospheric aberrations and their
limitation on the maximum source separation of order the isoplanatic angle for
a ground-based intensity interferometer; correcting for atmospheric turbulence
over angular separations significantly larger than the isoplanatic patch is
beyond the capabilities of present adaptive optics systems [115]. In the
future, should a lunar observatory be established, an EPIC observatory would
be an clear candidate that would benefit from the absence of an atmosphere;
study of a large opening angle design without atmospheric limitations is left
to future work.

We present several scientific applications that would experience major leaps
forward by leveraging the differential astrometric capabilities of EPIC. The
target sources of interest in the Milky Way are binary systems, accidental
double stars, and strong microlensing events. EPIC could detect exoplanets
astrometrically through the wobble they impart on their host stars with a
reach orders of magnitude better than Gaia’s projected DR5 [1]. An EPIC
interferometer targeting a binary system could determine all Keplerian angles
to a sub-degree level and, combined with radial velocity measurements,
determine the masses of the constituent stars and the line-of-sight distance
to the per-mil level or better. This level of precision would open up
possibilities of measurements of stellar mass loss and cosmic distance ladder
calibration. Observations of accidental double stars can be used to detect
astrometric microlensing and determine stellar masses to the sub-percent level
and, with later Phases of the setup, map the potential of the Milky Way by
directly measuring the induced stellar accelerations at the
$\mathrm{nas/yr^{2}}$ level, significantly surpassing Gaia sensitivity [15].
Finally, EPIC could target _extragalactic_ sources such as quasars (and
possibly supernovae). We study the specific case of measuring the astrometric
evolution of the images in a multiply-imaged quasar to characterize the
strong-lensing system and determine peculiar velocities at cosmological
distances.

Directions for future work are numerous. At the level of design, identifying a
simpler solution to the aforementioned spectroscopic suppression is important
to reducing both complexity and cost, potentially allowing the construction of
several baselines and many telescopes per observation site. Potential sources
of observational systematic uncertainties, such as contamination from ambient
light or other sources in a crowded field, should be characterized in more
detail. On the experimental side, it is crucial to demonstrate kilo-pixel
arrays of single-photon detection devices tailored to the requirements of
EPIC, including broadband sensitivity to a wide range of wavelengths.
Naturally, a laboratory demonstration of the EPIC’s path-extension system
would be the first step to large-scale implementation.

Further study of the full scientific potential of the applications is also of
interest. While we use real example systems to illustrate the power of our
technique, we have not performed an exhaustive survey of potential sources.
For instance, thousands of microlensing events are predicted to occur in the
next few decades based on the Gaia catalog alone [186], hundreds of which
could be targeted by EPIC. Gaia data can also be used to construct a larger
double-star catalog to map out the Milky Way’s potential via the relative
transverse angular accelerations it imparts on source pairs. Finding
sufficiently bright and distant binary systems of Cepheids and other giant
stars would aid in refining the first rung of the cosmic distance ladder.

We have also foretold the possibility to probe sub-solar-mass dark matter
substructure using astrometry on the multiple images of strongly lensed
quasars, which will be studied in detail in a future paper [128]. In addition,
we expect EPIC may test general relativity in the vicinity of Sagittarius A*
[203, 204, 205], with potentially sufficient sensitivity to measure the spin
of the central black hole and surrounding mass distribution [206]. Finally, it
would be interesting to generalize our analysis of multi-channel intensity
interferometry (with and without path extension) to other extended, high
surface brightness sources, such as the various types of novae and the cores
of globular clusters and galaxies, and investigate corresponding applications.

Extended-Path Intensity Correlation will enable new sets of astronomical
measurements at unprecedented precision. It has the potential to transform our
knowledge about stellar systems and the Milky Way, and may ultimately provide
a new window into the dark universe and physics Beyond the Standard Model.

## Acknowledgments

We thank Gordon Baym, Megan Bedell, Karl Berggren, Michael Blanton, Matteo
Cantiello, Cyril Creque-Sarbinowski, Liang Dai, Neal Dalal, Julianne
Dalcanton, Peter Graham, David Hogg, Miguel Morales, Oren Slone, and David
Spergel for valuable conversations and input, and Calvin Chen, Jacob Crawford,
Cyril Creque-Sarbinowski, David Dunsky, Marius Kongsore, and Jessie Yang for
comments on the manuscript. KVT thanks Jason Aufdenberg, Matthew Brown, James
Buckley, Dainis Dravins, David Kieda, Michael Lisa, Nolan Matthews, Andrei
Nomerotski, Ue-Li Pen, Naomi Vogel, Shiang-Yu Wang, and Luca Zampieri for
fruitful conversations, and Sebastian Karl for pointing out ref. [63], during
the 2023 Workshop on Stellar Intensity Interferometry at The Ohio State
University. KVT is supported by the National Science Foundation under Grant
PHY-2210551. NW is supported by NSF under award PHY-2210498, by the BSF under
grant 2018140, and by the Simons Foundation. MB is supported by the DOE Office
of Science under Award Number DE- SC0022348, University of Washington Royal
Research Fund and through the Department of Physics and College of Arts and
Science at the University of Washington.

The Center for Computational Astrophysics at the Flatiron Institute is
supported by the Simons Foundation. Research at Perimeter Institute is
supported in part by the Government of Canada through the Department of
Innovation, Science and Economic Development and by the Province of Ontario
through the Ministry of Colleges and Universities.

Part of this work was performed at the Aspen Center for Physics, which is
supported by National Science Foundation grant PHY-1607611 and PHY-2210452.
The participation of MB at the Aspen Center for Physics was supported by the
Simons Foundation. MB, MG, and KVT thank the Institute for Nuclear Theory at
the University of Washington for its kind hospitality and stimulating research
environment. The INT is supported in part by the U.S. Department of Energy
grant No. DE-FG02-00ER41132.

This work has made use of data from the European Space Agency (ESA) mission
Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing
and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. This research has made use
of the NASA Exoplanet Archive, which is operated by the California Institute
of Technology, under contract with the National Aeronautics and Space
Administration under the Exoplanet Exploration Program.

##  Appendix A Notation

In tables 2–7, we relist all variables used in this work in the first column,
along with their definitions in the second column.

Variable | Definition
---|---
$E_{p}$ |  Total electric field at telescope $p$
$E_{i\alpha}$ |  Electric field amplitude of $i$ and of energy $k_{\alpha}$
$I$, $\widetilde{I}_{j}$ |  Intensity, fractional intensity of emitter $j$
$k=2\pi/\lambda$, $\overline{k}$ | Wavenumber (in terms of wavelength), mean wavenumber
$\hat{\boldsymbol{\mathbf{\theta}}}_{j}$ |  Unit vector to $j$
$\boldsymbol{\mathbf{\theta}}_{ij}=\hat{\boldsymbol{\mathbf{\theta}}}_{i}-\hat{\boldsymbol{\mathbf{\theta}}}_{j}$ |  Relative angular distance of $i$ and $j$.
$r_{ip}$ |  Distance of $i$ from telescope $p$
$\phi^{\text{em}}_{i\alpha}$ |  Intrinsic phase of emitter $i$ and energy $k_{\alpha}$
$\widetilde{\phi}_{i\alpha}^{(p)}$ |  Aberration phase at telescope $p$
${\text{d}I}/({\text{d}k\text{d}\Omega})$ | Intensity (power per area) per unit wavenumber and unit solid angle
Table 2: Electric field variables. Variable | Definition
---|---
$\mathbb{d}$ | Baseline vector
$\mathbb{b}$ | Aperture vector
$A_{i}$ |  Area of aperture of telescope $i$
$D$ | Diameter of primary mirror
$f$, $f^{\prime}$, $f_{\text{c}}$ | Focal lengths of primary, secondary, and camera
$n_{p}^{\text{arr}}$ |  Number of telescopes at array site $p$
$\alpha_{nlm}$ | Optical aberration amplitude
$\ell_{p}$ |  Geometric path extension at telescope $p$
$\gamma$ | Angle of first beam-splitter
$\gamma^{\prime},\delta^{\prime}$ | Angles of wavefront correction mirror
$\epsilon$ | Transmission coefficient of beam splitter
$N$ | Number of illuminated lines on diffraction grating
$d_{\text{g}}$ | Periodicity of grating
$s$ | Slit size of grating
$\theta_{\text{F}}$ | Angle from camera to the pixel plane
$\mathcal{R}$ | Spectroscopic resolution
$\sigma_{k}$ | Range of wavenumbers on incident on pixel
$\tau$ | Timing offset
$\sigma_{t}/\sqrt{2}$ | Timing precision per photodetector
$t_{\text{obs}}$ | Observation time per epoch
$t_{\text{jitter}}$ | Time jitter of photodetector (FWHM)
$t_{\text{reset}}$ | Photodetector reset time
DCR | Dark count rate of photodetector
$N_{p}$ |  Total number of photons incident at $p$ after time $t_{\text{obs}}$
$\eta_{p}$ |  Photodetection efficiency of $p$
$C(\mathbb{d},\tau)$ | Intensity correlation
$\hat{X}$ |  Estimator of $X$
$\left<X\right>$ |  Expectation value of $X$ (averaged over mode phases)
$\overline{X(t)}$ |  Time average of $X(t)$ over observation time
$\sigma_{\hat{C}}$ |  Standard deviation of $\hat{C}$
$\sigma_{\theta_{\text{res}}}$ | Angular resolution of the interferometer
$\sigma_{\Delta\theta}$ | Angular dynamic range of the interferometer
$\sigma_{\hat{\boldsymbol{\theta}}}$ | Global astrometric precision
Table 3: Instrumental variables. Variable | Definition
---|---
$n[\mathbb{r}]$ |  Atmospheric index at 3D position $\mathbb{r}$
$D_{n}(r)$ |  Structure function, $r$-dependent
$h$ | Altitude in the atmosphere
$C_{n}(h)$ |  Structure function, dependence on altitude $h$
$l_{0}$ | Inner turbulence scale
$L_{0}$ | Outer turbulence scale
$B_{n}(r)$ | Correlation function of atmospheric index of refraction
$\theta_{0}(\overline{k},\gamma)$ |  Isoplanatic patch angle at an angle $\gamma$ away from zenith
$\sigma_{\phi}^{\text{atm}}$ | Standard deviation of differential phase fluctuation due to the atmosphere
$\theta_{\text{seeing}}$ | Seeing angle
$r_{0}$ | Fried parameter
Table 4: Atmospheric fluctuation variables from section 2.5, appendix B.2, and appendix C.3. Variable | Definition
---|---
$R_{\text{s}}$ | Stellar radius
$D_{\text{s}}$ | Distance to star/binary
$\theta_{\text{s}}$ | Angular radius of star
$M_{j}$ |  Mass of star $j$
$T_{j}$ |  Temperature of $j$
$\mathcal{F}_{j}$ |  Finite-size form factor of $j$
$T_{500}=2\pi\hbar c/500\,\mathrm{nm}$ |  Black body temperature corresponding to $500\,\mathrm{nm}$
$R_{\text{500}}$ |  Physical radius in quasar’s local frame with temperature $T_{\text{500}}$
$z_{\mathrm{q}}$ | Quasar redshift
$R_{\text{isco}}$ | Innermost stable circular orbit radius of quasar
$T_{\mathrm{q},0}$ | Quasar reference temperature
Table 5: Single source parameters. Variable | Definition
---|---
$f_{\text{rep}}$ | Observational repetition frequency
$G$ | Gravitational constant
$M$ | Total binary mass
$a$ | Semi-major axis
$e$ | Eccentricity
$i$ | Inclination
$\Omega$ | Longitude of ascending node
$\omega$ | Argument of periapsis
$\nu$ | True anomaly
$t_{0}$ | Reference epoch
$E$ | Eccentric anomaly
$\mathcal{M}$ | Mean anomaly
$P$ | Period of binary
$\dot{z}_{i}$ |  Radial velocity of star $i$
$K_{i}$ |  Radial velocity semi-amplitude of star $i$
Table 6: Binary system parameters from sections 4.1 and 4.2. Variable | Definition
---|---
$\mu_{+}$, $\mu_{-}$ |  Magnification factor of $+$ and $-$ images
$\boldsymbol{\mathbf{\mu}}$ | Angular velocity
$\boldsymbol{\mathbf{\alpha}}$ | Angular acceleration
$\varpi$ | Parallax
$\theta_{\text{E}}$ | Einstein angle
$\Phi$ | Gravitational potential
$\boldsymbol{\beta}$ | Source-plane location
$\boldsymbol{\varphi}$ | Lensing map
$A_{ij}$, $B_{ij}$ | Macro-lensing Jacobian and its inverse
Table 7: Microlensing parameters from sections 4.3 and 4.5.

##  Appendix B Derivations

###  B.1 Intensity correlation noise

In what follows, we derive expressions for the expectation values of various
moments of the intensity estimator $\hat{I}_{p}(t_{p})$ from eq. 2.43.
Firstly, its time average

$\displaystyle\overline{\hat{I}_{p}(t_{p})}\equiv\frac{1}{t_{\mathrm{obs}}}\int_{t_{p}-t_{\mathrm{obs}}/2}^{t_{p}+t_{\mathrm{obs}}/2}\mathrm{d}t^{\prime}\,\hat{I}_{p}(t^{\prime})$
(B.1)

is expected to equal the true intensity $\langle I_{p}\rangle$, since it is an
unbiased estimator:

$\displaystyle\left<\overline{\hat{I}_{p}(t_{p})}\right>=\frac{\hbar
c\overline{k}}{\eta_{p}A_{p}}\int_{t_{p}-t_{\text{obs}}/2}^{t_{p}+t_{\text{obs}}/2}\frac{\text{d}t^{\prime}}{t_{\text{obs}}}\sum_{i=1}^{N_{p}}\int\frac{\text{d}t_{i}}{t_{\text{obs}}}\frac{e^{\frac{-(t^{\prime}-t_{i})^{2}}{\sigma_{t}^{2}}}}{\sqrt{\pi}\sigma_{t}}=\frac{\hbar
c\overline{k}}{\eta_{p}A_{p}}\frac{N_{p}}{t_{\text{obs}}}\equiv\langle
I_{p}\rangle.$ (B.2)

Again, because $\hat{I}_{p}(t_{p})$ from eq. 2.43 is an unbiased estimator of
the local intensity field, it also faithfully reflects the intensity
_correlations_ at detectors $p=1$ and $p=2$:

$\displaystyle\left<\overline{\hat{I}_{1}(t_{1})\hat{I}_{2}(t_{2})}\right>$
$\displaystyle=\left<\overline{\hat{I}_{1}(t_{1})}\right>\left<\overline{\hat{I}_{2}(t_{2})}\right>\left[1+\underbrace{C(\boldsymbol{\mathbf{d}},\ell,t_{2}-t_{1})}_{\text{correlations}\ll
1}\right]$ (B.3)
$\displaystyle\simeq\frac{\hbar^{2}c^{2}\overline{k}^{2}}{\eta_{1}\eta_{2}A_{1}A_{2}}\frac{N_{1}N_{2}}{t_{\text{obs}}^{2}}=\langle
I_{1}\rangle\langle I_{2}\rangle.$ (B.4)

In the second line, we have approximated the excess intensity correlations to
be small ($C\ll 1$), which is usually a good approximation. Equation B.4
implies that $\hat{C}$ from eq. 2.44 is also an unbiased estimator of $C$.

The statistical error on the time-averaged intensities is negligibly small,
because:

$\displaystyle\left<\bigg{[}\overline{\hat{I}_{p}(t_{p})}\bigg{]}^{2}\right>\simeq\left(\frac{\hbar
c\overline{k}}{\eta_{p}A_{p}}\right)^{2}\left(\prod_{i=1}^{N_{p}}\int\frac{\text{d}t_{i}}{t_{\text{obs}}}\right)\int\frac{\text{d}t^{\prime}}{t_{\text{obs}}}\left(\sum_{i=1}^{N_{p}}\frac{e^{\frac{-(t^{\prime}-t_{i})^{2}}{\sigma_{t}^{2}}}}{\sqrt{\pi}\sigma_{t}}\right)\int\frac{\text{d}t^{\prime\prime}}{t_{\text{obs}}}\left(\sum_{i=1}^{N_{p}}\frac{e^{\frac{-(t^{\prime\prime}-t_{i})^{2}}{\sigma_{t}^{2}}}}{\sqrt{\pi}\sigma_{t}}\right)$
$\displaystyle\simeq\left(\frac{\hbar
c\overline{k}}{\eta_{p}A_{p}}\right)^{2}\left(\frac{N_{p}}{t_{\mathrm{obs}}}\right)^{2}=\langle
I_{p}\rangle^{2}.$ (B.5)

In other words, the fractional variance of $\overline{\hat{I}_{p}(t_{p})}$ is
of order $\mathcal{O}(1/N_{p})$, exceedingly small for very bright sources
with $N_{p}\gg 1$ needed to measure $C$ with high SNR. The (fractional)
statistical error on the intensity correlations is significantly larger, since
it receives contributions only from near-simultaneous photon arrivals. We
find:

$\displaystyle\left<\bigg{[}\overline{\hat{I}_{1}(t)\hat{I}_{2}(t+\tau)}\bigg{]}^{2}\right>=\frac{(\hbar
c\overline{k})^{4}}{(\eta_{1}\eta_{2}A_{1}A_{2})^{2}}\left(\prod_{i=1}^{N_{1}}\int\frac{\text{d}t_{i}^{(1)}}{t_{\text{obs}}}\right)\left(\prod_{m=1}^{N_{2}}\int\frac{\text{d}t_{m}^{(2)}}{t_{\text{obs}}}\right)\int\frac{\text{d}t^{\prime}}{t_{\text{obs}}}\int\frac{\text{d}t^{\prime\prime}}{t_{\text{obs}}}$
$\displaystyle\phantom{=}\times\left[\sum_{i=1}^{N_{1}}\frac{e^{\frac{-\left(t^{\prime}-t_{i}^{(1)}\right)^{2}}{\sigma_{t}^{2}}}}{\sqrt{\pi}\sigma_{t}}\right]\left[\sum_{m=1}^{N_{2}}\frac{e^{\frac{-\left(t^{\prime}-t_{m}^{(2)}\right)^{2}}{\sigma_{t}^{2}}}}{\sqrt{\pi}\sigma_{t}}\right]\left[\sum_{j=1}^{N_{1}}\frac{e^{\frac{-\left(t^{\prime\prime}-t_{j}^{(1)}\right)^{2}}{\sigma_{t}^{2}}}}{\sqrt{\pi}\sigma_{t}}\right]\left[\sum_{n=1}^{N_{2}}\frac{e^{\frac{-\left(t^{\prime\prime}-t_{n}^{(2)}\right)^{2}}{\sigma_{t}^{2}}}}{\sqrt{\pi}\sigma_{t}}\right]$
$\displaystyle=\frac{(\hbar
c\overline{k})^{4}}{(\eta_{1}\eta_{2}A_{1}A_{2})^{2}}\frac{1}{t_{\text{obs}}^{2}}\left(\prod_{i=1}^{N_{1}}\int\frac{\text{d}t_{i}^{(1)}}{t_{\text{obs}}}\right)\left(\prod_{m=1}^{N_{2}}\int\frac{\text{d}t_{m}^{(2)}}{t_{\text{obs}}}\right)\sum_{\begin{subarray}{c}i,j\\\
m,n\end{subarray}}\frac{e^{\frac{-1}{2\sigma_{t}^{2}}\left[\left(t_{i}^{(1)}-t_{m}^{(2)}\right)^{2}+\left(t_{j}^{(1)}-t_{n}^{(2)}\right)^{2}\right]}}{2\pi\sigma_{t}^{2}}$
(B.6) $\displaystyle\simeq\frac{(\hbar
c\overline{k})^{4}}{(\eta_{1}\eta_{2}A_{1}A_{2})^{2}}\frac{1}{t_{\text{obs}}^{2}}\left\\{\left[\sum_{i,m}\int\frac{\text{d}t_{i}^{(1)}}{t_{\text{obs}}}\int\frac{\text{d}t_{m}^{(2)}}{t_{\text{obs}}}\frac{e^{\frac{-1}{2\sigma_{t}^{2}}\left(t_{i}^{(1)}-t_{m}^{(2)}\right)^{2}}}{\sqrt{2\pi}\sigma_{t}}\right]^{2}\right.$
$\displaystyle\phantom{\simeq\frac{(\hbar
c\overline{k})^{4}}{(\eta_{1}\eta_{2}A_{1}A_{2})^{2}}\frac{1}{t_{\text{obs}}^{2}}\Bigg{\\{}\,}+\left.\sum_{i,m}\int\frac{\text{d}t_{i}^{(1)}}{t_{\text{obs}}}\int\frac{\text{d}t_{m}^{(2)}}{t_{\text{obs}}}\frac{e^{\frac{-1}{\sigma_{t}^{2}}\left(t_{i}^{(1)}-t_{m}^{(2)}\right)^{2}}}{2\pi\sigma_{t}^{2}}\right\\}$
$\displaystyle=\frac{(\hbar
c\overline{k})^{4}}{(\eta_{1}\eta_{2}A_{1}A_{2})^{2}}\left[\left(\frac{N_{1}N_{2}}{t_{\text{obs}}^{2}}\right)^{2}+\frac{N_{1}N_{2}}{t_{\text{obs}}^{3}\sqrt{4\pi}\sigma_{t}}\right]$

The latter term in the last equality is responsible for the statistical
variance on $\hat{C}$ in the low-SNR limit for each spectral channel. Putting
everything together, we find:

$\displaystyle\sigma_{\hat{C}}^{2}$
$\displaystyle\equiv\big{\langle}{\hat{C}^{2}}\big{\rangle}-\big{\langle}{\hat{C}}\big{\rangle}^{2}\simeq\frac{\left<\bigg{[}\overline{\hat{I}_{1}(t)\hat{I}_{2}(t+\tau)}\bigg{]}^{2}\right>-\bigg{[}\left<\overline{\hat{I}_{1}(t)\hat{I}_{2}(t+\tau)}\right>\bigg{]}^{2}}{\langle
I_{1}\rangle^{2}\langle
I_{2}\rangle^{2}}=\frac{t_{\text{obs}}}{\sqrt{4\pi}\sigma_{t}}\frac{1}{N_{1}N_{2}}.$
(B.7)

This corresponds to eq. 2.48 from the main text.

###  B.2 Atmospheric aberration

Figure 20: Structure constant $C_{n}^{2}$ for atmospheric refraction from eq.
2.58 as a function of altitude $h$ in the SLC night model.

We derive the mean square doubly-differential refraction phase from
atmospheric aberrations (eq. 2.60 in the main text). We find:

$\displaystyle\left(\sigma_{\phi}^{\text{atm}}\right)^{2}\equiv\Big{\langle}\left[\Delta\widetilde{\phi}^{(1,2)}(\overline{k},\hat{\boldsymbol{\mathbf{\theta}}}_{a},\hat{\boldsymbol{\mathbf{\theta}}}_{b})\right]^{2}\Big{\rangle}_{\mathrm{atm}}$
$\displaystyle=\overline{k}^{2}\iint\mathrm{d}s\,\mathrm{d}s^{\prime}\Big{\langle}\Big{(}n[\boldsymbol{\mathbf{x}}_{a1}(s)]+n[\boldsymbol{\mathbf{x}}_{b2}(s)]-n[\boldsymbol{\mathbf{x}}_{a2}(s)]-n[\boldsymbol{\mathbf{x}}_{b1}(s)]\Big{)}$
$\displaystyle\phantom{=k^{2}\iint\mathrm{d}s\,\mathrm{d}s^{\prime}\Big{\langle}}\times\Big{(}n[\boldsymbol{\mathbf{x}}_{a1}(s^{\prime})]+n[\boldsymbol{\mathbf{x}}_{b2}(s^{\prime})]-n[\boldsymbol{\mathbf{x}}_{a2}(s^{\prime})]-n[\boldsymbol{\mathbf{x}}_{b1}(s^{\prime})]\Big{)}\Big{\rangle}_{\mathrm{atm}}$
$\displaystyle=\overline{k}^{2}\iint_{0}^{L_{\mathrm{atm}}\sec\gamma}\mathrm{d}s\,\mathrm{d}s^{\prime}\,\frac{-C_{n}^{2}(s\cos\gamma)}{2}\Big{\\{}$
(B.8)
$\displaystyle+\big{|}\boldsymbol{\mathbf{x}}_{a1}(s)-\boldsymbol{\mathbf{x}}_{a1}(s^{\prime})\big{|}^{2/3}+\big{|}\boldsymbol{\mathbf{x}}_{b1}(s)-\boldsymbol{\mathbf{x}}_{b1}(s^{\prime})\big{|}^{2/3}+\big{|}\boldsymbol{\mathbf{x}}_{a2}(s)-\boldsymbol{\mathbf{x}}_{a2}(s^{\prime})\big{|}^{2/3}+\big{|}\boldsymbol{\mathbf{x}}_{b2}(s)-\boldsymbol{\mathbf{x}}_{b2}(s^{\prime})\big{|}^{2/3}$
$\displaystyle-\big{|}\boldsymbol{\mathbf{x}}_{a1}(s)-\boldsymbol{\mathbf{x}}_{b1}(s^{\prime})\big{|}^{2/3}-\big{|}\boldsymbol{\mathbf{x}}_{b1}(s)-\boldsymbol{\mathbf{x}}_{a1}(s^{\prime})\big{|}^{2/3}-\big{|}\boldsymbol{\mathbf{x}}_{a2}(s)-\boldsymbol{\mathbf{x}}_{b2}(s^{\prime})\big{|}^{2/3}-\big{|}\boldsymbol{\mathbf{x}}_{b2}(s)-\boldsymbol{\mathbf{x}}_{a2}(s^{\prime})\big{|}^{2/3}$
$\displaystyle-\big{|}\boldsymbol{\mathbf{x}}_{a1}(s)-\boldsymbol{\mathbf{x}}_{a2}(s^{\prime})\big{|}^{2/3}-\big{|}\boldsymbol{\mathbf{x}}_{b1}(s)-\boldsymbol{\mathbf{x}}_{b2}(s^{\prime})\big{|}^{2/3}-\big{|}\boldsymbol{\mathbf{x}}_{a2}(s)-\boldsymbol{\mathbf{x}}_{a1}(s^{\prime})\big{|}^{2/3}-\big{|}\boldsymbol{\mathbf{x}}_{b2}(s)-\boldsymbol{\mathbf{x}}_{b1}(s^{\prime})\big{|}^{2/3}$
$\displaystyle+\big{|}\boldsymbol{\mathbf{x}}_{a1}(s)-\boldsymbol{\mathbf{x}}_{b2}(s^{\prime})\big{|}^{2/3}+\big{|}\boldsymbol{\mathbf{x}}_{b1}(s)-\boldsymbol{\mathbf{x}}_{a2}(s^{\prime})\big{|}^{2/3}+\big{|}\boldsymbol{\mathbf{x}}_{a2}(s)-\boldsymbol{\mathbf{x}}_{b1}(s^{\prime})\big{|}^{2/3}+\big{|}\boldsymbol{\mathbf{x}}_{b2}(s)-\boldsymbol{\mathbf{x}}_{a1}(s^{\prime})\big{|}^{2/3}\Big{\\}}.$

To get to the second equality, we plugged in eq. 2.57. To get to the last
equality, we used eqs. 2.58 and 2.59. We integrate to the “height”
$L_{\mathrm{atm}}$ of the atmosphere, though technically $C_{n}^{2}(h)$ drops
to zero smoothly. In the last four lines, we have expanded the 16 correlation
functions into suggestive groups. The last two lines are not correct if $d\gg
L_{0}$, and we will ignore them in what follows; in that regime, the sum of
the last two lines (which involve cross-correlations between paths to
_different_ detectors) is very strongly suppressed. Each of the four terms in
the first line are equal to $[(s-s^{\prime})^{2}]^{1/3}$, while those in the
second line are of the form
$[(\hat{\boldsymbol{\mathbf{\theta}}}_{a}s-\hat{\boldsymbol{\mathbf{\theta}}}_{b}s^{\prime})^{2}]^{1/3}$.
Putting this all together gives:

$\displaystyle\left(\sigma_{\phi}^{\text{atm}}\right)^{2}$
$\displaystyle\simeq
2\overline{k}^{2}\iint_{0}^{L_{\mathrm{atm}}\sec\gamma}\mathrm{d}s\,\mathrm{d}s^{\prime}\,C_{n}^{2}(s\cos\gamma)\Big{\\{}(s^{2}-2ss^{\prime}\cos\theta_{ba}+s^{\prime
2})^{1/3}-(s-s^{\prime})^{1/3}\Big{\\}}$ $\displaystyle\simeq
2\overline{k}^{2}\int_{0}^{L_{\mathrm{atm}}\sec\gamma}\text{d}s\,C_{n}^{2}(s\cos\gamma)\underbrace{\int_{-s}^{+s}\text{d}x\,\Big{\\{}\big{[}x^{2}+s(s+x)\theta_{ba}^{2}\big{]}^{1/3}-\big{[}x^{2}\big{]}^{1/3}\Big{\\}}}_{\sqrt{\pi}\frac{\Gamma(-5/6)}{\Gamma(-1/3)}\theta_{ba}^{5/3}-2\theta_{ba}^{2}+\mathcal{O}(\theta_{ba}^{11/3})}.$
(B.9)

To get to the last line, we used the small-angle approximation, and used the
fact that most of the integral has support at $|x|=|s-s^{\prime}|\lesssim s$
to separate the integrals. This result, combined with the definition of the
isoplanatic patch angle in eq. 2.61, gives eq. 2.60. For reference, we show
the dependence of the structure constant $C_{n}^{2}(h)$ on altitude $h$ in
figure 20.

##  Appendix C Wave Optics

In this appendix, we carry out an analysis using scalar wave optics for the
propagation of starlight within each telescope. Emphasis is given on analytic
results to understand the requirements for an EPIC intensity interferometer.
We include finite telescope apertures and finite pixel sizes in our analysis,
as well as atmospheric and instrumental aberrations. Some simplifications and
idealizations are necessarily made, and we leave a detailed numerical study to
future work.

The structure of this appendix follows the depiction of our setup in figure 5:
section C.1 treats the primary and secondary mirrors, section C.2 the optical
elements of the geometric path-delay system, section C.3 the effect of
aberrations on image formation, and finally section C.4 presents the details
of the calculation of section 3.4 concerning the seeing angle suppression for
large telescopes and fast photodetectors resulting from the spectroscopic
instrument.

###  C.1 Primary and secondary mirrors

The focus of this section is the derivation of the spatial profile of the
electric field, after collimation from the secondary mirror (as depicted in
figure 5(a)), in the presence of mirror distortions. Our main result here is
the functional dependence of the electric field profile on the focal ratio of
the two mirrors, which magnifies both stellar separations and wavefront
aberrations induced by the primary mirror.

Let two point-like stars, $a$ and $b$, emit monochromatic spherical waves of
wavenumber $k$. Right before incidence on the parabolic mirror of a telescope,
the electric field is given by

$\begin{split}E_{k}(\mathbb{x},t)&=Ae^{i\omega
t}e^{-ik|\mathbb{r}_{a}+\mathbb{x}|+i\phi^{\text{em}}_{a}+i\tilde{\phi}^{\text{atm}}_{a}(\mathbb{x},t)}+Be^{i\omega
t}e^{-ik|\mathbb{r}_{b}+\mathbb{x}|+i\phi^{\text{em}}_{b}+i\tilde{\phi}^{\text{atm}}_{b}(\mathbb{x},t)}\\\
&\simeq e^{i\omega
t}\left(Ae^{-ik(r_{a}+\mathbb{x}\cdot\hat{\boldsymbol{\theta}}_{a})+i\phi^{\text{em}}_{a}+i\tilde{\phi}^{\text{atm}}_{a}(\mathbb{x},t)}+Be^{-ik(r_{b}+\mathbb{x}\cdot\hat{\boldsymbol{\theta}}_{b})+i\phi^{\text{em}}_{b}+i\tilde{\phi}^{\text{atm}}_{b}(\mathbb{x},t)}\right),\end{split}$
(C.1)

where $A$ and $B$ are the amplitudes of the electric field of each star on the
telescope, $\mathbb{r}_{i}$ is the physical location of star $i$ on the sky,
$\mathbb{x}$ a 2D vector on a plane parallel to the telescope’s aperture
(denoted $\boldsymbol{\mathbf{b}}$ in section 2), centered at the vertex of
the parabolic mirror, and
$\hat{\boldsymbol{\theta}}_{i}\equiv\mathbb{r}_{i}/r_{i}$ the unit vector
pointing to the position of star $i$. Because $\mathbb{x}$ is a 2D vector, we
can trade $\hat{\boldsymbol{\theta}}_{a}$ for its component parallel to
$\mathbb{x}$, which we denote $\boldsymbol{\theta}_{a}$ without the hat.
Clearly the telescope is pointed such that $|\boldsymbol{\theta}_{a}|\ll 1$.
The phase $\phi_{i}^{\text{em}}$ is the random phase of the emitter $i$ and
$\tilde{\phi}^{\text{atm}}_{i}(\mathbb{x},t)$ is the space- and time-dependent
phase the wave from emitter $i$ inherits from passing through turbulent layers
of the atmosphere. In the second line we neglected wavefront curvature effects
$\sim k|\mathbb{x}|^{2}/r\ll 1$. In what follows we will also drop the
explicit time-dependence $e^{i\omega t}$.

We take the primary mirror to have diameter $D$ and centered at
$\mathbb{x}=0$, with a surface given by

$z_{\text{primary}}(\mathbb{x})=\frac{|\mathbb{x}|^{2}}{4f}+\frac{1}{2}\tilde{z}^{\text{inst}}(\mathbb{x}),\quad|\mathbb{x}|\in[0,D],$
(C.2)

where the first term corresponds to an ideal spherical paraboloid of focal
length $f$ and $\tilde{z}^{\text{inst}}$ describes its deformations, which
induce instrumental (as opposed to atmospheric) aberrations on the wavefront.
We will assume that their effect is achromatic.

Using the Kirchhoff-Helmholtz diffraction integral [111], the electric field
from star $a$ at a distance $z\gg(D/2)^{2}/(4f)$ from the primary mirror will
be given by

$\begin{split}E_{a,k}(\mathbb{x}^{\prime},z)=&A\frac{-ik}{2\pi
z}\exp\left[-ikr_{a}+i\phi^{\text{em}}_{a}+ik\left(z+\frac{x^{\prime
2}}{2z}\right)\right]\\\
&\times\int\text{d}^{2}x\,\exp\left\\{ik\left[-\frac{x^{2}}{2}\left(\frac{1}{f}-\frac{1}{z}\right)-\mathbb{x}\cdot\left(\frac{\mathbb{x}^{\prime}}{z}+\boldsymbol{\theta}_{a}\right)+\tilde{z}^{\text{inst}}(\mathbb{x})\right]+i\tilde{\phi}_{a}^{\text{atm}}(\mathbb{x},t)\right\\},\end{split}$
(C.3)

where $\mathbb{x}^{\prime}$ is the 2D vector on a surface perpendicular to
propagation. The spatial profile, as described by the second line, is
intuitive: the first term describes smearing due to defocusing, which vanishes
at the focal plane $z=f$, while the second term would give the usual
diffraction limited electric field profile centered at
$-\boldsymbol{\theta}_{a}$ on the focal plane, in the absence of aberrations.

Analogous considerations apply for the secondary parabolic mirror that
collimates the beam. We find at a distance $z$ from the vertex of the
secondary mirror

$\begin{split}E_{a,k}^{\text{coll}}(\mathbb{x}^{\prime},z)&=A\exp\left[-ikr_{a}+i\phi^{\text{em}}_{a}+ik(z+f-f^{\prime})\right]\frac{f}{f^{\prime}}\frac{kF^{-1}(z)}{\pi}\\\
&\int\text{d}^{2}x\exp\left\\{ik\left[F^{-1}(z)\left(\mathbb{x}-\frac{f}{f^{\prime}}\mathbb{x}^{\prime}\right)^{2}+\mathbb{x}\cdot\boldsymbol{\theta}_{a}+\tilde{z}^{\text{inst}}(\mathbb{x})\right]+i\tilde{\phi}_{a}^{\text{atm}}(\mathbb{x},t)\right\\},\end{split}$
(C.4)

where
$F^{-1}(z)\equiv\frac{f^{\prime}}{2f}\frac{1}{z\frac{f}{f^{\prime}}+f-f^{\prime}}$
and $\mathbb{x}^{\prime}$ is, again, the 2D vector on a surface perpendicular
to propagation. In deriving this result we have assumed for simplicity that
the secondary mirror is an ideal paraboloid, large enough to collect all the
light from the primary and not induce any further diffraction effects. In what
follows, we suppress the atmospheric phase for brevity, as it can formally be
treated in the same way as instrumental aberrations.

We wish to gain some analytic understanding of the electric field profile that
exits the two-mirror system, in the presence of potentially large (compared to
the wavelength) aberrations, i.e. in the limit $k\tilde{z}^{\text{inst}}\gg
1$. This can be done with the stationary phase approximation.

We define the exponent to be
$g(\mathbb{x})\equiv\left[F^{-1}(z)\left(\mathbb{x}-\frac{f}{f^{\prime}}\mathbb{x}^{\prime}\right)^{2}+\mathbb{x}\cdot\boldsymbol{\theta}_{a}+\tilde{z}^{\text{inst}}(\mathbb{x})\right]$.
In the ideal mirror case $\tilde{z}^{\text{inst}}(\mathbb{x})=0$ and the
equations $\partial_{x}g=0$ and $\partial_{y}g=0$ have the simple solutions
(critical points) $x=x^{\prime}f/f^{\prime}$ and $y=y^{\prime}f/f^{\prime}$.
Since $|\mathbb{x}|\in[0,D/2]$, it follows that
$|\mathbb{x}^{\prime}|\in[0,Df^{\prime}/(2f)]$, so the collimated beam exists
over a smaller spatial extent, given by the geometric optics width of the beam
at the plane of the secondary mirror.

We may then approximate

$\int\text{d}^{2}x\,\exp[ikg(\mathbb{x})]\simeq\begin{cases}\pi
F(z)\exp\left(ikf\mathbb{x}^{\prime}\cdot\boldsymbol{\theta}_{a}/f^{\prime}\right)/k&\text{for~{}}|\mathbb{x}|\in\left[0,\frac{f}{f^{\prime}}\frac{D}{2}\right]\\\
0&\text{elsewhere}\end{cases},$ (C.5)

so that the collimated beam has modulus
$|E(\mathbb{x}^{\prime},z)|=Af/f^{\prime}$ and its extent is
$|\mathbb{x}^{\prime}|\in[0,Df^{\prime}/(2f)]$. The key takeaway is that, on
the $\mathbb{x}^{\prime}$ plane, the angle $\boldsymbol{\theta}_{a}$ appears
magnified by the magnification $f/f^{\prime}$, as expected from geometric
optics. Due to the prefactor $f/f^{\prime}$ of eq. C.4, the energy of the beam
is conserved.

In the presence of mirror imperfections $\tilde{z}^{\text{inst}}$, critical
points cannot always be found analytically, but we may work in the limit where
they induce small corrections to the shape of the collimated beam. As an
example, in this limit the critical points in the $x$ direction are

$\begin{split}x\simeq\frac{f}{f^{\prime}}x^{\prime}-\frac{F(z)}{2}\left[\theta_{a}^{(x)}+\left.\frac{\partial\tilde{z}^{\text{inst}}}{\partial
x}\right|_{x=\frac{f}{f^{\prime}}x^{\prime}}\right],\end{split}$ (C.6)

where we took $\partial_{x}\tilde{z}^{\text{inst}}$ to be a perturbation on
the zeroth order solution $x\simeq x^{(0)}=x^{\prime}f/f^{\prime}$ and
$\theta_{a}^{(x)}$ is the $x$-coordinate of $\boldsymbol{\theta}_{a}$. One can
work out optical tolerances for the primary mirror due to this approximation,
but they turn out to be less tight compared to those of section C.3 that
concern pixel sizes. Thus, we do not discuss them here.

Having established the above, we can approximate the collimated beam for both
stars as

$\begin{split}E_{k,\text{coll}}(\mathbb{x}^{\prime},z)&\simeq\frac{f}{f^{\prime}}\exp\left[ik(2(f-f^{\prime})-z)\right]\exp\left\\{ik\left[\frac{F(d)}{4}\left.\left(\nabla\tilde{z}^{\text{inst}}\right)^{2}\right|_{\mathbb{x}=\frac{f}{f^{\prime}}\mathbb{x}^{\prime}}+\tilde{z}^{\text{inst}}\left(\frac{f}{f^{\prime}}\mathbb{x^{\prime}}\right)\right]\right\\}\\\
&\bigg{[}A\exp\left[-ikr_{a}+i\phi^{\text{em}}_{a}+i\tilde{\phi}_{a}^{\text{atm}}\left(\frac{f}{f^{\prime}}\mathbb{x}^{\prime},t\right)\right]\exp\left\\{ik\left[\frac{f}{f^{\prime}}\mathbb{x^{\prime}}\cdot\boldsymbol{\theta}_{a}-\frac{F(d)\boldsymbol{\theta}_{a}^{2}}{4}\right]\right\\}\\\
&+B\exp\left[-ikr_{b}+i\phi^{\text{em}}_{b}+i\tilde{\phi}_{b}^{\text{atm}}\left(\frac{f}{f^{\prime}}\mathbb{x}^{\prime},t\right)\right]\exp\left\\{ik\left[\frac{f}{f^{\prime}}\mathbb{x^{\prime}}\cdot\boldsymbol{\theta}_{b}-\frac{F(d)\boldsymbol{\theta}_{b}^{2}}{4}\right]\right\\}\bigg{]},\end{split}$
(C.7)

where $|\mathbb{x}^{\prime}|\in[0,Df^{\prime}/(2f)]$, $z$ is the coordinate
distance with the primary mirror’s vertex being at $z=0$, and
$d=f-f^{\prime}-z$ is the (positive) distance from the secondary mirror. In
deriving this expression we have dropped the dependence of
$\tilde{z}^{\text{inst}}$ on
$\boldsymbol{\theta}_{a},\,\boldsymbol{\theta}_{b}$, as it is small. To
simplify notation hereafter we will denote the full wavefront aberration
imprinted by the primary-secondary mirror system as
$\tilde{z}^{\text{inst}}\left(f\mathbb{x}^{\prime}/f^{\prime}\right)$, instead
of the full exponent in the first line above, and drop the
$\propto\boldsymbol{\theta}^{2}$ terms.

The most important feature of our result, eq. C.7, is the appearance of the
magnification factor $f/f^{\prime}$ multiplying the physical positions of the
stars and the arguments of the aberration functions, both instrumental and
atmospheric. This is not a surprising result: just as relative stellar
separations $\boldsymbol{\theta}_{a}-\boldsymbol{\theta}_{b}$ are magnified by
$f/f^{\prime}$, so does the seeing angle due to the rescaled argument of
$\tilde{\phi}^{\text{atm}}$, and the mirror-induced wavefront aberrations due
to $\tilde{z}^{\text{inst}}$.

###  C.2 The EPIC optical system

The EPIC optical delay module has the form of a Mach-Zehnder interferometer of
unequal path lengths, achieved by positioning the first beamsplitter at an
angle $\gamma\neq 45^{\text{o}}$, as shown in figure 5(b) and described in
section 2. This system’s role is twofold: it both induces the geometric path
delay that expands the field of view, and corrects for the wavefront tilt
between the images of the two stars. The latter suppresses fringe contrast if
uncorrected, an effect known since the advent of intensity interferometry [34]
and relevant when each telescope _itself_ resolves the angular separation
being measured.

In section C.2 we study the induced phases by the multiple reflections within
an idealized beamsplitter and technical requirements when operating it with
broadband light. In section C.2 we generalize eqs. 2.17 and 2.26 for the
geometric path delay and wavefront correction for a general 2D relative
stellar separation and derive expressions for the required two-dimensional
wavefront correction angle. We conclude by commenting on tolerances relevant
to the EPIC module.

#### Fourier decomposition

The analysis in this section will be carried out with plane waves for clarity.
This retains all generality because any wave can be decomposed into plane
waves via a Fourier transform, and the effect of all optical elements is
linear so that no two wavelengths mix. As such, any electric field
$E(\mathbb{r})$, where $\mathbb{r}=(x,y,z)$, can be decomposed as

$\begin{split}E(\mathbb{r})\equiv\int\text{d}^{3}p\,\tilde{E}(\mathbb{p})e^{-i\mathbb{p}\cdot\mathbb{r}},\quad\tilde{E}(\mathbb{p})=\int\frac{\text{d}^{3}x}{(2\pi)^{3}}E(\mathbf{r})e^{+i\mathbb{p}\cdot\mathbb{r}}.\end{split}$
(C.8)

Each Fourier component has a fixed direction and so it is easier to calculate
how it transforms after propagation through the EPIC optical system of figure
5(b). Formally, the direction will be modified such that

$\tilde{E}(\mathbb{p})e^{-i\mathbb{p}\cdot\mathbb{r}}\to\tilde{E}(\mathbb{p})e^{-i\mathbb{M}(\mathbb{r)}\cdot\mathbb{p}},$
(C.9)

for some function
$\mathbb{M}(\mathbb{r})\equiv\left(M_{1}(x,y,z),M_{2}(x,y,z),M_{3}(x,y,z)\right)$
that depends on the optical instruments the wave encounters.We can resum these
plane waves into $E^{\prime}$ at some point after beam-recombination at the
exit of the EPIC module so that

$E^{\prime}(x,y,z)=E(M_{1},M_{2},M_{3}),$ (C.10)

which therefore corresponds to translations and distortions of the shape of
the incoming collimated beam.

#### Beamsplitters

Here we study an idealized beamsplitter, consisting of a single slab of
material with index of refraction $n_{2}$. Though this study is necessarily
incomplete, we wish to highlight some non-trivial requirements for these
optical elements when operated with broadband light. Our main result is that,
when accounting for multiple reflections, there are irreducible phases that
may extinguish the intensity of certain wavelengths. We defer more complete
modeling to future work.

Let a ray propagate in a medium of index of refraction $n_{1}$ (air) and be
incident on a semi-transparent plate of index of refraction $n_{2}$, at an
angle $\gamma$. If the angle of refraction is $\theta$, such that
$n_{1}\sin\gamma=n_{2}\sin\theta$ by Snell’s law, the reflection and
transmission coefficients at each surface of the plate for the parallel555We
will neglect polarization effects throughout the discussion here and derive
expressions only for the parallel component. component of the electric field
are [111]

$\text{$n_{1}\to n_{2}$:}\quad
t=\frac{2n_{1}\cos\gamma}{n_{2}\cos\gamma+n_{1}\cos\theta},\quad\quad
r=\frac{n_{2}\cos\gamma-n_{1}\cos\theta}{n_{2}\cos\gamma+n_{1}\cos\theta}$
(C.11) $\text{$n_{2}\to n_{1}$:}\quad
t^{\prime}=\frac{2n_{2}\cos\theta}{n_{1}\cos\theta+n_{2}\cos\gamma}=\frac{n_{2}\cos\theta}{n_{1}\cos\gamma}t,\quad\quad
r^{\prime}=\frac{n_{1}\cos\theta-
n_{2}\cos\gamma}{n_{1}\cos\theta+n_{2}\cos\gamma}=-r,$ (C.12)

where $r^{2}+tt^{\prime}=1$, but $r^{2}+t^{2}\neq 1$. Naively, the transmitted
beam suffers a $tt^{\prime}$ suppression, compared to the reflected which
suffers only an $r$ suppression. Consistency with energy conservation for
elements with similar $r$ and $t$ requires taking into account multiple
reflections within the beamsplitter [111]. If $E^{(i)}$ is the amplitude of
the incident field, then the reflected and transmitted amplitudes are,
respectively

$\begin{split}E^{(r)}=\sqrt{\mathbb{R}}\frac{1-e^{i\delta}}{1-\mathbb{R}e^{i\delta}}E^{(i)},\quad\quad
E^{(t)}=\frac{\mathbb{T}e^{i\Phi}}{1-\mathbb{R}e^{i\delta}}E^{(i)},\end{split}$
(C.13)

where $\mathbb{R}\equiv r^{2}$ and $\mathbb{T}\equiv tt^{\prime}$, so that
$\mathbb{R+T}=1$, and the phases accumulated are $\delta=2k\ell
n_{2}\cos\theta$, where $\ell_{\text{b}}$ is the width of the beamsplitter,
and $\Phi=k\ell_{\text{b}}(n_{2}\cos\theta-n_{1}\cos\gamma)$.

The reflected and transmitted intensities are

$\begin{split}\left|\frac{I^{(r)}}{I^{(i)}}\right|=\frac{F\sin^{2}\frac{\delta}{2}}{1+F\sin^{2}\frac{\delta}{2}}\equiv|\tilde{\mathbb{R}}|,\quad\quad\left|\frac{I^{(t)}}{I^{(i)}}\right|=\frac{1}{1+F\sin^{2}\frac{\delta}{2}}\equiv|\tilde{\mathbb{T}}|,\end{split}$
(C.14)

where $F\equiv\frac{4\mathbb{R}}{(1-\mathbb{R})^{2}}$, and
$|\tilde{\mathbb{R}}|+|\tilde{\mathbb{T}}|=1$ recovers energy conservation. In
the main text we denoted $|\mathbb{\tilde{T}}|\equiv\epsilon^{2}$ for
notational simplicity. The intensity of the recombined beam will scale as
$|\tilde{\mathbb{R}}|^{2}|\tilde{\mathbb{T}}|^{2}=\left[|\tilde{\mathbb{R}}|(1-|\tilde{\mathbb{R}}|)\right]^{2}$,
which is maximized for $|\tilde{\mathbb{R}}|=1/2$ as shown in the main text.

However, the intensity of the reflected and transmitted beams depends on the
wavelength through the phase $\delta$, so that $|\tilde{\mathbb{R}}|=1/2$ can
strictly be achieved for only one wavelength. In order to maximize the
intensity in as broad a range of wavelengths as possible, we require both
$|\tilde{\mathbb{R}}|^{2}|\tilde{\mathbb{T}}|^{2}=1/4$ at a particular
wavelength, and that the second derivative
$\frac{\text{d}^{2}}{\text{d}\delta^{2}}|\tilde{\mathbb{R}}|^{2}|\tilde{\mathbb{T}}|^{2}$
vanishes, so that the change around the preferred $\delta$ is as “slow” as
possible. This can be achieved formally for $F=1$, since the second derivative
scales as $\propto(1-F)$. As seen in the left panel of figure 21, for $F=1$
the final intensity is equal to $1/4$ for a factor of 2 in wavelength, but
decreases significantly outside this range. If we choose $\delta=\pi$ to
correspond to the central wavelength we are recording, e.g.
$\delta(450\text{nm})=\pi$, then we can cover the region $300-600$nm with
intensity loss $\lesssim 1.7\%$.

Figure 21: Left: Intensity of the main EPIC fringe normalized to total
incident intensity, when accounting for multiple reflections within each
beamsplitter, as a function of the phase
$\delta=2k\ell_{\text{b}}n_{2}\cos\theta$ accumulated after one crossing of
the beamsplitter. The product of the reflection and transmission coefficients
plotted on the vertical axis is another definition of the quantity
$\epsilon^{2}(1-\epsilon^{2})$ used in e.g. eq. 2.19. The lines, in order of
diminishing thickness, correspond to increasing values of the ratio $F\equiv
4\mathbb{R}/(1-\mathbb{R})^{2}$, as defined below eq. C.14. The vertical
dotted lines correspond to the angles $\pi\pm 1$, within which there is a
coverage of a factor of 2 in wavelength. An $F=1$ beamsplitter has optimal
response to broadband light, as proven in the text. Right: Required index of
refraction $n_{2}$ as a function of the incidence angle $\gamma$, for
different values of $F$. The index of refraction of air is $n_{1}=1$ and of
the beamsplitter $n_{2}$. The thickness grading is the same as in the left
figure.

Three comments are in order in light of the the above results. First, and most
importantly, the beamsplitter will be large compared to the wavelength, so
that $\ell_{\text{b}}\gg\lambda$. Figure 21 (left) is periodic in $\delta$
with period $2\pi$, so we can choose in general
$\delta(450\text{nm})=(2m+1)\pi$, $m\in\mathbb{N}$, and then the width of the
plate needs to be

$\ell_{\text{b}}=\frac{2m+1}{n_{2}\sqrt{1-\frac{n_{1}^{2}}{n_{2}^{2}}\sin^{2}\gamma}}\frac{\lambda}{4},\quad
m\in\mathbb{N},$ (C.15)

where we have neglected the wavelength dependence of the index of refraction.
If one chooses $m>1$ for $\lambda$, the above considerations no longer apply
as wavelengths $\lambda^{\prime}$ which are $\mathcal{O}(1)$ different from
$\lambda$ (as needed for a broadband setup) will acquire a $\mathcal{O}(1)$
different phase $\delta(\lambda^{\prime})$, which is not within $\pi$ of
$\delta(\lambda)$: for instance $\lambda/2$ acquires a phase $m\pi/2$ which
can be very much different from $m\pi-\pi/2$ for large $m$. The broadband
response of wide beamsplitters will have to be determined experimentally and
might entail some $\mathcal{O}(1)$ loss of light. In the limit
$\ell_{\text{b}}\gg\lambda$, approximately $83\%$ of the spectrum will exit
the beamsplitter with more than half its incident intensity and the total
intensity lost will be at most $35\%$. Alternatively, by energy conservation,
we can retain the full intensity if we made use of both readouts of the second
beamsplitter.

Second, the choice of $F$ imposes a specific combination of the incidence
angle and the index of refraction, since $F$ is a proxy for the reflection
coefficient. In particular, $r^{2}=1-\frac{2}{F}\left(\sqrt{1+F}-1\right),$
and so the required incidence angle as a function of the indices of refraction
and $F$ is given by

$\cos\gamma=\left[\frac{1-\frac{n_{1}^{2}}{n_{2}^{2}}}{\frac{n_{2}^{2}}{n_{1}^{2}}\left(\frac{1-r}{1+r}\right)^{2}-\frac{n_{1}^{2}}{n_{2}^{2}}}\right]^{1/2}.$
(C.16)

On the right panel of figure 21 we are plotting the required index of
refraction of the plate $n_{2}$ as a function of the incidence angle $\gamma$,
for different values of $F$ and assuming $n_{1}=1$ for air. From this figure
it is clear that, at least for this toy model, there might be a need for
materials with large indices of refraction.

Third, the above discussion concerned intensities, but it is also necessary to
know the $\lambda$-dependent phase shift of the electric fields due to each
beamsplitter, as it may induce smearing of the intensity fringe we are
observing. We find

$\begin{split}E^{(r)}&=\left(\frac{F\sin^{2}\frac{\delta}{2}}{1+F\sin^{2}\frac{\delta}{2}}\right)^{1/2}e^{i\Phi_{r}}E^{(i)},\quad\tan\Phi_{r}\equiv-\frac{(1-\mathbb{R})\sin\delta}{(1+\mathbb{R})(1-\cos\delta)}\\\
E^{(t)}&=\frac{e^{i(\Phi_{t,1}+\Phi_{t,2})}}{\left(1+F\sin^{2}\frac{\delta}{2}\right)^{1/2}}E^{(i)},\quad\tan\Phi_{t,1}\equiv\frac{\mathbb{R}\sin\delta}{1-\mathbb{R}\cos\delta}\quad\text{and}\quad\Phi_{t,2}=\frac{\delta}{2}-k\ell_{\text{b}}n_{1}\cos\gamma.\end{split}$
(C.17)

Assuming that the beamsplitters are identical, the electric field that exits
the EPIC module will have accumulated a phase
$\Phi_{\text{total}}=\Phi_{r}+\Phi_{t,1}+\Phi_{t,2}$. For $m>1$ this phase
oscillates rapidly for a finite range of $\lambda$ and may reduce fringe
contrast. The strictest tolerance concerns the size of the beamsplitter, which
needs to be such that $\sigma_{k}\ell_{\text{b}}\lesssim 1$, so that fringe
contrast is not lost within each pixel, and must be known with sub-wavelength
precision, so that the relative phase delay between pixels can be accounted
for. The former sets $\ell_{\text{b}}\sim$mm and the latter can be determined
experimentally.

To summarize, in this section we have studied the broadband response of a
beamsplitter made entirely of a single index of refraction. We have identified
an irreducible phase contribution that may reduce the intensity of the beam
exiting the EPIC module by a third at most, but can be recovered if both
readouts of the EPIC module are used. Given our simplistic modeling, we
anticipate that experimental work will be needed to optimize these optical
elements.

#### Geometric delay and wavefront correction

In this section, we derive the general expression for the geometric path delay
in the case of a relative 2D wavefront angle. We derive the expressions for
the required angles of the EPIC system’s mirrors, and discuss their
tolerances.

In order to achieve the wavefront correction described in section 2, we can
tune the upper mirror of the path-extension system at an angle
$\gamma^{\prime}\neq\gamma$, as shown in figure 5(b). Because the upper mirror
might need to correct for a 2D angle, we must include a small angle
$\delta^{\prime}$ with the axis that in figure 5(b) is perpendicular to the
page. We define as $\ell_{\perp}$ the distance between the first beamsplitter
and the mirror of the upper part of the EPIC delay system when
$\gamma^{\prime}=\gamma$ and $\delta^{\prime}=0$. This is an unambiguous
length convenient for this calculation. It can be shown that
$w_{p}=\ell_{\perp}\sin 2\gamma/\cos\gamma$, where $w_{p}$ is the length shown
on figure 5(b).

In the limit $|\gamma-\gamma^{\prime}|\equiv|\varepsilon|\ll 1$666Not to be
confused with the transmission coefficient $\epsilon^{2}$ we defined in
section 2. and $\delta^{\prime}\ll 1$ we find for the $\mathbb{M}(\mathbb{r})$
function of eq. C.9

$\begin{split}M_{1}&\simeq x+y\delta^{\prime}\sin
2\gamma+2z(\gamma-\gamma^{\prime})-2\ell_{\perp}\left[(\gamma-\gamma^{\prime})\cos\gamma+\sin\gamma\right]\\\
M_{2}&\simeq-x\delta^{\prime}\sin 2\gamma+y-z\delta^{\prime}(1+\cos
2\gamma)+2\ell_{\perp}\delta^{\prime}\cos\gamma\\\
M_{3}&\simeq-2(\gamma-\gamma^{\prime})x+y\delta^{\prime}(1+\cos
2\gamma)+z-2\ell_{\perp}\left[\cos\gamma-(\gamma-\gamma^{\prime})\sin\gamma\right].\end{split}$
(C.18)

These expressions have a simple geometric meaning: for a plane wave only in
the $z$-direction and for $\gamma=\gamma^{\prime}$ and $\delta^{\prime}=0$,
the center of the beam is shifted by $2\ell_{\perp}\sin\gamma$ in the
$x$-direction, and gains a phase because of the geometric path between points
with the same $z$-coordinate, $2\ell_{\perp}\cos\gamma$.

Neglecting instrumental and atmospheric aberrations, the recombined collimated
beam exiting the second beamsplitter is then

$\begin{split}E_{k}^{\text{recomb}}(\mathbb{x},z)\simeq&\tilde{\mathbb{R}}\tilde{\mathbb{T}}\frac{f}{f^{\prime}}\exp\left[ik(2(f-f^{\prime})-z)\right]\\\
&\bigg{[}Ae^{-ikr_{a}+i\phi^{\text{em}}_{a}}e^{ikx(2\varepsilon+\frac{f}{f^{\prime}}\theta_{a,x})}e^{iky(-2\delta^{\prime}\cos^{2}\gamma+\frac{f}{f^{\prime}}\theta_{a,y})}e^{2ik\ell_{\perp}(\cos\gamma-\varepsilon\sin\gamma-\frac{f}{f^{\prime}}\theta_{a,x}\sin\gamma)}\\\
&+Be^{-ikr_{b}+i\phi^{\text{em}}_{b}}e^{ikx(2\varepsilon+\frac{f}{f^{\prime}}\theta_{b,x})}e^{iky(-2\delta^{\prime}\cos^{2}\gamma+\frac{f}{f^{\prime}}\theta_{b,y})}e^{2ik\ell_{\perp}(\cos\gamma-\varepsilon\sin\gamma-\frac{f}{f^{\prime}}\theta_{b,x}\sin\gamma)}\bigg{]}\\\
&+e^{2ik\ell_{\perp}(\sin\gamma+\varepsilon\cos\gamma)}\bigg{[}Ae^{-ikr_{a}+i\phi^{\text{em}}_{a}}e^{ik\frac{f}{f^{\prime}}\mathbb{x}\cdot\boldsymbol{\theta}_{a}}e^{-2ik\ell_{\perp}\frac{f}{f^{\prime}}\theta_{a,x}\sin\gamma}\\\
&\hskip
99.58464pt+Be^{-ikr_{b}+i\phi^{\text{em}}_{b}}e^{ik\frac{f}{f^{\prime}}\mathbb{x}\cdot\boldsymbol{\theta}_{b}}e^{-2ik\ell_{\perp}\frac{f}{f^{\prime}}\theta_{b,x}\sin\gamma}\bigg{]},\end{split}$
(C.19)

where we have expanded the exponents to first order in $\varepsilon$ and
$\delta^{\prime}$. The overall coefficient
$\tilde{\mathbb{R}}\tilde{\mathbb{T}}$ has to do with the beamsplitters and
expressions for a simple model can be found in section C.2. We assume that
these are identical for both beamsplitters.

The principal EPIC fringe arises from the terms in the third and fourth line
of eq. C.19. The relative phase between these two terms is

$\Delta\phi=ik\left[x\left(2\epsilon+\frac{f}{f^{\prime}}\theta_{ba,x}\right)+y\left(-2\delta\cos^{2}\gamma+\frac{f}{f^{\prime}}\theta_{ba,y}\right)\right]+ik\ell_{p}$
(C.20)

where we defined
$\boldsymbol{\theta}_{ba}=(\theta_{ba,x},\theta_{ba,y})\equiv\boldsymbol{\theta}_{b}-\boldsymbol{\theta}_{a}$
and $\ell_{p}$ is the geometric delay, the tuning of which allows for the
expansion of the field of view (see section 2.2). It is given by

$\ell_{p}=\frac{\ell_{\perp}\cos(\gamma-\gamma^{\prime})}{\cos\gamma^{\prime}}\left[\cos
2\gamma+\cos 2(\gamma-\gamma^{\prime})-\sin 2\gamma-\sin
2(\gamma-\gamma^{\prime})\right]-2\ell_{\perp}\frac{f}{f^{\prime}}\theta_{{ba,x}}\sin\gamma,$
(C.21)

where the last term $\propto\theta_{ba,x}$ is in fact $\propto\varepsilon$, so
this expression reduces to eq. 2.17 for $\gamma=\gamma^{\prime}$. The
wavefront correction for large apertures and large stellar separations amounts
to canceling the $x$\- and $y$\- dependence of $\Delta\phi$ by choosing

$\varepsilon\equiv\gamma-\gamma^{\prime}=-\frac{f}{2f^{\prime}}\theta_{ba,x}\quad\text{and}\quad\delta^{\prime}=\frac{f}{2f^{\prime}\cos^{2}\gamma}\theta_{{ba,y}}.$
(C.22)

#### Wavefront correction and telescope pointing tolerances

The residual angle $\varphi$ from the aforementioned tuning has to be within
the diffraction limit of the telescope, so that
$\varphi\lesssim(kD/2)^{-1}\simeq 8\,$mas at $\lambda=500\,$nm and for a
$4\,$m diameter primary. The angle has to be both known and tuned to within
this precision. The former is satisfied by the first step of the source
localization protocol outlined in 3.1, while the latter sets an instrumental
tolerance which should be achievable.

Furthermore, we have found a term proportional to (a component of) the
magnified relative angle of the two stars in eq. C.21. We remind the reader
that $\ell_{p}$ must be tuned close to $\mathbb{d}\cdot\boldsymbol{\theta}$,
where $\mathbb{d}$ is the baseline vector and $\boldsymbol{\theta}$ the
relative stellar separation, to within a wavelength, as discussed in sections
2.2 and in particular 3.1 and 3.2. Even though the magnification factor
$f/f^{\prime}$ can be large, $\ell_{\perp}$ will only be on the order of tens
of cm. As a result, this factor does not induce any stricter tolerances than
what we have discussed so far.

Finally, we factored out an overall $\mathbb{x}$ dependent factor
$\exp\left(ikf\mathbb{x}\cdot\boldsymbol{\theta}_{a}/f^{\prime}\right)$. This
factor includes possible telescope pointing errors and, although similar is
form to the relative wavefront tilt, it does _not_ lead to fringe suppression.
Pointing of the telescope affects the light paths of both stars in the same
way, so, in fact, it is largely irrelevant: all it does is shift the overall
position of the image on the photodetector array, without any suppression of
correlation. The only tolerance on pointing is that any error remains within
the angle subtended by each pixel. Seeing fluctuations demand that this angle
be at least $1$ arcsec, as we will see in section C.4, and stable pointing
with such a precision is achievable even with commercial star-guiding systems
for amateur astronomy.

###  C.3 Optical aberrations

In this section, we derive tolerance conditions for our optical systems. We
first sketch a proof of the well-known fact that a periodic grating produces
an identical image of the source at every wavelength, so it is sufficient to
neglect the dispersion of the grating and, rather, study the effects of
aberrations on a single image of a point-like source. We derive tolerances on
primary aberrations analytically and show that the requirements for intensity
interferometry in general, and for EPIC in particular, can be orders of
magnitude looser than those of traditional imagers and amplitude
interferometers. Towards the end, we revisit atmospheric aberrations, which we
studied in sections 2.5 and B.2, this time focusing on the size of the finite
seeing disk on the pixel array. We defer a complete numerical study to future
work.

#### Multiple image formation by a diffraction grating

The presence of the diffraction grating complicates the analytic understanding
of the 2D image formation at the focal plane in the presence of optical
aberrations. Fortunately, we can leverage the well-known fact of spectroscopy
that a grating produces a copy of the image at every wavelength and the
angular separation of these images on the focal plane is the same as that of
the aberration-free images of point-like sources.

We may consider, for simplicity, a 1D grating of periodicity $d_{\text{g}}$
and linear size along the “spectroscopic direction” $w=Nd_{\text{g}}$, where
$N$ is the number of illuminated slits of width $s<d$, and a wave
$E_{k}\propto\exp\\{ik[z+\tilde{z}(x)]\\}$ at normal incidence with wavefront
aberrations $\tilde{z}$. Here, $z$ is the direction of propagation and $x$ is
the transverse coordinate, namely the one on the plane of the grating. The
outgoing light is focused on the photodetector array using a camera of focal
length $f_{\text{c}}$. The amplitude of the wave $E_{k}$ at the position
$\theta_{\text{F}}=x_{\text{F}}/f_{\text{c}}$ on the focal plane is

$\begin{split}E_{k}(\theta_{\text{F}})\propto&\sum_{n=-N/2+1}^{N/2}\int_{(2n-1)\frac{d_{\text{g}}}{2}-\frac{s}{2}}^{(2n-1)\frac{d_{\text{g}}}{2}+\frac{s}{2}}\text{d}\xi\,e^{-ik\xi\theta_{\text{F}}}e^{ik\tilde{z}(\xi)}\\\
&=\sum_{n}e^{ik\theta_{\text{F}}(2n-1)d_{\text{g}}/2}\int_{-\frac{s}{2}}^{\frac{s}{2}}\text{d}\xi^{\prime}\,e^{-ik\left[\xi^{\prime}\theta_{\text{F}}-\tilde{z}\left(\xi^{\prime}+(2n-1)d_{\text{g}}/2\right)\right]}.\end{split}$
(C.23)

Because $d>s$ we can ignore the $\xi^{\prime}$ dependence in the argument of
$\tilde{z}$ to get

$\begin{split}E_{k}(\theta_{\text{F}})\propto&s\frac{\sin\left(k\theta_{\text{F}}\frac{s}{2}\right)}{k\theta_{\text{F}}\frac{s}{2}}\sum_{n}e^{ik\theta_{\text{F}}(2n-1)d_{\text{g}}/2}e^{ik\tilde{z}\left((2n-1)d_{\text{g}}/2\right)}.\end{split}$
(C.24)

Neglecting the single-slit prefactor, which is effectively constant for
$\theta_{\text{F}}\ll 1/(ks)$, the remainder is periodic in
$\theta_{\text{F}}$ with period $T=4\pi/(kd_{\text{g}})$. Its modulus has half
that period, i.e $|E(\theta_{\text{F}}+mT/2)|=|E(\theta_{\text{F}})|$. Thus
the profile centered around any diffraction order $m\in\mathbb{Z}$ is the
same.

Within a period, we can take the limit $d_{\text{g}}\to 0$. Then we can
approximate $(2n-1)d_{\text{g}}/2\simeq x$, such that
$x\in[-Nd_{\text{g}}/2,Nd_{\text{g}}/2]=[-w/2,w/2]$, for $N\gg 1$. We
specialize to the first diffraction peak $m=1$, to find

$E_{k}(\theta_{\text{F}})\propto\int_{-w}^{w}\text{d}x\,e^{-ikx\left(\theta_{\text{F}}-\frac{2\pi}{kd_{\text{g}}}\right)}e^{ik\tilde{z}(x)}.$
(C.25)

This is the image of the source in the focal plane as seen through an aperture
of width $w$ which imprints aberrations $\tilde{z}$. The salient feature of
this expression is that the angle where each wavelength is focused at,
$2\pi/(kd_{\text{g}})$, depends on the wavelength as one expects from a
dispersion grating, but the functional form of the wave profile is the same
around that point. Prefactors can be determined from energy conservation.

The argument presented here disentangles the effects of aberrations from the
dispersion of the grating and allows us to consider only the aberrated
intensity distribution of a monochromatic wave in what follows.

#### Tolerances on primary aberrations

Intensity interferometry is not sensitive to the precise functional form of
the intensity distribution of an image. Instead, the only requirement is that
all the light is collected by the photodetector. In the simplest setup, we
allocate one pixel per $\sigma_{k}$ wavenumbers, so, according to the proof of
the previous subsection, we require the spread of a single image to be within
this one pixel. We provide a general analytic framework for deriving the
spread of the distorted image and calculate tolerances on Seidel aberrations.
Our discussion follows ref. [111] and is technical. We refer the reader
interested in how (primary) tolerances for EPIC compare to those of regular
diffraction-limited telescopes to figure 22 and especially table 8.

Let $\mathbb{x}=(x,y)$ be the 2D cartesian vector on the plane of the
photodetector, which lies at $z=z_{\text{p}}$ from a lens of width
$w=Nd_{\text{g}}$ and focal length $f_{\text{c}}$ centered at
$(\mathbb{x},z)=(\mathbb{0},0)$.777The lens here is a proxy for the grating.
In the case of a grating there is the additional distance from the grating to
the camera mirror/lens. The width of the lens is related to the aperture of
the primary by the magnification factor $w\equiv Df^{\prime}/f$, as we have
proven in section C.1, e.g. see eq. C.7. We define the corresponding polar
coordinates $(r,\phi)$ for the focal plane and $(\rho,\theta)$ for the lens.