# Precision speckle pattern reconstruction for high contrast imaging

Dotan Gazith  Barak Zackay  Dept. of Particle Phys. & Astrophys., Weizmann
Institute of Science, Rehovot 76100, Israel

###### Abstract

In High Contrast Imaging, a large instrumental, technological and algorithmic
effort is made to reduce residual speckle noise and improve the detection
capabilities. In this work, we explore the potential of using a precise
physical description of speckle images, in conjunction with the optimal
detection statistic to perform High Contrast Imaging. Our method uses short-
exposure speckle images, reconstructing the Point Spread Function (PSF) of
each image with phase retrieval algorithms. Using the reconstructed PSF’s we
calculate the optimal detection statistic for all images. We analyze the
arising bias due to the use of a reconstructed PSF and correct for it
completely up to its accumulation over $10^{4}$ images. We measure in
simulations the method’s sensitivity loss due to overfitting in the
reconstruction process and get to an estimated 5$\sigma$ detection limit of
$5\times 10^{-7}$ flux ratio at angular separations of
$0.1-0.5^{\prime\prime}$ for a $1h$ observation of Sirius A with a
2m-telescope.

High contrast techniques (2369), Astronomy data analysis (1858), Astronomical
seeing (92)

##  1 Introduction

High Contrast Imaging (HCI), the method to detect exoplanets in which the
exoplanet is seen directly as an additional source near its host star, is an
important method for current and future characterization of exoplanets
atmospheres (Currie et al. (2011), Konopacky et al. (2013)) as it measures
light from the exoplanet itself.

HCI requires separating the faint exoplanet from its bright nearby host star,
which is highly dependent on the telescope’s angular resolution. For ground-
based observatories, the atmosphere presents an additional constraint on the
telescope’s angular resolution. The turbulent flow of air in the atmosphere
creates index of refraction variations, that create rapidly changing phase
aberrations. The phase aberrations degrade the telescope’s angular resolution,
creating speckle images for short-exposure imaging and a broad seeing disk
image for long-exposure.

A way to overcome atmospheric seeing is the use of Adaptive Optics (AO)
systems that sense and correct, in real-time, the atmospheric phase
aberrations. Such systems bring the angular resolution from the seeing limit
close to the telescope’s diffraction scale and are used by the leading HCI
instruments on state-of-the-art telescopes (VLT-SPHERE Beuzit et al. (2019),
Subaru-CHARIS Groff et al. (2015) and MEC Walter et al. (2020), Gemini-GPI
Nielsen et al. (2019)).

The biggest challenge for those AO-based HCI instruments are residual
aberrations not corrected or sensed by the AO system Mawet et al. (2012),
those aberrations create quasi-static speckles, slowly varying speckles that
are confused as point sources. A wide range of efforts are made to tackle this
challenge, better wavefront sensing (eg. Baudoz et al. (2005), Skaf et al.
(2021), Skaf et al. (2022)), deformable mirror technological improvements (eg.
Madec (2012)), nulling coronagraphy (eg. Ruane et al. (2018)) and post-
processing methods (eg. Marois et al. (2006), Racine et al. (1999), Lafreniere
et al. (2007), Rodack et al. (2021), Frazin & Rodack (2021)).

Analyzing the idealized case of HCI with perfectly known, yet uncorrected,
Point Spread Function (PSF), we obtain that even a modest-sized telescope in
the seeing limited case can reach a fantastic contrast, comparable to the best
performance of state-of-the-art facilities, as can be seen in the dashed lines
in Figure 1.

Motivated by this computation, in this work, we propose a method to perform
HCI using short exposure speckle images, our scheme is illustrated in Figure
2, the measurements are both imaging and wavefront sensing, a sequence of
algorithms is used to reconstruct the atmospheric phase aberrations from the
measurements 4.1, the optimal detection statistic is then calculated 2.3 and
Parametric Bootstrap method is used to correct the bias arising in the
procedure 4.2. Analysis of the expected performance is presented in 2.3 and
4.3 and the results presented in Figure 1.

Figure 1: $5\sigma$ Detection limit at different separations. The detection
limit is calculated for a $2\mathrm{m}$-telescope, at seeing condition
$r_{0}=15\mathrm{cm}$, for total of $10^{14}$ photons (corresponds to
observing Sirius A for $\sim 1\mathrm{h}$ at 500-700nm with efficiency of
$\sim 70\%$). The colored-dotted lines are the optimal detection limit
described in Section 2, the colored-solid include overfitting losses described
in Section 4, the black line shows the performance of GPI from Ruffio et al.
(2017). The plot was created using the package by Bailey (2021).  Figure 2:
Schematic diagram of the HCI method. Schematically in our method, we take
short exposure images, reconstruct the phase aberrations for each image,
calculate the statistic using the image and reconstructed phase and sample to
correct for its bias.

##  2 High contrast imaging through the atmosphere

When imaging from the ground, the atmosphere presents a difficulty. As air
flows turbulently, different regions contain air of different velocities and
densities that cause spatial variations in the index of refraction. As light
passes through this medium with an in-homogeneous index of refraction, optical
path differences accumulate, and phase aberrations are presented on the pupil
plane of our telescope. Phase aberrations degrade the PSF of our telescope,
and the atmospheric phase aberrations change rapidly, the instantaneous
degraded PSF is typically comprised of diffraction-limited speckles, called
speckle image, and the time-averaged called the seeing disk.

