## 1.

## Introduction

The quest to characterize Earth-like planets was brought into sharp focus by the discovery of a terrestrial exoplanet orbiting the nearest star, Proxima Centauri.^{1} With a maximum projected separation of around 36 mas, Proxima b could in principle be characterized with current generation instruments,^{2} and will be readily studied with next-generation giant segmented mirror telescopes (GSMTs) such as the 25-m Giant Magelllan Telescope (GMT), 30-m Thirty Meter Telescope, and 39-m Extremely Large Telescope. Indeed, such observations are recognized as one of the key science goals of the GSMTs,^{3}^{,}^{4} and Proxima b will not be the only target. Results from the Kepler mission^{5} have made it clear that small planets occur frequently around low-mass stars.^{6}^{,}^{7} In addition to Proxima b, the recent discoveries of Ross 128 b,^{8} LHS 1140 b,^{9} and the seven terrestrial planets orbiting TRAPPIST-1^{10} likewise point to a high frequency of terrestrial planets in low-mass star systems near the Earth.

While some nearby planets, such as those orbiting TRAPPIST-1 and LHS 1140 b, will be amenable to atmospheric characterization with the transit spectroscopy technique, the low *a priori* transit probability means that the majority of nearby planets will have to be characterized with resolved direct imaging. Low-mass stars present unique challenges and opportunities. With lower stellar luminosity, planets with the same equilibrium temperature occur closer to the star. This results in a significant improvement in contrast between the planet and the star, but also results in needing large diameter telescopes to resolve planets at such small separations.^{4} This places such observations in the realm of ground-based telescopes equipped with adaptive optics (AO) and coronagraphs.

The analysis presented here is motivated by the question: what are the fundamental limits of ground-based high-contrast coronagraphic imaging of terrestrial exoplanets? Here we attack one part of this complicated issue: finding the fundamental limit imposed on coronagraphic contrast by residual atmospheric turbulence behind a closed-loop AO control system.

Much previous work has been done modeling AO performance in the spatial-frequency domain.^{11}^{–}^{14} This technique treats the AO system as a spatial filter, casting time and temporal-frequency domain processes (such as measurement noise and control filtering) into the the spatial domain as filters, and analyzing the effect of the combined AO system on the input turbulence power spectral density (PSD).

The previous analytic treatment of AO performance most closely related to our treatment here is Guyon.^{15} Using the correspondence between Fourier modes and speckle amplitude, coupled with straightforward scaling from the input PSD, a derivation of the postcoronagraph contrast and its dependence on various AO system parameters was developed. The resulting scaling laws are readily applicable to system design and analysis and have proven to be useful in placing results from RV and transit surveys in context.^{16}

The goal of this study is to present a way to analyze the achievable contrast in a closed-loop AO system, including both current generation telescopes and future GSMTs. To do this we develop a semianalytic framework in the temporal frequency domain. That is, we explicitly calculate temporal PSDs rather than using the spatial-filtering approach. At least in the frozen-flow regime, these should produce nearly identical results. Our exploration of this method was motivated by the fact that actual AO systems must be optimized and controlled in the time and temporal-frequency domains.

The paper is organized as follows. In Sec. 2, we state our entering assumptions. In Sec. 3, we lay out our notation and physical description of the problem, including the statistics of the Fourier basis. In Sec. 4, we derive a coronagraph model and deal with a subtle detail depending on how residual variance is calculated. In Sec. 5, we calculate the temporal PSDs of Fourier modes assuming von Kármán turbulence in frozen flow. In Sec. 6, we review the standard model of AO control with optimal gains and present a predictive controller based on the linear prediction formalism. Then, in Sec. 7, we use the previous results to predict the postcoronagraph contrasts for 6.5- and 25.4-m telescopes for a range of guide star magnitudes and compare optimal integrator and linear predictive controllers. We discuss our results in Sec. 8 and conclude in Sec. 9.

## 2.

## Assumptions

This analysis is concerned with AO in the “extreme” regime: high actuator counts ($\sim 2000$ on a 6.5-m aperture, $\sim \mathrm{21,000}$ on a 25.4 m) and high loop speeds ($\sim 2000\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{Hz}$ or more). Furthermore, we are mainly interested in the coronagraphic raw contrast, that is the ratio of the intensity of the point spread function (PSF) at some separation from the star, to the intensity at the peak of the star’s image or the PSF. Such extreme-AO (ExAO) systems are typically constructed for the purpose of high-contrast imaging of exoplanets and other narrow-angle circumstellar targets.

We do not seek to describe any particular ExAO system, nor address the details of design. Rather, our goal is to assess what can be achieved with an ideal system limited only by the reality of atmospheric turbulence. As such, our main assumptions are as follows:

1. We consider only on-axis guide stars and neglect anisoplanatism. This is of minor concern at the small angular separations we are concerned with for exoplanet imaging.

2. Our analysis will be monochromatic. We do not consider the effects of atmospheric dispersion or chromaticity of the air index of refraction or scintillation within the wavefront sensor (WFS) band. We also assume that the WFS and science wavelengths are the same, ignoring differential chromatic effects.

3. We ignore WFS aliasing. This is an approximation, but with justification: we analyze a pyramid WFS (PyWFS), which is less susceptible to aliasing than the Shack–Hartmann WFS.

^{17}Furthermore, a spatial filter^{18}could be incorporated in a PyWFS essentially eliminating aliasing from concern.4. We assume perfect knowledge of the AO system components. This means we neglect errors in calibration (e.g., WFS gain) and assume that we know the spatial and temporal transfer functions of all components.

5. We will ignore all non-Kolmogorov wavefront error sources. This means we do not consider static or noncommon path errors within the coronagraph nor analyze telescope error sources (e.g., vibrations). We also consider only free atmosphere turbulence, that is, we do not consider dome seeing. Such error sources are ultimately important in real systems; however, free atmosphere turbulence is the dominant term which must be dealt with first.

6. We will not consider changes in turbulence, e.g., changes in the wind velocity or changes in seeing.

7. We will not consider details of coronagraphic architecture. However, we do show that our idealized model is representative of real-world coronagraphs.

Thus, what follows can be taken as goal setting: we desire to know how well a coronagraph could perform if limited only by atmospheric turbulence.

## 3.

## Incoming Wavefront

We begin by introducing our notation and physical description of the problem. We summarize our notation in Table 1.

## Table 1

Summary of notation.

Symbol | Units | Parameter |
---|---|---|

$\overrightarrow{q}$ | m | Position vector in the pupil plane |

$u,v$ | m | Position coordinates in the pupil plane |

$\stackrel{^,}{u}\widehat{v}$ | Coordinate unit vectors in the pupil plane | |

$\overrightarrow{k}$ | ${\mathrm{m}}^{-1}$ | Spatial frequency vector in the pupil plane |

$\mathrm{\Delta}k$ | ${\mathrm{m}}^{-1}$ | Discrete spatial frequency sampling |

${k}_{u},{k}_{v}$ | ${\mathrm{m}}^{-1}$ | Spatial frequency components in the pupil plane |

$m,n$ | Integer indices of discrete spatial frequencies | |

$D$ | m | Telescope diameter |

$\mathcal{A}(\overrightarrow{q})$ | Aperture function | |

$\lambda $ | m | Wavelength of observation and wavefront sensing |

$\mathcal{P}$ | ${\mathrm{rad}}^{2}\text{\hspace{0.17em}}{\mathrm{m}}^{2}$ | Spatial power spectrum of the wavefront |

${\lambda}_{0}$ | m | Reference wavelength for turbulence parameters |

${r}_{0}$ | m | Fried parameter for atmospheric turbulence strength |

$Q$ | Denotes the Fourier transform of a mode | |

$z$ | m | Height above the observatory |

${C}_{n}^{2}({z}_{i})$ | Normalized turbulence strength for layer $i$ | |

${M}_{mn}$ | Denotes a Fourier mode | |

${h}_{mn}$ | m | Amplitude of a mode describing the wavefront phase |

${a}_{mn}$ | Amplitude of a mode describing the wavefront amplitude | |

$\mathrm{\varphi}(\overrightarrow{q},t)$ | rad | The phase of the wavefront |

$A(\overrightarrow{q},t)$ | The amplitude of the wavefront | |

${\sigma}_{mn}^{2}$ | ${\mathrm{rad}}^{2}$ | The variance of the amplitude of a mode |

$\mathcal{V}$ | $\mathrm{m}/\mathrm{s}$ | Wind speed |

$\mathrm{\Theta}$ | rad | Wind direction |

$\tau $ | s | Denotes a time delay |

$\overrightarrow{r}$ | $\lambda /D$ | Position vector in the image plane |

$I$ | $\text{photons}/\mathrm{s}/{(\lambda /D)}^{2}$ | Intensity in the image plane |

PSF | $\text{photons}/\mathrm{s}/{(\lambda /D)}^{2}$ | Point-spread function |

$S$ | Strehl ratio | |

$C$ | Contrast ratio | |

$f$ | Hz | Temporal frequency |

${\mathcal{T}}_{mn}$ | ${\mathrm{rad}}^{2}\text{\hspace{0.17em}\hspace{0.17em}}/\mathrm{Hz}$ | Temporal power spectrum of a single mode |

$H(s)$ | Denotes a transfer function |

## 3.1.

### Coordinate System

We will describe the incoming wavefront at the entrance pupil of the telescope. The position vector $\overrightarrow{q}$ in this plane is defined by the unit vectors $(\widehat{u},\widehat{v})$, such that

We will restrict our analysis to the unobstructed circular aperture of diameter $D$ defined by

