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

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

3.2.

Spatial Power Spectral Densities

The spatial PSD of the phase at the observation wavelength $\lambda$ in the von Kármán model¹⁹ is

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

So the piston subtracted phase PSD is

We model atmospheric turbulence as occurring in discrete layers. To account for the Fresnel propagation among these layers, we use the functions¹⁵

3.3.

Representing the Wavefront

3.4.

Statistics of the Fourier Basis

Following the recipe of Noll,²³ we can write the covariance of any two modes over the aperture as

We now make use of this expression to analyze the statistics of the Fourier basis under von Kármán turbulence.

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⟨{|h_{mn}^p|}^2⟩={\int }_0^{\infty }{\int }_0^{2\pi }{|Q_{mn}^p|}^2\mathcal{P}(k)\mathrm{d}\theta k\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}^{-})],$$

independent of $p$, it follows that $⟨{|h_{mn}^{+}|}^2⟩=⟨{|h_{mn}^{-}|}^2⟩$. From now on, we will drop the $p$ superscript except where it is explicitly needed.

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²⁶) 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.

Fig. 1

Modal variance versus spatial frequency. We compare the result of Eq. (25) (solid lines) with the simple estimate of ${\mathcal{P}}_{\mathrm{\varphi }}/D^2$ (dashed lines). The black curves are for uncorrected turbulence, and the red curves are for a simple treatment of AO correction. The discrepancy between the solid and dashed curves is due to the convolution on the aperture, which is neglected in the simple estimate.

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¹⁵ (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},$$

where $\mathcal{V}$ is the wind speed, ${\tau }_{\mathrm{tl}}$ is the total latency in the control system, and $d$ is the actuator spacing projected on the aperture. We next recalculated the modal variances with the corrected PSD, using $V=10\mathrm{m}/\mathrm{s}$, ${\tau }_{\mathrm{tl}}=1.25\mathrm{ms}$, and $d=13.5\mathrm{cm}$ on a $D=6.5\mathrm{m}$ aperture.

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.

Fig. 2

Correlations between the Fourier modes in uncorrected von Kármán turbulence. (a) Modes with $p^{\prime }=p$ and (b) modes with $p^{\prime }\ne p$. Each covariance point is normalized to the associated diagonal in (a), which is the variance of each mode. In uncorrected turbulence, the correlations are high relative to the modal variances.

Fig. 3

Same as Fig. 2, but for AO corrected von Kármán turbulence. After AO correction removes low-spatial-frequency power, the modal correlations are much smaller but not zero.

4.1.

Coronagraph Model

The intensity in an image plane is given as

The operator $\mathcal{F}\{·\}$ 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.²⁷ 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

In the case of the circular unobstructed aperture [Eq. (2)], this is the airy pattern

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³⁰

Macintosh et al.³⁰ only addressed phase, but the extension to amplitude is straightforward from their arguments. Dealing with the phase component first, we have

Evaluating a single term gives

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$)

Fig. 4

Illustration of the relationship between a phase aberration defined by a Fourier mode with spatial frequency ${\overrightarrow{k}}_{mn}$ and coefficient $h_{mn}$, to the long-exposure contrast of the resultant speckles in the postcoronagraph image plane. This relationship does not fully define the contrast in the presence of many speckles.

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

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

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, $⟨I_A(\overrightarrow{r})⟩$, 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, $⟨I_{\mathrm{NC}}(0)⟩$. This is related to the PSF by the long-exposure Strehl ratio $⟨S⟩$, 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.”³² 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,¹⁵ 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³⁵). 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\mathrm{m}$, an outer scale $L_0=25\mathrm{m}$, and windspeed of $V=10\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\mathrm{m}$ and actuator spacing to $d=13.5\mathrm{cm}$ (48 across) and assumed a wavelength of $\lambda =800\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.

Fig. 5

Numerical verification of coronagraph models. Here, we compare PSF contrasts for the perfect coronagraph,³² the coronagraph defined by Eq. (34), the APLCMC,³⁵ the contrasts $C_1+C_2$ from Guyon,¹⁵ and the variance of the Fourier modes as in Eq. (40). The left graph is for residual atmospheric phase aberrations only, and the right includes amplitude errors. The bottom graph shows the fractional errors of the model coronagraph w.r.t. the APLCMC. In general, the analytic coronagraph models match the physically realizable APLCMC well. This shows that such models are valid for analyzing the potential of ground-based AO-fed coronagraphs.

4.3.

Long-Exposure Contrast and the PSD

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

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.²⁷ 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.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$)

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

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

The temporal PSD for arbitrary wind direction is then simply

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³⁶ 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\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·{\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

Fig. 6

Calculated temporal PSDs in a seven-layer turbulence model at four different modes, corresponding to discrete spatial frequencies.

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¹⁵

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

Note that in the absence of detector noise, ${\mathcal{T}}_{\mathrm{ph},mn}(f)$ is independent of ${\tau }_{\mathrm{wfs}}$.

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.,³⁹ and see also the treatment by Madec.⁴⁰ Our goal is to calculate the error transfer function (ETF) of the control system, which is defined as

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

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

With the above component frequency responses, we can write the open-loop ETF as

From this it follows that the closed-loop ETF is

The closed-loop NTF is

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

The residual variance in this mode is just the integral of the residual PSD

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.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_{\max }$, which is the value of $g$ at which the system becomes unstable (this is sometimes called the “ultimate gain.”) Following Dessenne et al.,⁴⁵ 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_{\max }$ by finding the value of $g$ which causes the curve to intersect $(0,-1)$. This is illustrated in Fig. 7.

Fig. 7

Nyquist plot demonstrating the stability of our two control laws. The open-loop ETFs are plotted in the imaginary plane for both the SI and LP control laws. Note that frequency increases as the lines move out from the center. The controllers are stable so long as they do not encircle $(0,-1)$. The value of gain $g$ which causes them to intersect $(0,-1)$, and hence become unstable, is here called the maximum stable gain $g_{\max }$. The black circle has a radius of 1. We show the open-loop ETFs for each controller, with optimum gains for both a zeroth and an eighth magnitude star.

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

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.⁴⁷ The GPI instrument likewise uses optimal gains but with Fourier modes.³⁹

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).⁵¹ 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.⁵² 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,⁵³ 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.²⁴ 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.⁵⁶

While this paper was under review, Correia et al. presented an analysis of predictive control using the distributed Kalman filter method.¹⁴ 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,⁵⁷ 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⁵⁷

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.