###  2.1 Atmospheric seeing

The PSF is a characteristic of the optical system that can be calculated for
monochromatic illumination as the Fourier transform of the aperture function
($B$) and the phase aberration ($\phi$)

$P=\left|\mathcal{F}\left[Be^{i\phi}\right]\right|^{2}\,.$ (1)

In the case of perfect imaging from a circular aperture, the PSF will take the
shape of an Airy disk with an angular scale of the diffraction limit
$\lambda/D$. The phase aberrations caused by the atmosphere degrade the PSF
significantly and are a topic of many theoretical and observational studies.
As a model for the atmospheric phase aberrations, we will use the classical
Kolmogorov phase aberration structure function with a typical scale $r_{0}$
called the Fried parameter

$\displaystyle D_{\phi}\left(\Delta\vec{x}\right)$ $\displaystyle\coloneqq$
$\displaystyle\langle\left(\phi(\vec{r})-\phi(\vec{r}+\vec{\Delta
x})\right)^{2}\rangle$ (2) $\displaystyle=$ $\displaystyle
6.88\left(\frac{|\Delta\vec{x}|}{r_{0}}\right)^{5/3}\,.$

The Fried parameter can also be understood as the size of the telescope that
transitions from being diffraction limited to seeing limited.

To model imaging in a finite wavelength band, we integrate the single
wavelength PSFs incoherently and include the two leading effects that change
the PSF as a function of the wavelength:

  * •

The diffraction scale increases $\lambda_{0}/D\to\lambda/D$

  * •

The aberration in phase units decreases $\phi\to\phi\lambda_{0}/\lambda$

$P=\int_{\lambda_{l}}^{\lambda_{h}}E_{\lambda}\left|\mathcal{S}_{\lambda/\lambda_{0}}\mathcal{F}\left[Be^{i\phi\frac{\lambda_{0}}{\lambda}}\right]\right|^{2}d\lambda\,,$
(3)

where $E_{\lambda}$ is a term representing the spectrum of the source, and
$\mathcal{S}_{\lambda/\lambda_{0}}$ an operator that stretches the coordinates
as the diffraction scale increases.

###  2.2 Statistical model

To push the limits of direct imaging we have to stay as close as possible to
optimality, therefore we will start by mathematically formulating the high
contrast imaging detection problem rigorously.

Directly imaging exoplanets can be formulated as a hypothesis testing question
in which we want to distinguish between the null hypothesis of the images made
of a single star with some flux $f$

$H_{0}:\,T=f\delta_{0}\,,$ (4)

and the alternative hypothesis of an additional source with flux fraction
$\epsilon$ at relative position $\vec{q}$

$H_{1}:\,T=f\left((1-\epsilon)\delta_{0}+\epsilon\delta_{\vec{q}}\right)\,.$
(5)

The mapping between the true sky $T$ and the measured images is by the PSF,
with additive noise we will assume has a shot-noise component related to the
source, shot-noise related to the sky background and detector read noise,
taking the Gaussian approximation of the Poisson distribution (justified as we
have a high number of photons per pixel)

$I\sim\mathcal{N}\left(P\otimes T,P\otimes T+b^{2}\right)\,,$ (6)

where $\mathcal{N}$ is the normal distribution, $b$ is the noise from sky
background and read-noise, $\otimes$ is the convolution operator.

###  2.3 Detection statistic

To distinguish between the two hypotheses we will use the following statistic
(with further discussion in Appendix A)

$\displaystyle S[q]$ $\displaystyle=$
$\displaystyle\left(\overleftarrow{P_{\hat{\phi}}}\otimes\frac{I-fP_{\hat{\phi}}}{fP_{\hat{\phi}}+b^{2}}\right)_{\vec{q}}$
$\displaystyle-$
$\displaystyle\left(\overleftarrow{P_{\hat{\phi}}}\otimes\frac{I-fP_{\hat{\phi}}}{fP_{\hat{\phi}}+b^{2}}\right)_{0}\,;$
$\displaystyle\hat{\phi}$ $\displaystyle=$
$\displaystyle\mathrm{argmax}\\{P\left(I|H_{0},P_{\phi}\right)P(\phi)\\}\,$
(7)

where $\overleftarrow{P}$ is the reverse of $P$, and ${\hat{\phi}}$ is the
Maximum-A-Posteriori estimator of the atmospheric phase aberrations $phi$. We
can read this statistic as a filter matching the image’s deviation from a
point source, and subtracting the location of the primary source as the
unknown flux would create a flux deficit at its location. And the variance
associated with the image’s variance is

$\displaystyle V_{S}=$ $\displaystyle f^{2}$
$\displaystyle\Bigg{(}\overleftarrow{P^{2}_{\hat{\phi}}}\otimes\frac{1}{fP_{\hat{\phi}}+b^{2}}$
(8) $\displaystyle-$ $\displaystyle
2\overleftarrow{P^{2}_{\hat{\phi}}}\otimes\frac{P_{\hat{\phi}}}{fP_{\hat{\phi}}+b^{2}}+\left[\overleftarrow{P^{2}_{\hat{\phi}}}\otimes\frac{1}{fP_{\hat{\phi}}+b^{2}}\right]_{0}\Bigg{)}\,.$

We can estimate the $5\sigma$ detection threshold of this statistic for the
ideal case in which we know $\phi$ exactly