## Eq. (2)

$$\mathcal{A}(\overrightarrow{q})=\{\begin{array}{cc}\frac{4}{\pi {D}^{2}},& \text{if}\text{\hspace{0.17em}\hspace{0.17em}}q\le \frac{D}{2}\\ 0,& \text{if}\text{\hspace{0.17em}\hspace{0.17em}}q>\frac{D}{2}.\end{array}$$Spatial-frequency is a vector quantity

Wavefront control is inherently a discrete problem, and the aperture defines the discrete spatial-frequency sampling

Hence, we will often discretize spatial-frequency as

where $m$ and $n$ are integer indices.## 3.2.

### Spatial Power Spectral Densities

The spatial PSD of the phase at the observation wavelength $\lambda $ in the von Kármán model^{19} is

## Eq. (6)

$${\mathcal{P}}_{\mathrm{vK}}(\overrightarrow{k},\lambda )=0.0218{\left(\frac{{\lambda}_{0}}{\lambda}\right)}^{2}\frac{1}{{r}_{0}^{5/3}}\frac{1}{{({k}^{2}+{k}_{0}^{2})}^{11/6}}\text{\hspace{0.17em}\hspace{0.17em}}({\mathrm{rad}}^{2}/{\mathrm{m}}^{-2}),$$^{20}and ${\lambda}_{0}$ is the wavelength at which ${r}_{0}$ is reported. The outer scale, ${L}_{0}$, is included in the parameter

The mean value over the aperture, or the piston term, has no effect on our problem. However, the amplitude of the piston mode is nonzero in the above PSD, so we should subtract the power in piston from the PSD. The spatial PSD of a mode on the aperture is given as^{21}^{,}^{22}

## Eq. (8)

$${\mathcal{P}}_{\text{mode}}(\overrightarrow{k})={\mathcal{P}}_{\mathrm{o}}(\overrightarrow{k}){|{Q}_{\text{mode}}(\overrightarrow{k})|}^{2},$$^{23}

So the piston subtracted phase PSD is

## Eq. (10)

$${\mathcal{P}}_{\mathrm{vK},\mathrm{sub}}(\overrightarrow{k},\lambda )={\mathcal{P}}_{\mathrm{vK}}(\overrightarrow{k},\lambda )\{1-{[2\mathrm{Ji}(\pi Dk)]}^{2}\},$$We model atmospheric turbulence as occurring in discrete layers. To account for the Fresnel propagation among these layers, we use the functions^{15}

## Eq. (12)

$$X(k,\lambda )=\sum _{i}{C}_{n}^{2}({z}_{i}){\mathrm{cos}}^{2}(\pi {z}_{i}{k}^{2}\lambda ),\phantom{\rule{0ex}{0ex}}Y(k,\lambda )=\sum _{i}{C}_{n}^{2}({z}_{i}){\mathrm{sin}}^{2}(\pi {z}_{i}{k}^{2}\lambda ),$$## 3.3.

### Representing the Wavefront

## 3.3.1.

#### Fourier basis

Here, we develop a real-valued form of the Fourier basis. Consider the complex exponential form of the basis^{24}

## Eq. (16)

$${M}_{mn}=\mathrm{Re}\{h\}[\mathrm{cos}(2\pi {\overrightarrow{k}}_{mn}\xb7\overrightarrow{q})+i\text{\hspace{0.17em}}\mathrm{sin}(2\pi {\overrightarrow{k}}_{mn}\xb7\overrightarrow{q})]+\mathrm{Im}\{h\}[i\text{\hspace{0.17em}}\mathrm{cos}(2\pi {\overrightarrow{k}}_{mn}\xb7\overrightarrow{q})-\mathrm{sin}(2\pi {\overrightarrow{k}}_{mn}\xb7\overrightarrow{q})],$$## Eq. (17)

$${M}_{mn}^{p}(\overrightarrow{q})=\mathrm{cos}(2\pi {\overrightarrow{k}}_{mn}\xb7\overrightarrow{q})+p\text{\hspace{0.17em}}\mathrm{sin}(2\pi {\overrightarrow{k}}_{mn}\xb7\overrightarrow{q}),$$## Eq. (18)

$$\int \mathcal{A}(\overrightarrow{q}){[{M}_{mn}^{p}(\overrightarrow{q})]}^{2}\mathrm{d}\overrightarrow{q}=1,$$## Eq. (19)

$${Q}_{mn}^{p}(\overrightarrow{k})=(1+{i}^{p})Ji(\pi D{k}_{mn}^{+})+(1-{i}^{p})Ji(\pi D{k}_{mn}^{-}),$$## Eq. (20)

$$\begin{array}{c}{k}_{mn}^{+}=\sqrt{{({k}_{u}+\frac{m}{D})}^{2}+{({k}_{v}+\frac{n}{D})}^{2}},\\ {k}_{mn}^{-}=\sqrt{{({k}_{u}-\frac{m}{D})}^{2}+{({k}_{v}-\frac{n}{D})}^{2}}.\end{array}$$We have chosen this form of the Fourier basis primarily because it yields the useful result

This means that each spatial frequency can be described by a single spatial PSD, and, as we will show a single one-sided (positive frequencies only) temporal PSD.

## 3.3.2.

#### The Wavefront in the Fourier basis

We denote the wavefront phase as $\mathrm{\varphi}(\overrightarrow{q},t)$, which has units of radians. The coefficient of a mode is calculated as

## Eq. (22)

$${h}_{mn}^{p}(t)=\frac{\lambda}{2\pi}\int \mathcal{A}(\overrightarrow{q}){M}_{mn}^{p}(\overrightarrow{q})\mathrm{\varphi}(\overrightarrow{q},t)\mathrm{d}\overrightarrow{q}.$$These coefficients are real-valued. Poyneer and Veran^{25} showed that the Fourier basis can be used for analysis and synthesis, even though it is not orthogonal on an aperture. Since our slightly modified basis is simply a linear combination of the cosine and sine, their results apply to it as well. This means that we can write

## Eq. (23)

$$\mathrm{\varphi}(\overrightarrow{q},t)=\frac{2\pi}{\lambda}\sum _{mn}[{h}_{mn}^{+}(t){M}_{mn}^{+}(\overrightarrow{q})+{h}_{mn}^{-}(t){M}_{mn}^{-}(\overrightarrow{q})].$$Similarly, we represent the wavefront amplitude aberrations as

## 3.4.

### Statistics of the Fourier Basis

Following the recipe of Noll,^{23} we can write the covariance of any two modes over the aperture as

## Eq. (25)

$$\u27e8{h}_{mn}^{p*}{h}_{{m}^{\prime}{n}^{\prime}}^{{p}^{\prime}}\u27e9={\int}_{0}^{\infty}{\int}_{0}^{2\pi}{Q}_{mn}^{p*}{Q}_{{m}^{\prime}{n}^{\prime}}^{{p}^{\prime}}\mathcal{P}(k)\mathrm{d}\theta k\text{\hspace{0.17em}}\mathrm{d}k.$$We now make use of this expression to analyze the statistics of the Fourier basis under von Kármán turbulence.

## 3.4.1.

#### Ignoring the aperture

In Eq. (25), the Fourier transform $Q$ describes how the PSD is sampled by the modes. If we were to ignore the spatial effect of the aperture, we would have

where $\delta (x)$ is the Dirac delta distribution, which is the Fourier transform of the sine and cosine on an infinite domain. With this approximation, we have the variance of a mode## 3.4.2.

#### Modal variance

Including the aperture, the variance of a mode in radians is given as

## Eq. (28)

$${\sigma}_{mn}^{2}={\left(\frac{2\pi}{\lambda}\right)}^{2}\u27e8{|{h}_{mn}^{p}|}^{2}\u27e9={\int}_{0}^{\infty}{\int}_{0}^{2\pi}{|{Q}_{mn}^{p}|}^{2}\mathcal{P}(k)\mathrm{d}\theta k\text{\hspace{0.17em}}\mathrm{d}k.$$Now since

## Eq. (29)

$${|{Q}_{mn}^{p}|}^{2}={Q}_{mn}^{p*}{Q}_{mn}^{p}=2[{\mathrm{Ji}}^{2}(\pi D{k}_{mn}^{+})+{\mathrm{Ji}}^{2}(\pi D{k}_{mn}^{-})],$$In Fig. 1, we show the modal variance versus spatial frequency. This was calculated using numerical integration in `double` precision with the GNU Scientific Library (GSL^{26}) routines `gsl_integration_qagiu` (for the radial integral) and `gsl_integration_qag` (for the azimuthal integral). We found that an absolute tolerance of ${10}^{-10}$ and a relative tolerance of ${10}^{-4}$ gave good results.

Figure 1 also shows the naive result using Eq. (27). The simple estimate based on the PSD is significantly below the true variance.

One might guess that AO correction will reduce or eliminate this discrepancy. A simple way to estimate the effect of AO control on these statistics is to assume frozen flow and apply a scaling for the time-delay error. Following Guyon^{15} (as modified in Appendix B), the corrected PSD is approximated as

## Eq. (30)

$${\mathcal{P}}_{\mathrm{cor}}(\overrightarrow{k},\lambda )\approx \{\begin{array}{cc}{\mathcal{P}}_{\mathrm{\varphi}}(\overrightarrow{k}){(2\pi \mathcal{V}k)}^{2}{\tau}_{\mathrm{tl}}^{2},& k\le \frac{D}{2d}\\ {\mathcal{P}}_{\mathrm{\varphi}}(\overrightarrow{k}),& k>\frac{D}{2d}\end{array},$$The variance with and without the aperture is also shown in Fig. 1. Though the discrepancy is reduced, the simple estimate based on the PSD is still below the true variance. This means that the variance of a Fourier mode is not simply proportional to the PSD when it is sampled by the aperture.

The discrepancy is due to the aperture: there is a convolution inherent in Eq. (25) not addressed by Eq. (27). We show how to address this and discuss its implications for various coronagraph models, in Sec. 4.

## 3.4.3.

#### Correlation

We next calculated the covariance of pairs of Fourier modes according to Eq. (25). The results are shown in Fig. 2 for uncorrected von Kármán turbulence, and in Fig. 3 for the corrected case. In the figure, we show the covariances normalized by the modal variances, that is the correlation. There is significant correlation between various modes in the uncorrected case, but an AO system acts to suppress this. In the following section we show that this correlation between Fourier modes can be accounted for by convolving the PSF with the PSD.

## 4.

## Postcoronagraph Contrast

We next derive a PSD-based model of a coronagraph fed by an AO system in the long-exposure limit and assess its validity.

## 4.1.

### Coronagraph Model

The intensity in an image plane is given as

## Eq. (31)

$$I(\overrightarrow{r},t)={|\mathcal{F}\{\mathcal{A}(\overrightarrow{q})A(\overrightarrow{q},t){e}^{i\mathrm{\varphi}(\overrightarrow{q},t)}\}|}^{2}.$$The operator $\mathcal{F}\{\xb7\}$ describes the propagation of the wavefront from pupil to final image plane. From here on we assume this is the Fourier transform defined as in Appendix A.2, though it could be a series of transforms corresponding to various coronagraph components.^{27} We define the system PSF as the result of the above with no aberrations, that is ${h}_{mn}=0$ and ${a}_{mn}=0$ for all $m,n$, which is

## Eq. (32)

$$\mathrm{PSF}(\overrightarrow{r})={|\mathcal{F}\{\mathcal{A}(\overrightarrow{q})\}|}^{2}.$$In the case of the circular unobstructed aperture [Eq. (2)], this is the airy pattern

where we have used the image plane coordinate $\overrightarrow{r}=\lambda \overrightarrow{k}$. Telescope apertures are typically more complicated, consisting of secondary mirror obscuration and support structures, possibly segments, and can be noncircular. Furthermore, it could be position dependent due to the coronagraph. The result will be more difficult to analyze, usually requiring numerical calculation.Equation (31) can not be evaluated in closed form in the general case. However, by assuming small aberrations and an ideal coronagraph, we can use a second-order Taylor expansion^{28}^{,}^{29} to approximate the postcoronagraph PSF as^{30}

## Eq. (34)

$$I(\overrightarrow{r},t)\approx {\left|\mathcal{F}\right\{\mathcal{A}(\overrightarrow{q})\sum _{mn}{A}_{mn}(\overrightarrow{q},t)\left\}\right|}^{2}+{|\mathcal{F}\{\mathcal{A}(\overrightarrow{q})\mathrm{\varphi}(\overrightarrow{q},t)\}|}^{2}.$$Macintosh et al.^{30} only addressed phase, but the extension to amplitude is straightforward from their arguments. Dealing with the phase component first, we have

## Eq. (35)

$${I}_{\mathrm{\varphi}}(\overrightarrow{r},t)\approx {|\sum _{mn}\mathcal{F}\{\mathcal{A}(\overrightarrow{q}){\mathrm{\varphi}}_{mn}(\overrightarrow{q},t)\}|}^{2}.$$Evaluating a single term gives

## Eq. (36)

$$\mathcal{F}\{\mathcal{A}(\overrightarrow{q}){\mathrm{\varphi}}_{mn}(\overrightarrow{q},t)\}=\frac{2\pi}{\lambda}[{h}_{mn}^{+}(t){Q}_{mn}^{+}(\overrightarrow{k})+{h}_{mn}^{-}(t){Q}_{mn}^{-}(\overrightarrow{k})].$$Then taking the modulus squared and pulling out the terms with $mn={m}^{\prime}{n}^{\prime}$ gives (we are suppressing the $t$ dependence of $h$)

## Eq. (37)

$${I}_{\mathrm{\varphi}}(\overrightarrow{r},t)\approx {\left(\frac{2\pi}{\lambda}\right)}^{2}2\sum _{mn}\{[{({h}_{mn}^{+})}^{2}+{({h}_{mn}^{-})}^{2}][{\mathrm{Ji}}^{2}(\pi D{k}_{mn}^{+})+{\mathrm{Ji}}^{2}(\pi D{k}_{mn}^{-})]+8{h}_{mn}^{+}{h}_{mn}^{-}\mathrm{Ji}(\pi D{k}_{mn}^{+})\mathrm{Ji}(\pi D{k}_{mn}^{-})\phantom{\rule{0ex}{0ex}}+\sum _{\genfrac{}{}{0ex}{}{{m}^{\prime}{n}^{\prime}}{mn\ne {m}^{\prime}{n}^{\prime}}}\{[{h}_{mn}^{+}{h}_{{m}^{\prime}{n}^{\prime}}^{+}+{h}_{mn}^{-}{h}_{{m}^{\prime}{n}^{\prime}}^{-}][\mathrm{Ji}(\pi D{k}_{mn}^{+})\mathrm{Ji}(\pi D{k}_{{m}^{\prime}{n}^{\prime}}^{+})+\mathrm{Ji}(\pi D{k}_{mn}^{-})\mathrm{Ji}(\pi D{k}_{{m}^{\prime}{n}^{\prime}}^{-})]+[{h}_{mn}^{+}{h}_{{m}^{\prime}{n}^{\prime}}^{-}+{h}_{mn}^{-}{h}_{{m}^{\prime}{n}^{\prime}}^{+}][\mathrm{Ji}(\pi D{k}_{mn}^{+})\mathrm{Ji}(\pi D{k}_{{m}^{\prime}{n}^{\prime}}^{-})+\mathrm{Ji}(\pi D{k}_{mn}^{-})\mathrm{Ji}(\pi D{k}_{{m}^{\prime}{n}^{\prime}}^{+})]\}\},$$## Eq. (38)

$$\u27e8{I}_{\mathrm{\varphi},mn}(\overrightarrow{r})\u27e9=\frac{{\mathcal{P}}_{\mathrm{\varphi}}({\overrightarrow{k}}_{mn},\lambda )}{{D}^{2}}[\mathrm{PSF}(\overrightarrow{r}-{\overrightarrow{k}}_{mn}\lambda )+\mathrm{PSF}(\overrightarrow{r}+{\overrightarrow{k}}_{mn}\lambda )]$$^{31}We illustrate this in Fig. 4.

Finally, we perform the sum over all spatial frequencies. The complete intensity due to $\mathrm{\varphi}$ at $\overrightarrow{r}$ is

## Eq. (39)

$$\u27e8{I}_{\mathrm{\varphi}}(\overrightarrow{r})\u27e9=\sum _{mn}\u27e8{I}_{\mathrm{\varphi},mn}(\overrightarrow{r})\u27e9.$$This summation is an unnormalized convolution with the PSF and can be thought of as adding photons from nearby speckles. This step accounts for the effect of the aperture. In the AO corrected regime, it increases the intensity by as much as 50%, explaining the discrepancies in Fig. 1.

Now we apply the lesson of Sec. 3.4.2 above. If we instead use the variance calculated with Eq. (28), then we can skip the convolution and arrive directly at

## Eq. (40)

$$\u27e8{I}_{\mathrm{\varphi}}(\overrightarrow{r})\u27e9={\left(\frac{2\pi}{\lambda}\right)}^{2}\u27e8{|{h}_{mn}|}^{2}\u27e9[\mathrm{PSF}(\overrightarrow{r}-{\overrightarrow{k}}_{mn}\lambda )+\mathrm{PSF}(\overrightarrow{r}+{\overrightarrow{k}}_{mn}\lambda )].$$This is an assertion that the the complete evaluation of Eq. (37) in the long-exposure limit is equivalent to the modal variance calculation on the aperture with Eq. (28). Acknowledging that we have not actually proven this, we demonstrate its validity with numerical experiments in the following section.

The amplitude component of the postcoronagraph long-exposure intensity, $\u27e8{I}_{A}(\overrightarrow{r})\u27e9$, is evaluated in nearly identical fashion.

We define the “raw PSF contrast” as the ratio of the intensity at $\overrightarrow{r}$ to the peak intensity of a noncoronagraphic long-exposure image of the star, $\u27e8{I}_{\mathrm{NC}}(0)\u27e9$. This is related to the PSF by the long-exposure Strehl ratio $\u27e8S\u27e9$, defined by the relationship

The contrast is then

## 4.2.

### Numerical Verification

We conducted a set of numerical experiments to verify our long-exposure coronagraph model, compare it to other models, and to assess its relevance to real-world coronagraphy. Five distinct model coronagraphs were considered, in each case both for phase aberrations only and with the amplitude errors included.

The first comparison coronagraph model was the so-called “perfect coronagraph.”^{32} In this model, the complex plane wave that minimizes the energy in the pupil is subtracted, and the result is propagated to the image plane. It has been shown that analysis with the perfect coronagraph produces results that match the real-world four quadrant phase mask coronagraph.^{33}^{,}^{34}

The second coronagraph model consisted of propagating the complex wavefronts using Eq. (34) directly, employing the FFT. The phase and amplitude components were treated separately, then combined after taking the modulus squared.

The third coronagraph model corresponds to the measurement error and time delay limited contrast ${C}_{2}$ of Guyon,^{15} combined with the contrast from uncorrected amplitude ${C}_{1}$ due to atmospheric scintillation. In Appendix B, we present a derivation of ${C}_{2}$ utilizing our present notation and generalizing it for arbitrary parameters. This method uses the residual PSD, as in Eq. (27), and so requires a convolution. Hence, this tests the convolution solution for the discrepancy between the true variance and the simple PSD estimate highlighted above.

The fourth model consisted of directly measuring the variance of each modal coefficient. The coefficient of each mode was measured by inner product with the mean-subtracted phase and amplitude on the aperture for each randomly generated screen. In this case, the aperture is already taken into account so Eq. (40) describes the results.

The final comparison coronagraph was an apodized pupil Lyot complex mask coronagraph (APLCMC^{35}). We included this as an example of a coronagraph that can be implemented in the real world. In this case, the focal plane mask was 1.29 $\lambda /D$ in diameter at 800 nm, with a complex transmission of $-0.458$. The optimum pupil apodization gave a net intensity throughput of 0.587.

We generated a series of phase screens from the von Kármán PSD using Fourier-space convolution of the PSD with Gaussian white noise fields. Scintillation was included by multiplying the PSD by the $X$ and $Y$ functions, Eq. (12), resulting in phase and amplitude screens. AO correction of the phase was approximated as in Eq. (30). We used a Fried parameter of ${r}_{0}=0.2\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$, an outer scale ${L}_{0}=25\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$, and windspeed of $V=10\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}/\mathrm{s}$. The total time delay was set to 2.5 frames at 2000 Hz (1.25 ms). We set telescope diameter to $D=6.5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ and actuator spacing to $d=13.5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{cm}$ (48 across) and assumed a wavelength of $\lambda =800\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$. The phase screens were 2048 pixels across, and we used a circular unobstructed aperture 128 pixels across. We assumed an arbitrarily bright star so WFS measurement noise was ignored. The Strehl ratio of a noncoronagraphic image with these parameters was $\sim 94\%$.

10,000 realizations of the AO-corrected phase and uncorrected amplitude were generated. For each realization, we applied the five comparison coronagraph models. The results for these coronagraphs are shown in Fig. 5, where we see that, in general, these various coronagraph models match the APLCMC. We note that the “perfect coronagraph” and direct propagation using Eq. (34) are nearly identical to each other. The residual-PSD+convolution analysis also matches well, as does the direct calculation of modal variance.

## 4.3.

### Long-Exposure Contrast and the PSD

We can summarize the results of this section with the following points:

• The postcoronagraph contrast is not simply the residual PSD, due to the effect of the aperture. This is easily addressed by an unnormalized convolution with the PSF.

• If the variance is calculated with the aperture included, then the convolution is not needed.

• Ideal coronagraph models are valid in the regime of ground-based, well-corrected, AO-fed coronagraphy, and are useful proxies for real-world coronagraphs.

Our result that the postcoronagraph long-exposure image is a convolution of the PSF with the residual PSD is nearly identical to that developed by Herscovici-Schiller et al.^{27} It is fully equivalent if one uses a position-dependent PSF and includes internal aberrations, both details we neglect since we are analyzing an ideal system. The difference between this case, and the case when the aperture is included, is important when the modal variances are individually calculated as we do in the following sections.

## 5.

## Temporal Power Spectra

Next, we describe a process for calculating the temporal PSDs of the Fourier basis in wind-driven von Kármán turbulence.

## 5.1.

### Temporal PSD of Fourier Modes in Wind-Driven Turbulence

Restating Eq. (8), the spatial PSD of a Fourier mode on an unobstructed circular aperture is given as (as noted above, we are suppressing $p$)

## Eq. (43)

$${\mathcal{P}}_{mn}(\overrightarrow{k})={\mathcal{P}}_{o}(\overrightarrow{k}){|{Q}_{mn}(\overrightarrow{k})|}^{2}.$$We next invoke the Taylor or frozen-flow hypothesis, which for our purposes states that we can treat the time-domain behavior of turbulence as if fixed turbulent phase screens are blowing across the telescope in discrete layers. With this assumption we can calculate the temporal PSD of a turbulent mode with^{21}^{,}^{22}

## Eq. (44)

$${\mathcal{T}}_{mn}(f;\mathcal{V})=\frac{2}{\mathcal{V}}{\int}_{-\infty}^{\infty}{\mathcal{P}}_{mn}(\frac{f}{\mathcal{V}},{k}_{v})\mathrm{d}{k}_{v},$$Now to account for an arbitrary wind direction, let $\mathrm{\Theta}$ be the direction of the wind with respect to the $\widehat{u}$ axis. We then derotate the basis functions by this angle to align with the axes of the integral in Eq. (44) by defining

## Eq. (45)

$${m}^{\prime}=m\text{\hspace{0.17em}}\mathrm{cos}(\mathrm{\Theta})-n\text{\hspace{0.17em}}\mathrm{sin}(\mathrm{\Theta})\phantom{\rule{0ex}{0ex}}{n}^{\prime}=m\text{\hspace{0.17em}}\mathrm{sin}(\mathrm{\Theta})+n\text{\hspace{0.17em}}\mathrm{cos}(\mathrm{\Theta}).$$The temporal PSD for arbitrary wind direction is then simply

## Eq. (46)

$${\mathcal{T}}_{mn}(f;\mathcal{V},\mathrm{\Theta})=\frac{4}{\mathcal{V}}{\int}_{-\infty}^{\infty}{\mathcal{P}}_{{m}^{\prime}{n}^{\prime}}(\frac{f}{\mathcal{V}},{k}_{v})\mathrm{d}{k}_{v}.$$Finally, the temporal PSD of a Fourier mode for a turbulence profile is the ${C}_{n}^{2}$ weighted sum of the single layer profiles

## 5.2.

### Numerical Calculation

The improper integral in Eq. (46) does not in general have a simple closed-form solution, so we must calculate it numerically. In this work, we used the GSL `gsl_integration_qagi`^{36} routine. We found that a relative tolerance of ${10}^{-4}$ and an absolute tolerance of ${10}^{-10}$ in `double` precision gave good results. Though these choices noticeably extend calculation time, less demanding tolerances can produce numerically unsatisfying results. For example, modes with similar spatial frequencies should have similar total variance (integral of the PSD), but setting numerical tolerance too low can introduce occasional 50% jumps.

We also find that temporal frequency sampling is important, settling on $\mathrm{\Delta}f=0.1\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{Hz}$ as the largest sampling that should be used in general. Larger samplings tend to also produce jumps, caused by moving on and off sharp peaks associated with individual wind-layers. If such a peak falls between two frequencies and the sampling is too coarse, the power in the peak will be missed.

We illustrate calculations of multilayer turbulence using a seven layer turbulence model based on the GMT site survey at Las Campanas Observatory (LCO).^{37}^{,}^{38} The parameters of this profile are given in Table 2. The calculated temporal PSDs for several Fourier modes are shown in Fig. 6. Due to the identical wind velocities of several layers, we expect only four distinct peaks. The number of peaks seen at a given spatial frequency depends on ${f}_{pk,i}={\overrightarrow{\mathcal{V}}}_{i}\xb7{\overrightarrow{k}}_{m,n}$. If the frequency of two peaks is not well separated, they appear blended at the resolution of the calculations and plots.

## Table 2

LCO turbulence profile for median conditions assumed in this work. Based on the GMT site survey.37,38

Layer no. | Layer heightaz (m) | Layer strengthbCn2 | Wind DircΘ (deg) | Wind speed V (m/s) |
---|---|---|---|---|

1 | 250 | 0.42 | 60 | 10 |

2 | 500 | 0.03 | 60 | 10 |

3 | 1000 | 0.06 | 75 | 20 |

4 | 2000 | 0.16 | 75 | 20 |

5 | 4000 | 0.11 | 100 | 25 |

6 | 8000 | 0.10 | 110 | 30 |

7 | 16,000 | 0.12 | 100 | 25 |

## 5.3.

### Temporal PSD of Measurement Noise

We also need the PSD of the WFS measurement noise. The variance at a spatial frequency due to photon noise in the WFS is^{15}

## Eq. (48)

$${\sigma}_{\mathrm{ph},mn}^{2}=\frac{{\beta}_{p}^{2}({\overrightarrow{k}}_{mn})}{{N}_{\mathrm{ph}}}\text{\hspace{0.17em}\hspace{0.17em}}{(\mathrm{rad}}^{2}\text{\hspace{0.17em}}\mathrm{rms}),$$## Eq. (51)

$$S/N=\frac{{F}_{\gamma}{\tau}_{\mathrm{wfs}}}{\sqrt{{F}_{\gamma}{\tau}_{\mathrm{wfs}}+{n}_{\mathrm{px}}{F}_{\mathrm{bg}}{\tau}_{\mathrm{wfs}}+{n}_{\mathrm{px}}{\sigma}_{\mathrm{ron}}^{2}}},$$## Eq. (52)

$${\sigma}_{\mathrm{ph},mn}^{2}=\frac{{\beta}_{p}^{2}({\overrightarrow{k}}_{mn})}{{(S/N)}^{2}}\text{\hspace{0.17em}\hspace{0.17em}}({\mathrm{rad}}^{2}\text{\hspace{0.17em}}\mathrm{rms}).$$This variance is spread over the sampling bandwidth of the WFS, $1/2{\tau}_{\mathrm{wfs}}$. We can determine the temporal PSD of measurement noise from

## Eq. (53)

$${\int}_{0}^{1/2{\tau}_{\mathrm{wfs}}}{\mathcal{T}}_{\mathrm{ph},mn}(f)\mathrm{d}f={\sigma}_{\mathrm{ph},mn}^{2},$$Note that in the absence of detector noise, ${\mathcal{T}}_{\mathrm{ph},mn}(f)$ is independent of ${\tau}_{\mathrm{wfs}}$.

## 6.

## Closed-Loop Control

Now that we have established the temporal PSDs of the Fourier modes and of WFS measurement noise, we can use them to predict the residual wavefront variance in a closed-loop AO system. Here, we follow the standard model used throughout the AO literature and so only briefly introduce the concepts and our notation. Note that though here we apply the following to the frozen-flow turbulence PSDs we just derived, the techniques are general and could be applied to PSDs in any basis set and to PSDs including additional error sources such as telescope vibrations and nonfrozen-flow.

## 6.1.

### Transfer Functions

First, we need the error and noise transfer functions (NTF) of the AO control system. Here, we closely follow the development of Poyneer et al.,^{39} and see also the treatment by Madec.^{40} Our goal is to calculate the error transfer function (ETF) of the control system, which is defined as

## Eq. (55)

$${|\mathrm{ETF}(f)|}^{2}=\frac{{\mathcal{T}}_{mn,\text{out}}(f)}{{\mathcal{T}}_{mn,\text{in}}(f)}.$$That is, the ETF describes the AO system as a temporal filter which acts on the input PSD. The NTF is similarly defined, except that it acts on the measurement noise PSD. The ETF and NTF are constructed using the frequency responses, or transfer functions, of the various components of the system, which we will denote as $H$.

The WFS measurement is an average of the turbulence over the integration time. Likewise, the deformable mirror (DM) is assumed to move in discrete steps, holding its position between updates. The frequency response of such components is called a sample and hold, which has the functional form

where $T=1/{f}_{s}$, ${f}_{s}$ being the loop sampling frequency, and $s=i2\pi f$. There will be a delay $\tau $ due to sensor readout, data transfer, calculations, etc. It is common to discuss the “total delay” including 1/2 the WFS sample and hold (for the midpoint of the integration) and 1/2 the DM sample and hold. That is the total delay is $T+\tau $.^{40}The transfer function of the pure delay is

The system forms an estimate of the coefficient of a mode, $\tilde{h}$. In closed-loop, the WFS measures the change in the coefficient of a given mode, $\mathrm{\Delta}h$. This is combined with the current estimate of the total coefficient to form a new estimate of $\tilde{h}$, which is then applied to the DM. At time-step ${t}_{i}$, the coefficient is specified in terms of the current and previous $\mathrm{\Delta}h$ and $\tilde{h}$, by a linear combination

## Eq. (58)

$$\tilde{h}({t}_{i})=\sum _{j=1}^{J}{a}_{j}\tilde{h}({t}_{i-j})+g\sum _{l=0}^{L}{b}_{l}\mathrm{\Delta}h({t}_{i-l}),$$^{41}

## Eq. (59)

$${H}_{\mathrm{con}}(z;g)=\frac{g\sum _{l=0}^{L}{b}_{l}{z}^{-l}}{1+\sum _{j=1}^{J}{a}_{l}{z}^{-j}}$$With the above component frequency responses, we can write the open-loop ETF as

## Eq. (60)

$${\mathrm{ETF}}_{\mathrm{ol}}(s;g)={H}_{\mathrm{wfs}}(s){H}_{\mathrm{con}}(s;g){H}_{\tau}(s){H}_{\mathrm{dm}}(s).$$From this it follows that the closed-loop ETF is

The closed-loop NTF is

## Eq. (62)

$${\mathrm{NTF}}_{\mathrm{cl}}(s;g)={H}_{\mathrm{dm}}(s){H}_{\tau}(s){H}_{\mathrm{con}}(s;g){\mathrm{ETF}}_{\mathrm{cl}}(s;g).$$The residual PSD in closed-loop control is then the product of the modulus-squared ETF and the modal PSD, and the modulus-squared NTF and WFS noise PSD

## Eq. (63)

$${\mathcal{T}}_{\mathrm{cl},mn}(f;g)={\mathcal{T}}_{mn}(f){|{\mathrm{ETF}}_{\mathrm{cl}}(s;g)|}^{2}+{\mathcal{T}}_{\mathrm{ph},mn}(f){|{\mathrm{NTF}}_{\mathrm{cl}}(s;g)|}^{2}.$$The residual variance in this mode is just the integral of the residual PSD

## Eq. (64)

$${|{\sigma}_{mn}(g)|}^{2}={\int}_{0}^{{f}_{s}}{\mathcal{T}}_{\mathrm{cl},mn}(f;g)\mathrm{d}f.$$Note that, since we started with Eq. (43), this estimate of the variance explicitly includes the aperture and so directly describes the postcoronagraph intensity.

To employ these calculations, we next need to determine the parameters of Eq. (58) and then find the optimum gain.

## 6.2.

### Control Laws

Our next step in analyzing a closed-loop AO system is to choose the control law used to feed WFS measurements back to the DM. Here we analyze the simple integrator (SI) and the linear predictor (LP).

## 6.2.1.

#### Simple integrator

At time step ${t}_{i}$, the WFS measures the residual wavefront, $\mathrm{\Delta}{h}_{i}$. We combine this with the last estimate, ${\tilde{h}}_{i-1}$, which is the result of the previous cycle. The simplest estimate at the current step is the SI

This corresponds to a choice of $J=1$; ${a}_{1}=1$ and $L=0$; ${b}_{0}=1$ in Eq. (58). It is also common to choose ${a}_{1}<1$ (slightly), a so-called leak term, making this the “leaky integrator” control law. Here we will confine our analysis to the SI control law.

## 6.2.2.

#### Linear predictor

In general, the estimate formed by the SI controller will be slightly incorrect when actually applied due to the finite loop delay. Here, we consider a conceptually straightforward extension of the integrator to include a prediction. We note that more sophisticated methods exist, which we discuss after presenting our results.

One way to improve the accuracy of the estimate formed by the controller is with a filter of the form

We can determine a set of ${c}_{j}$ which is optimal in the least-squares sense from the temporal autocorrelation of the modal amplitude, ${\mathcal{R}}_{mn}(t)$, where $t$ is the lag, by solving the Normal or Yule–Walker equations^{41}^{,}^{42}

## Eq. (68)

$$\mathbf{R}=\left[\begin{array}{cccc}{\mathcal{R}}_{mn}(0)& {\mathcal{R}}_{mn}(1)& \cdots & {\mathcal{R}}_{mn}(N-1)\\ {\mathcal{R}}_{mn}(1)& {\mathcal{R}}_{mn}(0)& \vdots & {\mathcal{R}}_{mn}(N-2)\\ \vdots & \cdots & \ddots & \vdots \\ {\mathcal{R}}_{mn}(N-1)& {\mathcal{R}}_{mn}(N-2)& \cdots & {\mathcal{R}}_{mn}(0)\end{array}\right]$$## Eq. (69)

$$\overrightarrow{r}=\left[\begin{array}{c}{\mathcal{R}}_{mn}(1)\\ {\mathcal{R}}_{mn}(2)\\ \vdots \\ {\mathcal{R}}_{mn}(N)\end{array}\right].$$We can find the autocorrelation of the time-series of a mode from its PSD using the Wiener–Khinchin theorem,^{43} which states that the autocorrelation is the Fourier transform of the PSD

In Appendix C, we give a recipe for calculating the autocorrelation from a numerical PSD.

The solution vector $\overrightarrow{c}=[{c}_{0},{c}_{1},\cdots ,{c}_{N-1}]$ contains the optimal LP coefficients, which minimize the least-squares error of the output of Eq. (66). The matrix $\mathbf{R}$ is a symmetric Toeplitz matrix, and so the Normal equations can be solved very efficiently using Levinson’s recursion^{44}. A nice result of using Levinson’s recursion to solve the Yule–Walker equations is that the resultant filter is guaranteed to be stable.^{42}

To apply the LP in closed-loop, we use $\overrightarrow{c}$ as the coefficients in Eq. (58). That is, $L=N-1$, ${b}_{l}={c}_{l}$ and $J=N$, ${a}_{j}={c}_{j-1}$. We can then to determine the optimal loop gain $g$.

## 6.3.

### Gain Optimization

All that remains is to describe how to choose the optimum value of the gain. The first step is to identify the maximum stable gain, ${g}_{\mathrm{max}}$, which is the value of $g$ at which the system becomes unstable (this is sometimes called the “ultimate gain.”) Following Dessenne et al.,^{45} we evaluate stability of the control law using a Nyquist plot of the open-loop ETF, which plots the real and imaginary components of Eq. (60). In this plane, as long as the ETF does not encircle the point $(0,-1)$, then the system is stable—a result known as Nyquist’s stability criterion. We then can find the value of ${g}_{\mathrm{max}}$ by finding the value of $g$ which causes the curve to intersect $(0,-1)$. This is illustrated in Fig. 7.

From here we follow the standard recipe^{25}^{,}^{46} for using the input (or “open-loop”) PSD to optimize the closed-loop gain. The per-mode variance-minimizing optimum gain is given by solving

## Eq. (71)

$${g}_{\mathrm{opt}}=\underset{0\le g<{g}_{\mathrm{max}}}{\mathrm{arg}\text{\hspace{0.17em}}\mathrm{min}}{|{\sigma}_{mn}(g)|}^{2}.$$This is a straightforward optimization problem, which we solve using Brent’s method as implemented in the boost c++ library.

Integrator control with optimal gains is implemented in the SPHERE SAXO system, which uses a Karhunen–Loéve basis with optimal modal gains.^{47} The GPI instrument likewise uses optimal gains but with Fourier modes.^{39}

## 7.

## Closed-Loop Contrast Calculations

Given the temporal PSD, the ETF and NTF, and the optimal gain, Eq. (64) describes the residual variance at a given spatial frequency in a closed-loop AO system. We showed how this is related directly to the PSF contrast behind a coronagraph in Sec. 4. So, calculating the closed-loop long-exposure postcoronagraph contrast consists of the following steps:

1. Calculate the temporal PSDs at each $(m,n)$ according to Secs. 5.1 and 5.2. Due to the Hermitian symmetry of the Fourier basis, this needs to be done in only half the plane.

2. Calculate the temporal PSD of measurement noise according to Sec. 5.3.

3. Determine the optimum controller coefficients per Sec. 6.2, and the associated optimum gains per Sec. 6.3.

4. Determine the ETF and NTF per Sec. 6.1

5. Populate the residual PSD as follows:

6. Calculate the intensity:

7. Calculate the long-exposure Strehl ratio as the sum of variance over all spatial frequencies to find ${\sigma}_{\mathrm{tot}}^{2}$ and use the extended Marechal approximation

8. Calculate the contrast according to Eq. (42).

The first four steps are illustrated in Fig. 8. The PSD of a Fourier mode is shown in Fig. 8(a). The gain-optimized ETFs and NTFs, for both SI and LP controllers, are shown in Fig. 8(b). Figure 8(a) also shows the residual PSDs after applying the ETFs and NTFs to the input turbulence PSD.

We illustrate the complete framework in Fig. 9. The frequency response of each component is applied to the input PSD, and to the noise PSD, according to the ETF and NTF. The output residual PSD is then input to the coronagraph.

We applied our method to the cases of a 6.5-m circular aperture, here motivated by the in-development MagAO-X system.^{48} We also consider a 25.4-m telescope, for which we used the collecting area of the GMT (seven 8.4 m segments), but ignored details of the aperture in PSD and contrast calculations. The parameters of our calculations for these telescopes are listed in Table 3.

## Table 3

Parameters.

Symbol | Parameter | Value | Units | Notes |
---|---|---|---|---|

${f}_{s}$ | Max sampling frequency | 2000 | Hz | |

$\tau $ | Loop delay | 1.5 | Frames | |

$d$ | Actuator spacing | 0.135 | m | 48 across 6.5 m |

${\beta}_{p}$ | WFS sensitivity | $\sqrt{2}$ | For unmodulated pyramid | |

$\lambda $ | WFS wavelength | 800 | nm | Based on the MagAO WFS^{49} |

WFS filter width | 357 | nm | Based on the MagAO WFS^{49} | |

WFS throughput | 0.1 | ratio | ||

${\tau}_{\mathrm{wfs}}$ | WFS exposure time | 0.0005 | s | 2 kHz |

${\sigma}_{\mathrm{ron}}$ | WFS readout noise | 0.1 | Electrons | |

0 mag flux density @ 800 nm | $5.0\times {10}^{7}$ | $\text{Photons}/\mathrm{sec}/{\mathrm{m}}^{2}/\mathrm{nm}$ | Vega based | |

Sky background @ 800 nm | 20 | ${\mathrm{mags}/\mathrm{arcsec}}^{2}$ | ||

${r}_{0}$ | Fried param. | 0.16 | m | Median 0.62 in. at $0.5\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$^{38} |

$\overline{\mathcal{V}}$ | Weighted mean wind speed | 18.7 | m/s | $5/3$ moment^{19} |

Magellan clay | ||||

$D$ | Telescope diameter | 6.5 | m | |

${F}_{\gamma},0$ | 0 mag photon flux | $5.9\times {10}^{10}$ | $\text{Photons}/\mathrm{s}$ | Unobstructed aperture |

${n}_{\mathrm{px}}$ | WFS pixels | 7280 | Pixels | Four quadrants with $D/d$ across |

${F}_{\mathrm{bg}}$ | Background flux | 0.22 | $\text{Photons}/\mathrm{pix}/\mathrm{s}$ | Assumes 2 in. FOV for WFS |

GMT | ||||

$D$ | Telescope diameter | 25.4 | m | Edge to edge |

${F}_{\gamma},0$ | 0 mag photon flux | $6.6\times {10}^{11}$ | $\text{Photons}/\mathrm{s}$ | Vega, using GMT aperture |

${n}_{\mathrm{px}}$ | WFS pixels | 111,208 | Pixels | Four quadrants with $D/d$ across |

${F}_{\mathrm{bg}}$ | Background flux | 0.01 | $\text{Photons}/\mathrm{pix}/\mathrm{s}$ | Assumes 2 in. FOV for WFS |

Figure 10 shows the Strehl ratio versus guide star magnitude predicted for the 6.5-m AO system at 800 nm. Only residual turbulence is accounted for, including the uncorrected power (a.k.a. fitting error).

In Fig. 11, we show the contrast maps for the 6.5-m system observing 5th, 8th, and 10th magnitude stars, for both controllers at the WFS wavelength. To estimate the PSF contrast at other wavelengths, all of these results can be scaled by $1/{\lambda}^{2}$. The figure shows, as expected, deeper contrasts on brighter stars. The SI controller (left-hand column) exhibits the well known “wind butterfly,” and we see that the gain optimization process turns off modes on fainter stars, changing the shape of the dark hole. The LP controller does not show the “wind butterfly,” rather predictive control essentially eliminates it. It also results in deeper contrasts. Figure 12 shows the PSF contrast profiles both along and across the wind direction.

The same calculations were made for the 25.4 GMT-like telescope, shown in Figs. 13 and 14. The results are qualitatively similar, but show a nearly $\times 10$ improvement due to the larger aperture. This can be understood as the larger diameter causing the spatial-frequency sampling interval, $1/{D}^{2}$, to impose less power per $\lambda /D$, resulting in less residual power per mode.

For both the 6.5- and the 25.4-m telescopes, the gain in contrast at small separations provided by the LP controller is remarkable. On the 6.5 m observing a 0 mag star, the LP contrast is over 200 times better at 1 $\lambda /D$ and is nearly 10 times better for an eighth mag star. On the 25.4 m, the LP contrast is over 1400 times better at 1 $\lambda /D$ on a 0-mag star and is over 30 times better for an eighth mag star.

## 8.

## Discussion

## 8.1.

### Predictive Control

Predictive control has been discussed in many other studies. Here, we used a relatively simple form of it, rather than more sophisticated options. The LP controller provides a conceptually simple extension of the integrator controller and provides an efficient and easy to understand method for calculating its coefficients since we know the open-loop PSD. However, it is not the most sophisticated method proposed to-date in the literature and our intent here is not to propose it as the best one.

Dessenne et al.^{45}^{,}^{50} presented a method using a linear filter, essentially identical to the technique we used. Their method for optimizing the parameters was based on an explicit least squares minimization using past data, rather than the PSD, and allowed the coefficients $a$ and $b$ to differ. This technique was demonstrated on sky.

We recently presented a related predictive control strategy, based on empirical orthogonal functions (EOFs).^{51} EOFs are essentially a time-and-space Karhunen–Loève basis, and this method allows for including non-WFS measurements, such as accelerometers. Filter optimization is done on past time-domain data. We found essentially the same level of improvement in simulations.

Johnson et al.^{52} developed a control strategy based on predicting the wind-driven translation of turbulent layers. With simulations incorporating on-sky data from the ViLLaGEs AO system,^{53} they showed significant increase in Strehl on fainter stars and predicted large gains in Strehl for larger telescopes.

Much recent work has been in the area of linear quadratic Gaussian (LQG) controllers. Poyneer et al.^{24} proposed a Fourier-mode by Fourier-mode optimization of a Kalman filter controller, and using end-to-end simulations showed large gains in contrast similar to what we have shown. They presented their method as a control filter much as we did here. They pointed out that any filter could be cast in the direct form as in Eq. (59), though it is then subject to errors due to numerical precision if actually used to implement the filter.

More broadly, the determination of optimum control filters is of great interest in the AO-related literature. In the context of AO, this has been studied extensively for LQG controllers.^{54}^{,}^{55} See also the review of various system identification strategies presented by Kulcsár et al.^{56}

While this paper was under review, Correia et al. presented an analysis of predictive control using the distributed Kalman filter method.^{14} Their analysis employed the spatial filtering method (as opposed to our temporal frequency method) and showed postcoronagraph contrast gains for the Keck AO system at least as large as we show here.

Any of these control methods can be represented by their temporal transfer functions, and so could be included in further studies using our technique. In summary, the LP is not the only, nor the most advanced, technique we have available, but it is apparent that predictive control offers significant gains for high-contrast imaging of exoplanets.

## 8.2.

### Sensitivity to Frozen Flow

Our model of input turbulence was based on the frozen flow model. However, our semianalytic method does not require frozen flow. That is, any source of the temporal PSDs could be used as input to the analysis. We could incorporate boiling or disturbances derived from telescope models (e.g., vibrations due to wind shake). This could also be extended to assess the impact of time-variable turbulence profiles.

It should be expected that inclusion of additional sources of error will reduce projected performance. In the case of, e.g., telescope vibrations this is simply because there is more variance in the system. In the case of time-variable turbulence, i.e., changing wind and seeing, the reduced performance will result from suboptimality of the filter coefficients.

## 8.3.

### Impact on Exoplanet Characterization

Our motivation for this study was to analyze the limits of ground-based coronagraphic exoplanet characterization, and to assess the impact that predictive control will have on it. Consider the coupling of an ExAO-fed coronagraph to a high-dispersion spectrograph,^{57} a technique recently dubbed “high-dispersion coronagraphy” (HDC).^{58}^{,}^{59} We focus on this technique because it is, in principle, photon-noise limited, meaning we do not have to consider speckle noise. Neglecting detector and background noise, the signal-to-noise ratio ($S/N$) in the HDC technique is given as^{57}

## Eq. (74)

$$S/N=\frac{{F}_{\mathrm{p}}\mathrm{\Delta}t}{\sqrt{{F}_{*}\mathrm{\Delta}t\u27e8C\u27e9}}\sqrt{{N}_{\text{lines}}}\text{\hspace{0.17em}},$$So, in terms of required integration times, the potential benefits of predictive control are quite profound. Using just the relatively simple form of it considered in this study, we find that predictive control of atmospheric turbulence could make a GSMT up to 30 times more efficient at exoplanet characterization around a star as bright as Proxima Centauri, and up to 1400 times more efficient on brighter stars.

## 8.4.

### Code Availability

The code we developed for this analysis is freely available in a Github repository: https://github.com/jaredmales/aoSystem. It is written in c++, is highly parallelized, and compiles to a command line application which accepts configuration files to set up the analysis. It depends on the library available in a Github repository: https://github.com/jaredmales/mxlib, as well as several easily available open-source libraries. No licenses are needed.

## 9.

## Conclusion

We have developed a semianalytic framework for predicting the raw PSF contrast in an AO-fed coronagraphic image under closed-loop control. Instead of the usual spatial-filter analysis, we explicitly calculated temporal PSDs and analyzed AO performance in the temporal-frequency domain. Our technique takes into account closed-loop dynamics as well as a multilayer frozen-flow turbulence model. We also analyzed performance of ideal coronagraphs, showing how to account for the subtle difference between the residual PSD and the residual modal variance on an aperture. Using our model, we showed that predictive control will improve raw contrast by up to a factor of 1400 on bright stars, and by a factor of 30 on eighth magnitude stars, compared to the SI control law. Assuming a photon-noise limited observation, the exposure time required to characterize an exoplanet will decrease by the same large factors. This has significant implications for the potential of the coming ground-based GSMTs to search for life outside our solar system and makes predictive control one of the key techniques that must be perfected for exoplanet characterization.

## Appendices

## Appendix A:

### Mathematics of the Fourier Basis

We derive the normalizations and Fourier transforms of the Fourier basis on a circular aperture.

## A.1.

#### Normalization

The ${L}^{2}$ norm of the cosine mode on the aperture is

## Eq. (76)

$${|\mathrm{cos}(2\pi {\overrightarrow{k}}_{mn}\xb7\overrightarrow{q})|}^{2}=\int \mathcal{A}(\overrightarrow{q}){\mathrm{cos}}^{2}(2\pi {\overrightarrow{k}}_{mn}\xb7\overrightarrow{q})\mathrm{d}\overrightarrow{q}.$$For a symmetric aperture, we can restrict ourselves to $n=0$ with no loss of generality, which simplifies the integral to

## Eq. (77)

$${|\mathrm{cos}(2\pi {\overrightarrow{k}}_{mn}\xb7\overrightarrow{q})|}^{2}=\frac{4}{\pi {D}^{2}}{\int}_{0}^{D/2}{\int}_{0}^{2\pi}{\mathrm{cos}}^{2}\left(2\pi \frac{m}{D}q\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\theta \right)\mathrm{d}\theta q\text{\hspace{0.17em}}\mathrm{d}q.$$Using the exponential form of the cosine

## Eq. (78)

$${|\mathrm{cos}(2\pi {\overrightarrow{k}}_{mn}\xb7\overrightarrow{q})|}^{2}={\int}_{0}^{D/2}{\int}_{0}^{2\pi}\frac{1}{2}+\frac{1}{4}({e}^{4\pi \frac{m}{D}q\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\theta}+{e}^{-4\pi \frac{m}{D}q\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\theta})\mathrm{d}\theta q\text{\hspace{0.17em}}\mathrm{d}q.$$With the following identity:

## Eq. (79)

$$2\pi {J}_{0}(z)={\int}_{0}^{2\pi}{\mathrm{e}}^{\pm iz\text{\hspace{0.17em}}\mathrm{cos}(\varphi )}\mathrm{d}\varphi .$$(cf. Ref. 60 8.411 #7), the integral evaluates as

## Eq. (80)

$${|\mathrm{cos}(2\pi {\overrightarrow{k}}_{mn}\xb7\overrightarrow{q})|}^{2}=\frac{1}{2}+\frac{4}{{D}^{2}}{\int}_{0}^{D/2}{J}_{0}\left(4\pi \frac{m}{D}q\right)q\mathrm{d}q.$$Next using the identity

(cf. Ref. 60 5.52 #1) and a change of variables we have

## Eq. (82)

$${|\mathrm{cos}(2\pi {\overrightarrow{k}}_{mn}\xb7\overrightarrow{q})|}^{2}=\frac{1}{2}+\mathrm{Ji}(2\pi {k}_{mn}D).$$Similarly, the ${L}^{2}$ norm of the sine mode on the aperture is

and the norm of the product is## A.2.

#### Fourier Transforms on an Aperture

We need the Fourier transform of the Fourier modal basis over the aperture of the telescope, that is

## Eq. (85)

$${Q}_{mn}^{c}(\overrightarrow{k})={\int}_{-\infty}^{\infty}{\int}_{-\infty}^{\infty}\mathcal{A}(u,v)\mathrm{cos}(2\pi {\overrightarrow{k}}_{mn}\xb7\overrightarrow{q}){\mathrm{e}}^{i2\pi ({k}_{u}u+{k}_{v}v)}\mathrm{d}u\text{\hspace{0.17em}}\mathrm{d}v.$$The unobstructed circular aperture is described by Eq. (2).

Now switching to polar coordinates and collecting terms we have

## Eq. (86)

$${Q}_{mn}^{c}(\overrightarrow{k})=\frac{4}{2\pi {D}^{2}}{\int}_{0}^{D/2}{\int}_{0}^{2\pi}\{{\mathrm{e}}^{i2\pi q\left[\right({k}_{u}+\frac{m}{D})\mathrm{cos}\text{\hspace{0.17em}\hspace{0.17em}}\theta +({k}_{v}+\frac{n}{D}\left)\mathrm{sin}\text{\hspace{0.17em}}\theta \right]}+{\mathrm{e}}^{-i2\pi q\left[\right({k}_{u}-\frac{m}{D})\mathrm{cos}\text{\hspace{0.17em}}\theta +({k}_{v}-\frac{n}{D}\left)\mathrm{sin}\text{\hspace{0.17em}}\theta \right]}\}\mathrm{d}\theta q\text{\hspace{0.17em}}\mathrm{d}q,$$## Eq. (87)

$${k}_{mn}^{+}=\sqrt{{({k}_{u}+\frac{m}{D})}^{2}+{({k}_{v}+\frac{n}{D})}^{2}}\phantom{\rule{0ex}{0ex}}{k}_{mn}^{-}=\sqrt{{({k}_{u}-\frac{m}{D})}^{2}+{({k}_{v}-\frac{n}{D})}^{2}}$$## Eq. (88)

$$\mathrm{cos}\text{\hspace{0.17em}}\alpha =({k}_{u}+\frac{m}{D})/{k}_{mn}^{+}\phantom{\rule[-0.0ex]{1em}{0.0ex}}\mathrm{sin}\text{\hspace{0.17em}}\alpha =({k}_{v}+\frac{n}{D})/{k}_{mn}^{+}\phantom{\rule{0ex}{0ex}}\mathrm{cos}\text{\hspace{0.17em}}\beta =({k}_{u}-\frac{m}{D})/{k}_{mn}^{-}\phantom{\rule[-0.0ex]{1em}{0.0ex}}\mathrm{sin}\text{\hspace{0.17em}}\beta =({k}_{v}-\frac{n}{D})/{k}_{mn}^{-}.$$The transform can now be written

## Eq. (89)

$${Q}_{mn}^{c}(\overrightarrow{k})=\frac{2}{\pi {D}^{2}}{\int}_{0}^{D/2}{\int}_{0}^{2\pi}[{\mathrm{e}}^{i2\pi q{k}_{mn}^{+}\text{\hspace{0.17em}}\mathrm{cos}(\theta -\alpha )}+{\mathrm{e}}^{-i2\pi q{k}_{mn}^{-}\text{\hspace{0.17em}}\mathrm{cos}(\theta -\beta )}]\mathrm{d}\theta q\text{\hspace{0.17em}}\mathrm{d}q.$$With the following identity:

## Eq. (90)

$$2\pi {J}_{0}(z)={\int}_{0}^{2\pi}{\mathrm{e}}^{\pm iz\mathrm{cos}(\varphi )}\mathrm{d}\varphi ,$$## Eq. (91)

$${Q}_{mn}^{c}(\overrightarrow{k})=\frac{4}{{D}^{2}}{\int}_{0}^{D/2}[{J}_{0}(2\pi {k}_{mn}^{+}q)+{J}_{0}(2\pi {k}_{mn}^{-}q)]q\mathrm{d}q,$$## Eq. (93)

$${Q}_{mn}^{c}(\overrightarrow{k})=[\frac{{J}_{1}(\pi D{k}_{mn}^{+})}{\pi D{k}_{mn}^{+}}+\frac{{J}_{1}(\pi D{k}_{mn}^{-})}{\pi D{k}_{mn}^{-}}].$$Similarly for the sine modes

## Appendix B:

### Contrast *C*_{2} of Guyon

Here we present a derivation of the contrast due to time delay and measurement error from Ref. 15, there called ${C}_{2}$. We extend the analysis of Ref. 15 to include finite loop delays, detector, and background noise, and make it consistent with the results of Sec. 4. A spatial frequency component of the wavefront phase at time $t$ can be written as

## Eq. (95)

$${\mathrm{\varphi}}_{mn}(\overrightarrow{q},t)=\frac{2\pi}{\lambda}{h}_{mn}^{\u2020}(t)\mathrm{cos}[{\overrightarrow{k}}_{mn}\xb7\overrightarrow{q}+{\varphi}_{mn}^{\u2020}(t)],$$## Eq. (96)

$${\mathrm{\varphi}}_{mn}(\overrightarrow{q},t+{\tau}_{\mathrm{tl}})\approx {\mathrm{\varphi}}_{mn}(\overrightarrow{q},t)+\frac{\partial {\mathrm{\varphi}}_{mn}(\overrightarrow{q},t)}{\partial t}{\tau}_{\mathrm{tl}}.$$The residual variance after applying the wavefront correction will be

## Eq. (97)

$${\sigma}_{\mathrm{tl},mn}^{2}\approx \u27e8{\left[\frac{\partial {\mathrm{\varphi}}_{mn}(\overrightarrow{q},t)}{\partial t}\right]}^{2}\u27e9{\tau}_{\mathrm{tl}}^{2}\text{\hspace{0.17em}},$$## Eq. (98)

$${\sigma}_{\mathrm{tl},mn}^{2}\approx {\left(\frac{2\pi}{\lambda}\right)}^{2}\{\u27e8{\left[\frac{d{h}_{mn}^{\u2020}(t)}{dt}\right]}^{2}\u27e9+\u27e8{[{h}_{mn}^{\u2020}(t)\frac{d{\varphi}_{mn}^{\u2020}(t)}{dt}]}^{2}\u27e9\}{\tau}_{\mathrm{tl}}^{2}.$$We next employ the frozen flow hypothesis, assuming that the amplitude does not change and that the phase changes only due to wind-driven flow. This gives

## Eq. (99)

$$\frac{d{h}_{mn}^{\u2020}(t)}{dt}=0\phantom{\rule{0ex}{0ex}}\frac{d{\varphi}_{mn}^{\u2020}(t)}{dt}=2\pi {\overrightarrow{k}}_{mn}\xb7\overrightarrow{\mathcal{V}},$$## Eq. (100)

$${\sigma}_{\mathrm{tl},mn}^{2\u2020}\approx {\left(\frac{2\pi}{\lambda}\right)}^{2}\u27e8{h}_{mn}^{\u20202}(t)\u27e9{(2\pi {\overrightarrow{k}}_{mn}\xb7\overrightarrow{\mathcal{V}})}^{2}{\tau}_{\mathrm{tl}}^{2}.$$Now the variance of each mode will be given by the spatial PSD according to

## Eq. (101)

$$\u27e8{h}_{mn}^{\u23e72}(t)\u27e9=\frac{\mathcal{P}({\overrightarrow{k}}_{mn})}{{D}^{2}}{\left(\frac{{\lambda}_{0}}{2\pi}\right)}^{2}\text{\hspace{0.17em}\hspace{0.17em}}({\mathrm{m}}^{2}\text{\hspace{0.17em}}\mathrm{rms}).$$Employing the wavefront variance from measurement noise given by Eq. (52), the total residual variance at the science wavelength will be

## Eq. (102)

$${\sigma}_{mn}^{\u20202}={\left(\frac{{\lambda}_{0}}{\lambda}\right)}^{2}\frac{\mathcal{P}(\overrightarrow{k})}{{D}^{2}}{(2\pi {\overrightarrow{k}}_{mn}\xb7\overrightarrow{\mathcal{V}})}^{2}{\tau}_{\mathrm{tl}}^{2}+\frac{{\beta}_{p}^{2}({\overrightarrow{k}}_{mn})}{S/{N}^{2}}{\left(\frac{{\lambda}_{\mathrm{wfs}}}{\lambda}\right)}^{2}\text{\hspace{0.17em}\hspace{0.17em}}({\mathrm{rad}}^{2}\text{\hspace{0.17em}}\mathrm{rms}).$$In this open-loop framework, the equivalent time lag is ${\tau}_{\mathrm{tl}}={\tau}_{\mathrm{wfs}}+\tau $ where $\tau $ is the closed-loop delay defined in Sec. 6.1, and ${\tau}_{\mathrm{wfs}}$ accounts for the sample-and-hold of the WFS and DM.^{40} We calculate the optimum frame rate that minimizes the variance by differentiating Eq. (102) with respect to ${\tau}_{\mathrm{wfs}}$, which yields a quartic equation

## Eq. (103)

$${\lambda}_{0}^{2}\frac{\mathcal{P}({\overrightarrow{k}}_{mn})}{{D}^{2}}{(2\pi {\overrightarrow{k}}_{mn}\xb7\overrightarrow{\mathcal{V}})}^{2}[{\tau}_{\mathrm{wfs}}^{4}+\tau {\tau}_{\mathrm{wfs}}^{3}]-\frac{{\lambda}_{\mathrm{wfs}}^{2}{\beta}_{p}^{2}({\overrightarrow{k}}_{mn})}{{F}_{\gamma}^{2}}[({F}_{\gamma}+{n}_{\mathrm{px}}{F}_{\mathrm{bg}}){\tau}_{\mathrm{wfs}}+2{n}_{\mathrm{px}}{\sigma}_{\mathrm{ron}}^{2}]=0.$$Each mode will have its own optimum ${\tau}_{\mathrm{wfs}}$, which is somewhat analogous to modal gains in closed-loop control. Since $\u27e8{h}_{mn}^{p2}\u27e9=\u27e8{h}_{mn}^{c\u20202}\u27e9+\u27e8{h}_{mn}^{s\u20202}\u27e9=\u27e8{h}_{mn}^{\u20202}\u27e9$, we can use the results of Sec. 4 to write

## Eq. (104)

$$\u27e8{I}_{\mathrm{\varphi},mn}^{\u2020}(\overrightarrow{r})\u27e9={\sigma}_{mn}^{\u20202}[\mathrm{PSF}(\overrightarrow{r}-{\overrightarrow{k}}_{mn}\lambda )+\mathrm{PSF}(\overrightarrow{r}+{\overrightarrow{k}}_{mn}\lambda )].$$The summation in Eq. (39) defines $\u27e8{I}_{\mathrm{\varphi}}^{\u2020}(r)\u27e9$, which leads to the expression for the contrast ${C}_{2}$

## Eq. (105)

$${C}_{2}(\overrightarrow{r})=\frac{\u27e8{I}_{\mathrm{\varphi}}^{\u2020}(\overrightarrow{r})\u27e9}{\u27e8S\u27e9\mathrm{PSF}(0)}.$$The Strehl ratio can be calculated by summing the total WFE after the convolution and then employing the Marechal approximation

It is straightforward to derive the other contrasts, ${C}_{\mathrm{0,1},\mathrm{3,4},\mathrm{5,6}}$, of Ref. 15 using similar arguments.

## Appendix C:

### Calculating the Autocorrelation

The autocorrelation $\mathcal{R}$ can be calculated from a numerical PSD using the following recipe:

1. Start with temporal PSD ${\mathcal{T}}_{mn}(f)$ at $N$ discrete frequencies ${f}_{j}$, where ${f}_{1}=\mathrm{\Delta}f$ and ${f}_{N}=\frac{1}{2}{f}_{s}$.

2. Form the two-sided temporal PSD of length $2N$ by augmenting with 0 at the beginning, and augmenting with the reverse at the end

3. Calculate the DFT of the augmented PSD using the FFT algorithm, giving the autocorrelation at lag $\tau $

Note that this is will be a circular autocorrelation wrapping around after $N$ points. $N$ should be chosen to be larger (by at least a factor of 2) than the needed number of points of $\mathcal{R}$.

4. The lags of the points are given by

where only the first $N+1$ points are unique.## Eq. (109)

$${t}_{{j}^{\prime}}=\{\begin{array}{cc}\frac{{j}^{\prime}-1}{2N{f}_{s}},& {j}^{\prime}=1\dots N+1\\ \frac{{j}^{\prime}-2N-1}{2N{f}_{s}},& {j}^{\prime}=N+2\dots 2N\end{array},$$

## Acknowledgments

This work was motivated by a question from Marcos van Dam. We thank Richard Frazin for reviewing an early draft of this manuscript. We thank the two anonymous referees for their careful and insightful reviews, which significantly improved this manuscript. Early steps in this work were performed under contract with the California Institute of Technology (Caltech)/Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. The authors acknowledge support from NSF awards 1506818 and 1625441.

## References

## Biography

**Jared R. Males** is an astronomer at the University of Arizona’s Steward Observatory. He received his PhD from the University of Arizona, and he was a NASA Sagan Fellow there from 2013 to 2016. His research is focused on exoplanet characterization and high contrast imaging instrumentation. He is the principal investigator of the NSF-funded MagAO-X extreme adaptive optics system.

**Olivier Guyon** is an astronomer at the Subaru Telescope, the University of Arizona, and the Japanese Astrobiology Center. He specializes in development of optical techniques to detect and study exoplanets. He develops high contrast imaging techniques combining coronagraphs and extreme adaptive optics to image nearby exoplanets with space and ground telescopes. His research interests also include advanced adaptive optics control techniques and calibration for high contrast imaging.