$\min{\epsilon}\,:\,5\leq\frac{E\left[S|H_{1}\right]}{\sqrt{E\left[V_{S}\right]}}\,,$
(9)

by simulating atmospheric phase screens from the power spectrum calculated by
Noll (1976), the resulting contrast thresholds as a function of separation for
different bandwidths are shown in the colored dashed lines in Figure 1.

##  3 Simulations

To test our method we used end-to-end numerical simulations of the proposed
measurement instruments.

We generate atmospheric aberrations by sampling the Fourier modes of the
Kolmogorov spectrum as derived by Noll (1976). Images are calculated with the
assumed true sky image and PSF as calculated in Equation 3, and then sampled
according to Poisson with additional Gaussian noise. The simulation of the
wavefront sensor neglects chromatic behavior, as only the centroid position
will be used.

Reference values for simulation parameters are listed in Table 1.

Parameter | Value
---|---
Pixel grid | $256\times 256$
Aperture diameter (physical) | $2m$
Aperture diameter (pixels) | $85\frac{1}{3}px$
Nyquist oversampling ratio | $1.5$
Fried parameter | $15cm$
Flux (per image, finite bandwidth) |  $10^{7}$ @ $500-700nm$
Background/read noise | 1 photon per pixel
Table 1: Reference parameters for simulations.

##  4 Methods

In order to use this statistic, we first need to show a method to calculate
$\hat{\phi}$, and even though we know how the statistic distributes given the
correct $\phi$, we need to examine carefully how it distributes when using
$\hat{\phi}$.

###  4.1 Recovering the atmospheric phase aberrations

In order to calculate the statistic from Equation 2.3 we have to find
$\mathrm{argmax}\\{P\left(I|H_{0},P_{\phi}\right)P(\phi)\\}$ , the Maximum-a-
Posteriori (MAP) estimator for the atmospheric phase aberrations.

The general problem of phase retrieval, recovering a phase from its Fourier
modulus in the presence of noise is known to be hard and a topic of many
studies (for example the reviews Fienup (2013), Shechtman et al. (2015)).
Therefore we employ a simple yet powerful sequence of algorithms. Starting
from direct measurement of the phase using a Shack-Hartmann WaveFront Sensor
(SHWFS, Platt & Shack (2001)), then improving that estimator using the
Gerchberg-Saxton (GS, Gerchberg & Saxton (1972)) algorithm and finally using
Gradient Descent to ensure we get the MAP.

This procedure converges to the correct MAP for most instances in the case of
enough photons ($\sim 10^{6}$ per image), reasonable atmosphere
($r_{0}=15\mathrm{cm}@D=2\mathrm{m}$), but only for imaging in a single
wavelength, further discussion, and elaborate results are presented in
Appendix B.

###  4.2 Detection in the presence of learned PSF - $H_{0}$

When applying the statistic from Equation 2.3 on _learned PSF_ we have to deal
with the effect of its inevitable errors. From simulations we learn, as
expected, that our phase estimator can be modeled as distributing normally
with some bias $\mu$ and covariance $\Sigma$ around the correct atmospheric
phase aberration

$\hat{\phi}\sim\mathcal{N}\left(\phi+\mu,\Sigma\right)\,.$ (10)

From Equations 3 and 2.3 we can see that our statistic is a non-linear
function of the phase estimator. Therefore, our estimator for the statistic,
calculated using our estimator of the phase, is slightly _biased_

$\Delta S_{\phi}\equiv E\left[S_{\hat{\phi}}-S_{\phi}\right]\neq 0\,,$ (11)

but as it is a bias shared by all images, when we accumulate the statistic
over many images it accumulates too and put an upper limit on the number of
images we can use once it becomes significant.

As the notation suggests, this bias depends on $\phi$, and for simplicity, we
will assume it varies slowly relative to our phase estimation error and
therefore can be expanded in a series

$\displaystyle\Delta S_{\hat{\phi}}[q]=\Delta S_{\phi}[q]$ $\displaystyle+$
$\displaystyle\frac{d\Delta S_{\phi}[q]}{d\phi}\Delta\phi$ (12)
$\displaystyle+$ $\displaystyle\frac{1}{2}\Delta\phi\frac{d^{2}\Delta
S_{\phi}[q]}{d\phi^{2}}\Delta\phi\,,$

and its expectation value

$\displaystyle E_{\hat{\phi}}\left[\Delta S_{\hat{\phi}}[q]\right]=\Delta
S_{\phi}[q]$ $\displaystyle+$ $\displaystyle\frac{d\Delta
S_{\phi}[q]}{d\phi}\mu$ (13) $\displaystyle+$
$\displaystyle\frac{1}{2}\mu\frac{d^{2}\Delta S_{\phi}[q]}{d\phi^{2}}\mu$
$\displaystyle+$ $\displaystyle\frac{1}{2}\mathrm{tr}\left(\frac{d^{2}\Delta
S_{\phi}[q]}{d\phi^{2}}\Sigma\right)\,.$

We extend the method of Parametric Bootstrap Dekking (2005) to estimate this
bias. By sampling images and solving their corresponding phases based on the
recovered phase estimator we can sample the distribution in Equation 10

$\hat{\hat{\phi}}_{1},\hat{\hat{\phi}}_{2}\sim\mathcal{N}\left(\hat{\phi}+\mu,\Sigma\right)\,,$
(14)

and calculate the following linear combination that has the expectation value
as the bias we want to correct (detailed calculation in Appendix C)

$E\left[\Delta S_{2\hat{\phi}-\hat{\hat{\phi}}_{1}}-\Delta
S_{2\hat{\hat{\phi}}_{1}-\hat{\phi}}+\Delta
S_{\hat{\hat{\phi}}_{1}+\hat{\hat{\phi}}_{2}-\hat{\phi}}\right]=\Delta
S_{\phi}\,.$ (15)

Figure 3: Detection statistic (a) without and (b) with bias correction. The
images were simulated for a $2\mathrm{m}$-telescope, at seeing condition
$r_{0}=15\mathrm{cm}$, for a bandwidth of 10%, and the statistic was
accumulated over $10^{4}$ images.

In practice, we calculate this term a few times to get better convergence to
the mean and have a handle on its variance, which is empirically smaller than
the inherent variance of the statistic, that is calculated in Equation 8, and
can be reduced as $1/\sqrt{N}$ with more simulations.

###  4.3 Signal loss due to over fitting

After we made sure we are recovering null detections for images sampled from
$H_{0}$, we need to test the expected performance of the method for signals
from $H_{1}$.

An inevitable tension is present in our phase aberration estimation process.
The posterior probability of the phase aberration, $\phi$, is dominated by the
deviation of its corresponding PSF, $P_{\phi}$, from the image. This will
lead, in the presence of an additional source, to _overfitting_ of the
secondary as part of the primary. Detecting using an overfitted $\hat{\phi}$
subtracts part of the secondary and therefore loses some of the signal
available with perfect knowledge of $\phi$. To quantify this effect, we
perform injection-recovery simulations for different separations and different
brightness.

Figure 4: Injection-Recovery simulation. (a) The recovered ratio for different
injected SNR per image, injection at $|q|=0.2^{\prime\prime}$, we can see a
plateau as we go for weaker signals, (b) recovery ratio for different
separations at $\mathrm{SNR_{inj}=1}$ per image.

In general, we expect the recovered SNR to be some function of the injected
SNR and the separation from the star at which the signal is injected. In the
weak signal limit at which we work, the injected SNR per image is small
$\mathrm{SNR_{inj}}\ll 1$, and we expect the recovered SNR to be proportional
to the injected SNR, as can be seen in Figure 4a, and the dependence on the
location to be some smooth function that takes into account the assumed phase
aberrations spectrum and the degeneracy between the phase aberrations and the
additional source, as can be seen in Figure 4b.

$\mathrm{SNR_{inj}}=\alpha(|q|)\mathrm{SNR_{inj}}\,.$ (16)

We use the values of $\alpha(|q|)$ from Figure 4 to convert the optimal dotted
lines in Figure 1 to the more realistic solid lines. In the future, we will
investigate using several wavelength bands in tandem to mitigate this loss,
using the a-chromatic position of the secondary and the chromatic effect of
the phase aberrations.

##  5 Conclusion and outlook

We presented a method to reconstruct the atmospheric phase aberrations using
phase retrieval algorithms, and use its corresponding PSFs to calculate the
optimal statistic to detect the presence of a secondary faint source in a
sequence of short-exposure speckle images. We presented a method to correct
the arising bias in our statistic, showing its complete correction for stacks
of up to $10^{4}$ images, and a pathway to increase it by expanding Equation
12 beyond second order, which is crucial for the use of the detection method.

The leading idea for our method is to optimize the PSF of the image through
parameters of a physically constrained model, use the optimal statistic, and
correct for its arising bias or other artifacts. The analysis presented in
this work treats the case of imaging with no AO, but the method isn’t limited
to such observations. A probable avenue for improvement is to use AO to some
extent which will reduce the photon noise of the host star, and concentrate
the light of the planet to a diffraction-limited spot, thus enabling even
better detection of the secondary source.

To advance this method to achieve high contrast on sky, further work is
required, including:

  * •

Better phase retrieval scheme needs to be devised, to enable the retrieval
from images taken in a narrow wavelength band.

  * •

Demonstration of the method in-lab is essential and will ensure we can
precisely describe speckle images.

  * •

Generalizing the method for simultaneous imaging in a few bands, which is
important for collecting enough light, and might make the phase retrieval more
regular and therefore easier.

## Acknowledgements

BZ is supported by a research grant from the Center for New Scientists at the
Weizmann Institute of Science and a research grant from the Ruth and Herman
Albert Scholarship Program for New Scientists. BZ and DG are supported by the
Israeli Council for Higher Education: Competitive Program for Data Science &
AI Research Centers.

##  Appendix A Approximated Statistic

As shown by Neyman & Pearson (1933), when testing hypotheses the likelihood
ratio test achieves maximal detection power for a given false positive
probability, therefore we will employ it for our detection problem

$\Lambda\left(I\right)=\frac{P\left(I|H_{1}\right)}{P\left(I|H_{0}\right)}\,,$
(A1)

which we can marginalize over the atmospheric aberrations

$\Lambda\left(I\right)=\frac{\int
P\left(I|H_{1},P_{\phi}\right)P(\phi)d\phi}{\int
P\left(I|H_{0},P_{\phi}\right)P(\phi)d\phi}\,.$ (A2)

Calculating explicitly both marginalized likelihoods is impractical so we
employ the Gaussian approximation. At the limit of well-measured phase
aberrations, in which we are interested, the posterior distribution of the
phase aberrations is well approximated as a Gaussian.

Under the Gaussian approximation, the integrals reduce to the value at maximum
times the square root of its covariance determinant. We further approximate
the MAP estimate under $H_{1}$ as the MAP estimate under $H_{0}$. This
approximation is justified as in a single image we cannot favor $H_{1}$, so
specifically, we cannot detect a difference in the MAP.

The covariance also doesn’t contribute, as it presents a logarithmic
correction while being nearly identical for the null and alternative
hypotheses.

$\displaystyle\Lambda\left(I\right)$ $\displaystyle\approx$
$\displaystyle\frac{P\left(I|H_{1},P_{\hat{\phi}}\right)P(\hat{\phi})}{P\left(I|H_{0},P_{\hat{\phi}}\right)P(\hat{\phi})}\,;$
$\displaystyle\hat{\phi}$ $\displaystyle=$
$\displaystyle\mathrm{argmax}\\{P\left(I|H_{0},P_{\phi}\right)P(\phi)\\}\,.$
(A3)

and for the case of a companion detection with known PSF, as we assume here,
the likelihood-ratio test can be simplified to the expression in Equation 2.3,
as calculated by Nir et al. (2019).

##  Appendix B Phase recovery procedure

###  B.1 Shack-Hartmann Wavefront Sensor

The SHWFS is a device that images small areas of the aperture separately. For
sufficiently small areas (of the order of $r_{0}$) the atmospheric phase
aberrations can be estimated as simply linear. When imaging a wavefront with
some constant slope the image is an off-centered point, with the tip-tilt of
the point related linearly to the slope of the wavefront.

In our simulation, we used an array of $26\times 26$ sub-apertures which for
our reference seeing conditions and telescope size resulted in each sub-
aperture being of the size $\sim r_{0}/2\times r_{0}/2$, to recover the phase
aberrations from the WFS we used a maximum-a-posteriori estimator in the
Gaussian regime and using a measured design matrix, as done by Clare (2004).

Figure 5: Comparison between a (a) sampled PSF and the (b) the PSF
corresponding to the phase measured by the SHWFS. An obvious, order-unity
difference is visible between the two PSFs, therefore the measured phase is
not accurate enough to be used for detection.

The phase aberrations recovered by the WFS have significant errors, with the
PSF calculated using them having order-unity deviations from the input PSF as
can be seen in Figure 5, therefore cannot be used to calculate the statistic.
To improve the phase estimator we use further optimization techniques.

###  B.2 Gerchberg-Saxton algorithm

The Gerchberg-Saxton algorithm is an algorithm that was designed and proved to
solve the phase retrieval problem under the $l_{2}$ norm by Gerchberg & Saxton
(1972). We started the algorithm with the initial guess of the pupil function
and the phase measured by the WFS, from that the algorithm transforms the
guessed field to the image plane, constraining the amplitude by the measured
image, then returning to the pupil plane constraining the pupil function again
and repeating until convergence or some maximal number of iterations.

To quantify the convergence of the algorithm we compare an observable,
likelihood-based criteria

$\sum\frac{\left(I-fP\right)^{2}}{fP+b^{2}}\lessgtr N_{\text{pixel}}\,,$ (B1)

and an unobservable, phase error criteria, we simulated images in the single
wavelength case for different seeing and flux conditions we examine the
performance of the algorithm in the observable criteria and its agreement with
the unobservable one.

$r_{0}[\mathrm{cm}]$ $N_{\mathrm{photons}}$ | $10^{5}$ | $3\times 10^{5}$ | $10^{6}$ | $3\times 10^{6}$
---|---|---|---|---
15 |  100 |  90 |  67 |  65
18 |  100 |  98 |  93 |  89
21 |  100 |  100 |  98 |  96
24 |  100 |  100 |  100 |  100
Table 2: Success rate percent for different seeing conditions and photons per
speckle image (in the used sub-band). The success rate was determined
according to the likelihood threshold for $10^{3}$ attempts.

In Table 2 we show the fraction of phase aberrations that are successful in
converging to a solution that passes the likelihood threshold, we see that the
algorithm converges more for smaller phase aberrations or flux. However, in
Table 3 we see the phase error of runs that converged successfully in the
likelihood sense is high for the small flux case, which we can understand as
not constraining enough. To ensure the gradient-based method will work we want
the phase error to be smaller than a radian so the problem will be close to
quadratic or some low-order polynomial.

$r_{0}[\mathrm{cm}]$ $N_{\mathrm{photons}}$ | $10^{5}$ | $3\times 10^{5}$ | $10^{6}$ | $3\times 10^{6}$
---|---|---|---|---
15 |  0.79 |  0.45 |  0.16 |  0.05
18 |  0.81 |  0.43 |  0.18 |  0.06
21 |  0.71 |  0.41 |  0.19 |  0.07
24 |  0.80 |  0.40 |  0.19 |  0.08
Table 3: Phase root-mean-square for different seeing conditions and photons
per image. The phase root-mean-square was averaged for runs that passed the
likelihood threshold.

In future works we will expand this optimization step to allow convergence for
finite bandwidth imaging.

###  B.3 Gradient descent

To go from the $l_{2}$ solution to the MAP estimator we optimize a function
proportional to the log-posterior

$\mathcal{L}(\phi)\coloneqq\sum_{x,y}\frac{\left(I-fP_{\phi}\right)^{2}}{fP_{\phi}+b^{2}}+\phi^{T}C^{-1}\phi\,,$
(B2)

using gradient descent with the method to rapidly calculate the gradient as
done by Fienup (1999), this method takes advantage of the Fourier transform in
Equation 3 to calculate it with only 2 FFT operations

$\frac{\partial\mathcal{L}}{\partial\phi_{u}}=\left[2\phi^{T}C^{-1}\right]_{u}+\int_{\lambda_{l}}^{\lambda_{h}}\Re{\left\\{-iB\frac{\lambda_{0}}{\lambda}e^{-i\phi\frac{\lambda_{0}}{\lambda}}\mathcal{F}^{-1}\left[2f\left(\mathcal{S}^{-1}_{\lambda/\lambda_{0}}\frac{\partial\mathcal{L}_{I}}{\partial
P}\right)\mathcal{F}\left[Be^{i\phi\frac{\lambda_{0}}{\lambda}}\right]\right]\right\\}_{u}}d\lambda\,.$
(B3)

##  Appendix C Second order bias estimate

To prove Equation 15 we start from $\hat{\hat{\phi}}_{1}$ and
$\hat{\hat{\phi}}_{2}$ that are estimates for $\hat{\phi}$ and i.i.d according
to $\mathcal{N}\left(\hat{\phi}+\mu,\Sigma\right)$, examining each term
separately, first an unbiased phase estimate

$\displaystyle 2\hat{\phi}-\hat{\hat{\phi}}_{1}$ $\displaystyle=$
$\displaystyle\hat{\phi}+(\hat{\phi}-\hat{\hat{\phi}}_{1})$
$\displaystyle\sim$ $\displaystyle\mathcal{N}\left(\phi,2\Sigma\right)$
$\displaystyle E_{\hat{\phi},\hat{\hat{\phi}}_{1}}\left[\Delta
S(2\hat{\phi}-\hat{\hat{\phi}}_{1})\right]$ $\displaystyle=$
$\displaystyle\Delta S(\phi)+\mathrm{tr}\left(\Sigma\frac{\partial^{2}\Delta
S}{\partial\phi^{2}}|_{\phi}\right)\,,$ (C1)

and two biased phase estimates with different variance

$\displaystyle 2\hat{\hat{\phi}}_{1}-\hat{\phi}$ $\displaystyle=$
$\displaystyle 2(\hat{\hat{\phi}}_{1}-\hat{\phi})+\hat{\phi}$
$\displaystyle\sim$ $\displaystyle\mathcal{N}\left(\phi+3\mu,5\Sigma\right)$
$\displaystyle E_{\hat{\phi},\hat{\hat{\phi}}_{1}}\left[\Delta
S(2\hat{\hat{\phi}}_{1}-\hat{\phi})\right]$ $\displaystyle=$
$\displaystyle\Delta S(\phi)+3\mu\frac{\partial\Delta
S}{\partial\phi}|_{\phi}+\frac{9}{2}\mu\frac{\partial^{2}\Delta
S}{\partial\phi^{2}}|_{\phi}\mu+\frac{5}{2}\mathrm{tr}\left(\Sigma\frac{\partial^{2}\Delta
S}{\partial\phi^{2}}|_{\phi}\right)\,,$ (C2)
$\displaystyle\hat{\hat{\phi}}_{1}+\hat{\hat{\phi}}_{2}-\hat{\phi}$
$\displaystyle=$
$\displaystyle(\hat{\hat{\phi}}_{1}-\hat{\phi})+(\hat{\hat{\phi}}_{2}-\hat{\phi})+\hat{\phi}$
$\displaystyle\sim$ $\displaystyle\mathcal{N}\left(\phi+3\mu,3\Sigma\right)$
$\displaystyle
E_{\hat{\phi},\hat{\hat{\phi}}_{1},\hat{\hat{\phi}}_{2}}\left[\Delta
S(\hat{\hat{\phi}}_{1}+\hat{\hat{\phi}}_{2}-\hat{\phi})\right]$
$\displaystyle=$ $\displaystyle\Delta S(\phi)+3\mu\frac{\partial\Delta
S}{\partial\phi}|_{\phi}+\frac{9}{2}\mu\frac{\partial^{2}\Delta
S}{\partial\phi^{2}}|_{\phi}\mu+\frac{3}{2}\mathrm{tr}\left(\Sigma\frac{\partial^{2}\Delta
S}{\partial\phi^{2}}|_{\phi}\right)\,.$ (C3)

Inserting those relations to the LHS of Equation 15 easily leads to the RHS.

## References

  * Bailey (2021) Bailey, V. 2021, CGI-flux-ratio-plot.  https://github.com/nasavbailey/DI-flux-ratio-plot
  * Baudoz et al. (2005) Baudoz, P., Boccaletti, A., Baudrand, J., & Rouan, D. 2005, Proceedings of the International Astronomical Union, 1, 553
  * Beuzit et al. (2019) Beuzit, J.-L., Vigan, A., Mouillet, D., et al. 2019, Astronomy & Astrophysics, 631, A155, doi: 10.1051/0004-6361/201935251
  * Clare (2004) Clare, R. M. 2004, doi: 10.26021/1762
  * Currie et al. (2011) Currie, T., Burrows, A., Itoh, Y., et al. 2011, The Astrophysical Journal, 729, 128, doi: 10.1088/0004-637X/729/2/128
  * Dekking (2005) Dekking, M. 2005, A modern introduction to probability and statistics: understanding why and how (London: Springer).  http://site.ebrary.com/id/10228809
  * Fienup (1999) Fienup, J. R. 1999, Journal of the Optical Society of America A, 16, 1831, doi: 10.1364/JOSAA.16.001831
  * Fienup (2013) —. 2013, Applied optics, 52, 45
  * Frazin & Rodack (2021) Frazin, R. A., & Rodack, A. T. 2021, arXiv preprint arXiv:2105.06590
  * Gerchberg & Saxton (1972) Gerchberg, R. W., & Saxton, W. O. 1972, Optik, 35, 237
  * Groff et al. (2015) Groff, T. D., Kasdin, N. J., Limbach, M. A., et al. 2015, in Techniques and Instrumentation for Detection of Exoplanets VII, Vol. 9605 (SPIE), 457–466
  * Konopacky et al. (2013) Konopacky, Q. M., Barman, T. S., Macintosh, B. A., & Marois, C. 2013, Science, 339, 1398
  * Lafreniere et al. (2007) Lafreniere, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, . 2007, The Astrophysical Journal, 660, 770
  * Madec (2012) Madec, P.-Y. 2012, in Adaptive Optics Systems III, Vol. 8447 (SPIE), 22–39
  * Marois et al. (2006) Marois, C., Lafreniere, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, The Astrophysical Journal, 641, 556
  * Mawet et al. (2012) Mawet, D., Pueyo, L., Lawson, P., et al. 2012, in Space Telescopes and Instrumentation 2012: Optical, Infrared, and Millimeter Wave, Vol. 8442 (SPIE), 62–82
  * Neyman & Pearson (1933) Neyman, J., & Pearson, E. S. 1933, Philosophical Transactions of the Royal Society of London Series A, 231, 289, doi: 10.1098/rsta.1933.0009
  * Nielsen et al. (2019) Nielsen, E. L., De Rosa, R. J., Macintosh, B., et al. 2019, The Astronomical Journal, 158, 13
  * Nir et al. (2019) Nir, G., Zackay, B., & Ofek, E. O. 2019, The Astronomical Journal, 158, 70, doi: 10.3847/1538-3881/ab27c7
  * Noll (1976) Noll, R. J. 1976, JOsA, 66, 207
  * Platt & Shack (2001) Platt, B. C., & Shack, R. 2001, History and principles of Shack-Hartmann wavefront sensing, SLACK Incorporated Thorofare, NJ
  * Racine et al. (1999) Racine, R., Walker, G. A., Nadeau, D., Doyon, R., & Marois, C. 1999, Publications of the Astronomical Society of the Pacific, 111, 587
  * Rodack et al. (2021) Rodack, A. T., Frazin, R. A., Males, J. R., & Guyon, O. 2021, arXiv preprint arXiv:2105.06589
  * Ruane et al. (2018) Ruane, G., Riggs, A., Mazoyer, J., et al. 2018, in Space Telescopes and Instrumentation 2018: Optical, Infrared, and Millimeter Wave, Vol. 10698 (SPIE), 851–869
  * Ruffio et al. (2017) Ruffio, J.-B., Macintosh, B., Wang, J. J., et al. 2017, The Astrophysical Journal, 842, 14, doi: 10.3847/1538-4357/aa72dd
  * Shechtman et al. (2015) Shechtman, Y., Eldar, Y. C., Cohen, O., et al. 2015, IEEE signal processing magazine, 32, 87
  * Skaf et al. (2021) Skaf, N., Guyon, O., Boccaletti, A., et al. 2021, arXiv:2110.15361 [astro-ph].  http://arxiv.org/abs/2110.15361
  * Skaf et al. (2022) Skaf, N., Guyon, O., Gendron, ., et al. 2022, Astronomy & Astrophysics, 659, A170
  * Walter et al. (2020) Walter, A. B., Fruitwala, N., Steiger, S., et al. 2020, Publications of the Astronomical Society of the Pacific, 132, 125005, doi: 10.1088/1538-3873/abc60f


contiguous region $A$ and the (infinite) complement $\bar{A}$, as a function
of the distance $x$ between the qubit and the boundary of $A$. Numerical data
from stabilizer simulations of the spacetime-dual of a random Clifford circuit
(showing only $x<|A|/2$); annealed average over $2\times 10^{4}$-$5\times
10^{5}$ realizations depending on size. At $|A|=1024$ the best fit to a power
law decay is $\sim x^{-1.2}$ (dotted line), though there is a substantial
finite-size drift towards larger exponents. Also shown for comparison are
power-laws $x^{-1}$ (dashed line) and $x^{-1.5}$ (dot-dashed line).

To compute this quantity in our setting, we consider a finite subsystem $A$
embedded in an infinite system, initialized in a Bell pair (pure) product
state, and evolved under the dual of a random unitary circuit for a long time;
in the thermodynamic limit, tracing out $\bar{A}$ yields a cancellation
similar to the one we used for the pure-state entanglement calculation in Fig.
3, which elides all of the circuit except the portion contained inside a
“light cone” terminating at $A$. One is left with a finite-size calculation on
subsystem $A$ only, with a fully mixed state and “light-cone” boundary
conditions as in Ref. [30]. It is then possible to compute $f(x)$ from the
resulting mixed state on $A$ by exploiting the fact that $A\cup\bar{A}$ is in
a pure state and using

$\displaystyle f(x)$
$\displaystyle=I_{2}([x]:\bar{A})=S_{[x]}+S_{\bar{A}}-S_{[x]\cup\bar{A}}$
$\displaystyle=S_{[x]}+S_{A}-S_{A\setminus[x]}$ (13)

All three entropies in the expression are computable from $\rho_{A}$ without
reference to $\bar{A}$. Hence our calculation implicitly takes the
$\bar{A}\to\infty$ limit.

The numerical results for subsystems $A$ of up to 1024 qubits, shown in Fig.
12, suggest that the best fit to a power law $\sim x^{-\lambda}$ is
$\lambda\simeq 1.2$. This is in remarkable agreement with the results of Ref.
[98], which finds an exponent $\lambda\simeq 1.25$ in an almost identical
quantity, strongly suggesting that the two volume-law phases (in generic
monitored circuits and in duals of chaotic circuits) are in fact the same.

##  Appendix D Random-walk entanglement computation for Griffiths circuits

In this Appendix we present results on entanglement in spacetime-duals of the
“Griffiths circuits” discussed in Sec. V based on an alternative method, which
also allows us to analyze more carefully the role of corrections to scaling.

We consider the random-walk computation of the annealed average of the
entanglement entropy $S^{(a)}$ for Haar-random circuits, as discussed in Sec.
IV. The idea there is to take advantage of the edge decoherence picture and
adapt a method developed for entanglement growth in unitary circuits, where a
recursive formula (Eq. (4)) can be used to compute the purity exactly. In the
Griffiths circuit model, as long as the non-identity gates are sampled from
the Haar measure, the same approach goes through, up to a small modification
to the recursion: the purity $\mathcal{P}(x,t)$ is either updated as in Eq.
(4) (if a gate $U\neq\mathbb{1}$ acts at bond $x$ at time $t$), or left
unchanged (if no gate acts). These two occurrences are weighted according to
the rate $\gamma_{x}$ of nontrivial gates acting at bond $x$, giving

$\displaystyle\mathcal{P}(x,t)$
$\displaystyle=\gamma_{x}\frac{q^{-v_{E}}}{2}[\mathcal{P}(x-1,t-1)+\mathcal{P}(x+1,t-1)]$
$\displaystyle\qquad+(1-\gamma_{x})\mathcal{P}(x,t-1)$ (14)

away from the decohering edge, and an analogous modification for the recursion
at the edge.

We remark that evaluating this recursive formula yields the average of
$\mathcal{P}$ over both the _location_ (whether a gate $U\neq\mathbb{1}$ is
present at bond $x$ at time $t$) and _choice_ (which $U\in{\rm U}(q^{2})$) of
Haar-random gates, given the (quenched) rates $\gamma_{x}$. The average
$\overline{\mathcal{P}}$ obtained in this way then gives an entropy,
$S=-\log_{q}\overline{\mathcal{P}}$, which can then be averaged over the
choice of rates $\\{\gamma_{x}\\}$ out of the distribution
$P(\gamma)=(a+1)\gamma^{a}$. The resulting average of $S$ is thus partly
annealed and partly quenched222222One could also average
$\overline{\mathcal{P}}$ over the choice of rates, getting the overall
annealed average of $S$, but this is statistically problematic in the model at
hand, where fluctuations of purity over different choices of rates are large..

Figure 13: (Top) Entropy $\tilde{S}$ vs subsystem size $\tilde{\ell}$ for the
Griffiths circuit model with $a=1$, with Haar-random gates, at depth
$\tilde{T}=1024$ (sufficient for saturation to $\tilde{S}_{\infty}$ at the
sizes shown). The data is obtained from the recursive method, for qubits
($q=2$) and 10-dimensional qudits ($q=10$). (Bottom) Logarithmic derivative
reveals a drift of the data towards the predicted value. The corrections to
scaling are large for qubits, but get much smaller at large $q$.

Results obtained by numerically evaluating the recursion formula on
$\tilde{T}=1024$ qubits evolving under edge decoherence are shown in Fig. 13,
for qubits (the relevant case to compare to the Clifford numerics of Fig. 8)
as well as systems with large local Hilbert spaces ($q=10$), compared to the
analytical prediction $\tilde{S}_{\infty}\sim\tilde{\ell}^{\alpha}$,
$\alpha=(a+1)/(a+2)$. A careful look at the size dependence of $\tilde{S}$,
e.g. via the logarithmic derivative $d\log(\tilde{S})/d\log(\tilde{\ell})$,
reveals that the asymptotic scaling is consistent with the predicted one, but
that corrections to scaling are substantial for qubits, so that the best fit
for the exponent in $\tilde{S}\sim\tilde{\ell}^{\alpha}$ is significantly off
even at large sizes, $\tilde{\ell}\gtrsim 10^{3}$. (The data shown is for
$a=1$, where the predicted scaling is $\alpha=2/3$ and the discrepancy with
the finite-size data is worst.) This suggests that the data in Fig. 8 (from
stabilizer simulations) is consistent with the analytical prediction. A more
thorough investigation of subleading corrections, aiming to decouple the
effects of (i) unitary dynamics vs edge decoherence, (ii) average over the
Haar measure vs the Clifford group, and (iii) quenched vs annealed vs mixed
averages, is an interesting question for future studies. In particular, it
would be exciting to detect the logarithmic correction predicted for this
model, though the presence of other power-law corrections due to randomness
likely makes this task very hard.