Translator Disclaimer
6 February 2018 Ground-based adaptive optics coronagraphic performance under closed-loop predictive control
Author Affiliations +
The discovery of the exoplanet Proxima b highlights the potential for the coming generation of giant segmented mirror telescopes (GSMTs) to characterize terrestrial—potentially habitable—planets orbiting nearby stars with direct imaging. This will require continued development and implementation of optimized adaptive optics systems feeding coronagraphs on the GSMTs. Such development should proceed with an understanding of the fundamental limits imposed by atmospheric turbulence. Here, we seek to address this question with a semianalytic framework for calculating the postcoronagraph contrast in a closed-loop adaptive optics system. We do this starting with the temporal power spectra of the Fourier basis calculated assuming frozen flow turbulence, and then apply closed-loop transfer functions. We include the benefits of a simple predictive controller, which we show could provide over a factor of 1400 gain in raw point spread function contrast at 1 λ/D on bright stars, and more than a factor of 30 gain on an I=7.5  mag star such as Proxima. More sophisticated predictive control can be expected to improve this even further. Assuming a photon-noise limited observing technique such as high-dispersion coronagraphy, these gains in raw contrast will decrease integration times by the same large factors. Predictive control of atmospheric turbulence should therefore be seen as one of the key technologies that will enable ground-based telescopes to characterize terrestrial planets.



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 mission5 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-110 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.1114 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.



This analysis is concerned with AO in the “extreme” regime: high actuator counts (2000 on a 6.5-m aperture, 21,000 on a 25.4 m) and high loop speeds (2000  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 filter18 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.


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.

qmPosition vector in the pupil plane
u,vmPosition coordinates in the pupil plane
u^,v^Coordinate unit vectors in the pupil plane
km1Spatial frequency vector in the pupil plane
Δkm1Discrete spatial frequency sampling
ku,kvm1Spatial frequency components in the pupil plane
m,nInteger indices of discrete spatial frequencies
DmTelescope diameter
A(q)Aperture function
λmWavelength of observation and wavefront sensing
Prad2m2Spatial power spectrum of the wavefront
λ0mReference wavelength for turbulence parameters
r0mFried parameter for atmospheric turbulence strength
QDenotes the Fourier transform of a mode
zmHeight above the observatory
Cn2(zi)Normalized turbulence strength for layer i
MmnDenotes a Fourier mode
hmnmAmplitude of a mode describing the wavefront phase
amnAmplitude of a mode describing the wavefront amplitude
ϕ(q,t)radThe phase of the wavefront
A(q,t)The amplitude of the wavefront
σmn2rad2The variance of the amplitude of a mode
Vm/sWind speed
ΘradWind direction
τsDenotes a time delay
rλ/DPosition vector in the image plane
Iphotons/s/(λ/D)2Intensity in the image plane
PSFphotons/s/(λ/D)2Point-spread function
SStrehl ratio
CContrast ratio
fHzTemporal frequency
Tmnrad2  /HzTemporal power spectrum of a single mode
H(s)Denotes a transfer function


Coordinate System

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

Eq. (1)


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

Eq. (2)

A(q)={4πD2,if  qD20,if  q>D2.

Spatial-frequency is a vector quantity

Eq. (3)


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

Eq. (4)


Hence, we will often discretize spatial-frequency as

Eq. (5)

where m and n are integer indices.


Spatial Power Spectral Densities

The spatial PSD of the phase at the observation wavelength λ in the von Kármán model19 is

Eq. (6)

PvK(k,λ)=0.0218(λ0λ)21r05/31(k2+k02)11/6  (rad2/m2),
where r0 is Fried’s parameter20 and λ0 is the wavelength at which r0 is reported. The outer scale, L0, is included in the parameter

Eq. (7)


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 as21,22

Eq. (8)

where Po(k) is the input PSD ignoring the aperture, and Qmode denotes the Fourier transform of the mode on the aperture. As a Zernike polynomial, piston is simply Zpiston=1, and its Fourier transform on A(q) is23

Eq. (9)


So the piston subtracted phase PSD is

Eq. (10)

where the Jinc function is

Eq. (11)

and J1 is the cylindrical Bessel function of the first kind.

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

Eq. (12)

where Cn2(z) is the normalized turbulence strength profile as a function of altitude z. With these corrections, we have the PSD of the phase

Eq. (13)

and the PSD of the amplitude

Eq. (14)

at the aperture of the telescope.


Representing the Wavefront


Fourier basis

Here, we develop a real-valued form of the Fourier basis. Consider the complex exponential form of the basis24

Eq. (15)

where we have included a complex-valued coefficient h to illustrate the form of this function. Decomposing this into sines and cosines, we have

Eq. (16)

where Re{·} and Im{·} denote the real and imaginary parts, respectively. This suggests a real-valued Fourier basis defined as

Eq. (17)

where p=±1 (in the superscript we will use + and ). These modes are neither odd nor even and are normalized

Eq. (18)

although not orthogonal, on the aperture. For completeness, we show how the normalization is derived in Appendix A.1. The Fourier transforms on the aperture are (see Appendix A.2)

Eq. (19)

where for notational simplicity, we have defined

Eq. (20)


We have chosen this form of the Fourier basis primarily because it yields the useful result

Eq. (21)


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.


The Wavefront in the Fourier basis

We denote the wavefront phase as ϕ(q,t), which has units of radians. The coefficient of a mode is calculated as

Eq. (22)


These coefficients are real-valued. Poyneer and Veran25 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)


Similarly, we represent the wavefront amplitude aberrations as

Eq. (24)



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)


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


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

Eq. (26)

where δ(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

Eq. (27)

σmn2=Pϕ(kmn,λ)(Δk)2=Pϕ(kmn,λ)D2  (rad).


Modal variance

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

Eq. (28)


Now since

Eq. (29)

independent of p, it follows that |hmn+|2=|hmn|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 (GSL26) routines gsl_integration_qagiu (for the radial integral) and gsl_integration_qag (for the azimuthal integral). We found that an absolute tolerance of 1010 and a relative tolerance of 104 gave good results.

Fig. 1

Modal variance versus spatial frequency. We compare the result of Eq. (25) (solid lines) with the simple estimate of Pϕ/D2 (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 Guyon15 (as modified in Appendix B), the corrected PSD is approximated as

Eq. (30)

where V is the wind speed, τ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  m/s, τtl=1.25  ms, and d=13.5  cm on a D=6.5  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.



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=p and (b) modes with pp. 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.



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.


Coronagraph Model

The intensity in an image plane is given as

Eq. (31)


The operator 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.27 We define the system PSF as the result of the above with no aberrations, that is hmn=0 and amn=0 for all m,n, which is

Eq. (32)


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

Eq. (33)

where we have used the image plane coordinate r=λ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 expansion28,29 to approximate the postcoronagraph PSF as30

Eq. (34)


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)


Evaluating a single term gives

Eq. (36)


Then taking the modulus squared and pulling out the terms with mn=mn gives (we are suppressing the t dependence of h)

Eq. (37)

which is the instantaneous intensity at a position r and time t in the postcoronagraph image plane. Now, we seek the long exposure, or expected, value of this expression. Though this appears quite intractable, we see that this will require evaluating covariances of the form hmnphmnp. So, to proceed we make the ansatz that we can temporarily ignore the aperture as in Sec. 3.4.1. The effect of this assumption is that all covariances are 0, which leads directly to

Eq. (38)

for a single term in the sum on the first line of Eq. (37). This is the classic result of a Fourier aberration producing a pair of symmetric speckles.31 We illustrate this in Fig. 4.

Fig. 4

Illustration of the relationship between a phase aberration defined by a Fourier mode with spatial frequency kmn and coefficient hmn, 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 ϕ at r is

Eq. (39)


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)


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, IA(r), is evaluated in nearly identical fashion.

We define the “raw PSF contrast” as the ratio of the intensity at r to the peak intensity of a noncoronagraphic long-exposure image of the star, INC(0). This is related to the PSF by the long-exposure Strehl ratio S, defined by the relationship

Eq. (41)


The contrast is then

Eq. (42)



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 C2 of Guyon,15 combined with the contrast from uncorrected amplitude C1 due to atmospheric scintillation. In Appendix B, we present a derivation of C2 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 (APLCMC35). 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 λ/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 r0=0.2  m, an outer scale L0=25  m, and windspeed of V=10  m/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  m and actuator spacing to d=13.5  cm (48 across) and assumed a wavelength of λ=800  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 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,32 the coronagraph defined by Eq. (34), the APLCMC,35 the contrasts C1+C2 from Guyon,15 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.



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.


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.


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)


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 with21,22

Eq. (44)

where f is temporal frequency, V is the wind-speed, and Tmn denotes the temporal PSD of the mode specified by mn. With no loss of generality, we have chosen coordinates such that the wind vector V=Vu^. We have included a factor of 2 since we will consider only f>0.

Now to account for an arbitrary wind direction, let Θ be the direction of the wind with respect to the 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)


The temporal PSD for arbitrary wind direction is then simply

Eq. (46)


Finally, the temporal PSD of a Fourier mode for a turbulence profile is the Cn2 weighted sum of the single layer profiles

Eq. (47)



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_qagi36 routine. We found that a relative tolerance of 104 and an absolute tolerance of 1010 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 Δf=0.1  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 fpk,i=Vi·km,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 strengthbCn2Wind DircΘ (deg)Wind speed V (m/s)

aHeight above the observatory.


cAngle east of north.

Fig. 6

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



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 is15

Eq. (48)

σph,mn2=βp2(kmn)Nph  (rad2rms),
where the factor βp(kmn) quantifies the sensitivity of the WFS to photon noise. Nph is the total number of photons available for wavefront sensing

Eq. (49)

where Fγ is the photon rate (photons s1) sensed by the WFS, and τwfs is the integration time of the WFS. Fγ depends on star brightness through

Eq. (50)

where mag is the Vega magnitude of the star in the WFS bandpass, and Fγ,0 is the photon rate of a 0 mag star in the WFS bandpass including total system throughput. Here, we expand this to include background and readnoise. The signal-to-noise ratio is

Eq. (51)

where npx is the number of pixels involved in measuring Nph, each pixel having readout noise with variance σron2, and Fbg is the photon rate per pixel from the background. Now the WFE due to photon noise at the WFS wavelength can be written as

Eq. (52)

σph,mn2=βp2(kmn)(S/N)2  (rad2rms).

This variance is spread over the sampling bandwidth of the WFS, 1/2τwfs. We can determine the temporal PSD of measurement noise from

Eq. (53)

assuming white noise (a flat PSD) we have

Eq. (54)


Note that in the absence of detector noise, Tph,mn(f) is independent of τwfs.


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.


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)


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

Eq. (56)

where T=1/fs, fs being the loop sampling frequency, and s=i2πf. There will be a delay τ 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+τ.40 The transfer function of the pure delay is

Eq. (57)


The system forms an estimate of the coefficient of a mode, h˜. In closed-loop, the WFS measures the change in the coefficient of a given mode, Δh. This is combined with the current estimate of the total coefficient to form a new estimate of h˜, which is then applied to the DM. At time-step ti, the coefficient is specified in terms of the current and previous Δh and h˜, by a linear combination

Eq. (58)

where g is the loop gain. Below we consider in depth how to choose the parameters of this control filter. The frequency response of this general linear filter is41

Eq. (59)

and we can map from the z-domain to the s-domain by the substitution zesT.

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

Eq. (60)


From this it follows that the closed-loop ETF is

Eq. (61)


The closed-loop NTF is

Eq. (62)


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)


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

Eq. (64)


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.


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


Simple integrator

At time step ti, the WFS measures the residual wavefront, Δhi. We combine this with the last estimate, h˜i1, which is the result of the previous cycle. The simplest estimate at the current step is the SI

Eq. (65)


This corresponds to a choice of J=1; a1=1 and L=0; b0=1 in Eq. (58). It is also common to choose a1<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.


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

Eq. (66)


We can determine a set of cj which is optimal in the least-squares sense from the temporal autocorrelation of the modal amplitude, Rmn(t), where t is the lag, by solving the Normal or Yule–Walker equations41,42

Eq. (67)


Eq. (68)


Eq. (69)


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

Eq. (70)


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

The solution vector c=[c0,c1,,cN1] contains the optimal LP coefficients, which minimize the least-squares error of the output of Eq. (66). The matrix R is a symmetric Toeplitz matrix, and so the Normal equations can be solved very efficiently using Levinson’s recursion44. 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 c as the coefficients in Eq. (58). That is, L=N1, bl=cl and J=N, aj=cj1. We can then to determine the optimal loop gain g.


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, gmax, 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 gmax 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 gmax. 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 recipe25,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)


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


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:

    Eq. (72)

    σmn2={|σmn(gopt)|2,if  gopt>0.Pmn(k)/D2,if  gopt=0.

  • 6. Calculate the intensity:

    • a. For the corrected modes, the variance is properly sampled and directly describes the intensity according to Eq. (40).

    • b. For the uncorrected modes, we must use Eq. (38) and then Eq. (39). Note that we include the corrected-mode variances in the convolution.

  • 7. Calculate the long-exposure Strehl ratio as the sum of variance over all spatial frequencies to find σtot2 and use the extended Marechal approximation

    Eq. (73)


  • 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.

Fig. 8

(a) The input turbulent PSD for Fourier mode (m=12,n=12) (black), along with the residual PSDs in closed-loop control for the SI (blue) and LP (red) controllers. The WFS noise PSD (dashed) corresponds to an eighth mag star on a 6.5-m telescope. (b) The ETF (solid) and NTF (dashed) for both SI (blue) and LP (red) controllers, optimized for the input PSDs shown in (a).


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.

Fig. 9

Illustration of our framework, using one input PSD as an example. The input wavefront, described by its temporal PSD, is corrected by the DM. The WFS measures the residual, and noise is added as described by a temporal PSD. A controller, represented as a linear filter, determines an estimate of the wavefront at the next time-step. After a delay, this estimate determines the new shape of the DM. The control loop acts on the input PSD according to the ETF and on the noise PSD according to the NTF. The resultant residual temporal PSD determines the residual variance of this mode, which is then input to the coronagraph. Blue curves are for the SI and red curves are for the LP.


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


fsMax sampling frequency2000Hz
τLoop delay1.5Frames
dActuator spacing0.135m48 across 6.5 m
βpWFS sensitivity2For unmodulated pyramid
λWFS wavelength800nmBased on the MagAO WFS49
WFS filter width357nmBased on the MagAO WFS49
WFS throughput0.1ratio
τwfsWFS exposure time0.0005s2 kHz
σronWFS readout noise0.1Electrons
0 mag flux density @ 800 nm5.0×107Photons/sec/m2/nmVega based
Sky background @ 800 nm20mags/arcsec2
r0Fried param.0.16mMedian 0.62 in. at 0.5  μm38
V¯Weighted mean wind speed18.7m/s5/3 moment19
Magellan clay
DTelescope diameter6.5m
Fγ,00 mag photon flux5.9×1010Photons/sUnobstructed aperture
npxWFS pixels7280PixelsFour quadrants with D/d across
FbgBackground flux0.22Photons/pix/sAssumes 2 in. FOV for WFS
DTelescope diameter25.4mEdge to edge
Fγ,00 mag photon flux6.6×1011Photons/sVega, using GMT aperture
npxWFS pixels111,208PixelsFour quadrants with D/d across
FbgBackground flux0.01Photons/pix/sAssumes 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).

Fig. 10

Strehl ratio versus guide star magnitude with optimized gains for the two controllers at the WFS wavelength (800 nm). Only residual modal variance and fitting error are included in this plot. The LP controller improves Strehl at all magnitudes, with the benefit being greatest between 8th and 10th magnitude.


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/λ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.

Fig. 11

Postcoronagraph long-exposure contrast maps for a 6.5-m telescope at 800 nm. Here, we compare results for the SI and LP controllers, for 5th, 8th, and 10th magnitude stars. As we do not address chromaticity here, the WFS and observation wavelengths are both 800 nm.


Fig. 12

Contrast profiles for a 6.5-m telescope observing at 800 nm, using the parameters shown in Tables 2 and 3. As we do not address chromaticity here, the WFS and observation wavelengths are both 800 nm. (a) is along the main wind direction (up-and-down in Fig. 11), and (b) is perpendicular. The dashed lines are for the SI controller, standard in current AO systems. The solid lines are for the LP controller, the form of predictive control we analyze in this work. Since we are conducting a mode-by-mode analysis in discrete spatial frequencies, the lowest spatial frequency we analyze is 1/D, so the smallest separation we show is 1 λ/D.


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 ×10 improvement due to the larger aperture. This can be understood as the larger diameter causing the spatial-frequency sampling interval, 1/D2, to impose less power per λ/D, resulting in less residual power per mode.

Fig. 13

Postcoronagraph long-exposure contrast maps for a 25.4-m telescope at 800 nm, with collecting area defined by the GMT aperture. Here, we compare results for the SI and LP controllers, for 5th, 8th, and 10th magnitude stars. Note that the color scale is ×10 lower than in Fig. 11. As we do not address chromaticity here, the WFS and observation wavelengths are both 800 nm.


Fig. 14

Contrast profiles for a 25.4-m telescope with the collecting area of the GMT aperture, observing at 800 nm, using the parameters shown in Tables 2 and 3. As we do not address chromaticity here, the WFS and observation wavelengths are both 800 nm. (a) is along the main wind direction (up-and-down in Fig. 13), and (b) is perpendicular. The dashed lines are for the SI controller, standard in current AO systems. The solid lines are for the LP controller, the form of predictive control we analyze in this work. Since we are conducting a mode-by-mode analysis in discrete spatial frequencies, the lowest spatial frequency we analyze is 1/D, so the smallest separation we show is 1 λ/D.


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 λ/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 λ/D on a 0-mag star and is over 30 times better for an eighth mag star.




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.


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.


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 as57

Eq. (74)

where Fp is the photon rate from the planet and F* is the photon rate from the star (both in one spectral resolution element), and Nlines is the number of resolved lines in the planet spectrum used for crosscorrelation template matching. C is the raw PSF contrast as we have derived here. Δt is the integration time needed to reach a desired S/N. This leads directly to

Eq. (75)


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.


Code Availability

The code we developed for this analysis is freely available in a Github repository: 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:, as well as several easily available open-source libraries. No licenses are needed.



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.


Appendix A:

Mathematics of the Fourier Basis

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



The L2 norm of the cosine mode on the aperture is

Eq. (76)


For a symmetric aperture, we can restrict ourselves to n=0 with no loss of generality, which simplifies the integral to

Eq. (77)


Using the exponential form of the cosine

Eq. (78)


With the following identity:

Eq. (79)


(cf. Ref. 60 8.411 #7), the integral evaluates as

Eq. (80)


Next using the identity

Eq. (81)


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

Eq. (82)


Similarly, the L2 norm of the sine mode on the aperture is

Eq. (83)

and the norm of the product is

Eq. (84)



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)


The unobstructed circular aperture is described by Eq. (2).

Now switching to polar coordinates and collecting terms we have

Eq. (86)

Qmnc(k)=42πD20D/202π{ei2πq[(ku+mD)cos  θ+(kv+nD)sinθ]+ei2πq[(kumD)cosθ+(kvnD)sinθ]}dθqdq,
where we have made use of the complex exponential form of the cosine. We next introduce the notational simplification

Eq. (87)

and the associated angles defined by

Eq. (88)


The transform can now be written

Eq. (89)


With the following identity:

Eq. (90)

(cf. Ref. 60 8.411 #7) the transform is equivalent to

Eq. (91)

where Jn(x) are the cylindrical Bessel functions of the first kind. With the identity

Eq. (92)

(cf. Ref. 60 5.52 #1) we finally have

Eq. (93)


Similarly for the sine modes

Eq. (94)


Appendix B:

Contrast C2 of Guyon

Here we present a derivation of the contrast due to time delay and measurement error from Ref. 15, there called C2. 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)

where the coefficient hmn(t) is the combined amplitude of the Fourier mode, and ϕmn(t) is its phase. The † superscript is to differentiate this unnormalized basis from the normalized Fourier modes presented in Sec. 3.3.1. Now at a time τtl later the wavefront will be

Eq. (96)


The residual variance after applying the wavefront correction will be

Eq. (97)

which becomes

Eq. (98)


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)

which leads to the residual variance due to time lag

Eq. (100)


Now the variance of each mode will be given by the spatial PSD according to

Eq. (101)

hmn2(t)=P(kmn)D2(λ02π)2  (m2rms).

Employing the wavefront variance from measurement noise given by Eq. (52), the total residual variance at the science wavelength will be

Eq. (102)

σmn2=(λ0λ)2P(k)D2(2πkmn·V)2τtl2+βp2(kmn)S/N2(λwfsλ)2  (rad2rms).

In this open-loop framework, the equivalent time lag is τtl=τwfs+τ where τ is the closed-loop delay defined in Sec. 6.1, and τ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 τwfs, which yields a quartic equation

Eq. (103)


Each mode will have its own optimum τwfs, which is somewhat analogous to modal gains in closed-loop control. Since hmnp2=hmnc2+hmns2=hmn2, we can use the results of Sec. 4 to write

Eq. (104)


The summation in Eq. (39) defines Iϕ(r), which leads to the expression for the contrast C2

Eq. (105)


The Strehl ratio can be calculated by summing the total WFE after the convolution and then employing the Marechal approximation

Eq. (106)


It is straightforward to derive the other contrasts, C0,1,3,4,5,6, of Ref. 15 using similar arguments.

Appendix C:

Calculating the Autocorrelation

The autocorrelation R can be calculated from a numerical PSD using the following recipe:

  • 1. Start with temporal PSD Tmn(f) at N discrete frequencies fj, where f1=Δf and fN=12fs.

  • 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

    Eq. (107)


  • 3. Calculate the DFT of the augmented PSD using the FFT algorithm, giving the autocorrelation at lag τ

    Eq. (108)


    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 R.

  • 4. The lags of the points are given by

    Eq. (109)

    where only the first N+1 points are unique.


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.



G. Anglada-Escudé et al., “A terrestrial planet candidate in a temperate orbit around Proxima Centauri,” Nature, 536 437 –440 (2016). Google Scholar


C. Lovis et al., “Atmospheric characterization of Proxima b by coupling the SPHERE high-contrast imager to the ESPRESSO spectrograph,” Astron. Astrophys., 599 A16 (2017). AAEJAF 0004-6361 Google Scholar


O. Guyon et al., “How ELTs will acquire the first spectra of rocky habitable planets,” Proc. SPIE, 8447 84471X (2012). PSISDG 0277-786X Google Scholar


J. R. Males et al., “Direct imaging of exoplanets in the habitable zone with adaptive optics,” Proc. SPIE, 9148 914820 (2014). PSISDG 0277-786X Google Scholar


W. J. Borucki et al., “Kepler planet-detection mission: introduction and first results,” Science, 327 977 –980 (2010). SCIEAS 0036-8075 Google Scholar


C. D. Dressing and D. Charbonneau, “The occurrence of potentially habitable planets orbiting M dwarfs estimated from the full Kepler dataset and an empirical measurement of the detection sensitivity,” Astrophys. J., 807 45 (2015). ASJOAB 0004-637X Google Scholar


G. D. Mulders, I. Pascucci and D. Apai, “A stellar-mass-dependent drop in planet occurrence rates,” Astrophys. J., 798 112 (2015). ASJOAB 0004-637X Google Scholar


X. Bonfils et al., “A temperate exo-earth around a quiet M dwarf at 3.4 parsecs,” Astron. Astrophys., (2017). Google Scholar


J. A. Dittmann et al., “A temperate rocky super-earth transiting a nearby cool star,” Nature, 544 333 –336 (2017). Google Scholar


M. Gillon et al., “Seven temperate terrestrial planets around the nearby ultracool dwarf star TRAPPIST-1,” Nature, 542 456 –460 (2017). Google Scholar


F. J. Rigaut, J.-P. Veran and O. Lai, “Analytical model for Shack-Hartmann-based adaptive optics systems,” Proc. SPIE, 3353 1038 –1048 (1998). PSISDG 0277-786X Google Scholar


B. L. Ellerbroek, “Linear systems modeling of adaptive optics in the spatial-frequency domain,” J. Opt. Soc. Am. A, 22 310 –322 (2005). JOAOD6 0740-3232 Google Scholar


L. Jolissaint, J.-P. Véran and R. Conan, “Analytical modeling of adaptive optics: foundations of the phase spatial power spectrum approach,” J. Opt. Soc. Am. A, 23 382 –394 (2006). JOAOD6 0740-3232 Google Scholar


C. M. Correia et al., “Modeling astronomical adaptive optics performance with temporally filtered Wiener reconstruction of slope data,” J. Opt. Soc. Am. A, 34 1877 –1887 (2017). JOAOD6 0740-3232 Google Scholar


O. Guyon, “Limits of adaptive optics for high-contrast imaging,” Astrophys. J., 629 592 –614 (2005). ASJOAB 0004-637X Google Scholar


I. J. M. Crossfield, “On high-contrast characterization of nearby, short-period exoplanets with giant segmented-mirror telescopes,” Astron. Astrophys., 551 A99 (2013). AAEJAF 0004-6361 Google Scholar


C. Vérinaud et al., “Adaptive optics for high-contrast imaging: pyramid sensor versus spatially filtered Shack-Hartmann sensor,” Mon. Not. R. Astron. Soc., 357 L26 –L30 (2005). MNRAA4 0035-8711 Google Scholar


L. A. Poyneer and B. Macintosh, “Spatially filtered wave-front sensor for high-order adaptive optics,” J. Opt. Soc. Am. A, 21 810 –819 (2004). JOAOD6 0740-3232 Google Scholar


J. W. Hardy, Adaptive Optics for Astronomical Telescopes, Oxford University Press, Oxford (1998). Google Scholar


D. L. Fried, “Statistics of a geometric representation of wavefront distortion,” J. Opt. Soc. Am., 55 1427 –1431 (1965). JOSAAH 0030-3941 Google Scholar


F. Roddier et al., “One-dimensional spectra of turbulence-induced Zernike aberrations: time-delay and isoplanicity error in partial adaptive compensation,” J. Opt. Soc. Am. A, 10 957 –965 (1993). JOAOD6 0740-3232 Google Scholar


J.-M. Conan, G. Rousset and P.-Y. Madec, “Wave-front temporal spectra in high-resolution imaging through turbulence,” J. Opt. Soc. Am. A, 12 1559 –1570 (1995). JOAOD6 0740-3232 Google Scholar


R. J. Noll, “Zernike polynomials and atmospheric turbulence,” J. Opt. Soc. Am., 66 207 –211 (1976). JOSAAH 0030-3941 Google Scholar


L. A. Poyneer, B. A. Macintosh and J.-P. Véran, “Fourier transform wavefront control with adaptive prediction of the atmosphere,” J. Opt. Soc. Am. A, 24 2645 –2660 (2007). JOAOD6 0740-3232 Google Scholar


L. A. Poyneer and J.-P. Véran, “Optimal modal Fourier-transform wavefront control,” J. Opt. Soc. Am. A, 22 1515 –1526 (2005). JOAOD6 0740-3232 Google Scholar


M. Galassi et al., GNU Scientific Library Reference Manual, 3rd ed.2016). Google Scholar


O. Herscovici-Schiller, L. M. Mugnier and J.-F. Sauvage, “An analytic expression for coronagraphic imaging through turbulence. Application to on-sky coronagraphic phase diversity,” Mon. Not. R. Astron. Soc., 467 L105 –L109 (2017). MNRAA4 0035-8711 Google Scholar


A. Sivaramakrishnan et al., “Speckle decorrelation and dynamic range in speckle noise-limited imaging,” Astrophys. J. Lett., 581 L59 –L62 (2002). AJLEEY 0004-637X Google Scholar


M. D. Perrin et al., “The structure of high Strehl ratio point-spread functions,” Astrophys. J., 596 702 –712 (2003). ASJOAB 0004-637X Google Scholar


B. Macintosh et al., “Speckle lifetimes in high-contrast adaptive optics,” Proc. SPIE, 5903 59030J (2005). PSISDG 0277-786X Google Scholar


F. Malbet, J. W. Yu and M. Shao, “High-dynamic-range imaging using a deformable mirror for space coronography,” Publ. Astron. Soc. Pac., 107 386 –398 (1995). PASPAU 0004-6280 Google Scholar


C. Cavarroc et al., “Fundamental limitations on earth-like planet detection with extremely large telescopes,” Astron. Astrophys., 447 397 –403 (2006). AAEJAF 0004-6361 Google Scholar


J.-F. Sauvage et al., “Analytical expression of long-exposure adaptive-optics-corrected coronagraphic image first application to exoplanet detection,” J. Opt. Soc. Am. A, 27 A157 –A170 (2010). JOAOD6 0740-3232 Google Scholar


D. Rouan et al., “The four-quadrant phase-mask coronagraph. I. Principle,” Publ. Astron. Soc. Pac., 112 1479 –1486 (2000). PASPAU 0004-6280 Google Scholar


O. Guyon et al., “High performance Lyot and PIAA coronagraphy for arbitrarily shaped telescope apertures,” Astrophys. J., 780 171 (2014). ASJOAB 0004-637X Google Scholar


G. Prieto et al., “Giant Magellan telescope site testing seeing and turbulence statistics,” Proc. SPIE, 7733 77334O (2010). PSISDG 0277-786X Google Scholar


J. E. Thomas-Osip et al., “Giant Magellan telescope site testing summary,” (2011). Google Scholar


L. A. Poyneer et al., “Performance of the Gemini Planet Imager’s adaptive optics system,” Appl. Opt., 55 323 –340 (2016). APOPAI 0003-6935 Google Scholar


P.-Y. Madec, “Control techniques,” Adaptive Optics in Astronomy, 131 –154 Cambridge University Press, Cambridge (1999). Google Scholar


A. Oppenheim and R. Schafer, Discrete-Time Signal Processing, Pearson Education, Prentice Hall, Upper Saddle River, New Jersey (2011). Google Scholar


P. Vaidyanathan, The Theory of Linear Prediction, Morgan and Claypool, San Rafael, California (2008). Google Scholar


J. Honerkamp, Stochastic Dynamical Systems: Concepts, Numerical Methods, Data Analysis, Wiley, New York (1993). Google Scholar


W. H. Press et al., Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press, Cambridge (1992). Google Scholar


C. Dessenne, P.-Y. Madec and G. Rousset, “Optimization of a predictive controller for closed-loop adaptive optics,” Appl. Opt., 37 4623 –4633 (1998). APOPAI 0003-6935 Google Scholar


E. Gendron and P. Lena, “Astronomical adaptive optics. 1: modal control optimization,” Astron. Astrophys., 291 337 –347 (1994). AAEJAF 0004-6361 Google Scholar


C. Petit et al., “SPHERE eXtreme AO control scheme: final performance assessment and on sky validation of the first auto-tuned LQG based operational system,” Proc. SPIE, 9148 91480O (2014). PSISDG 0277-786X Google Scholar


J. R. Males et al., “The path to visible extreme adaptive optics with MagAO-2K and MagAO-X,” Proc. SPIE, 9909 990952 (2016). PSISDG 0277-786X Google Scholar


K. M. Morzinski et al., “MagAO: status and science,” Proc. SPIE, 9909 990901 (2016). PSISDG 0277-786X Google Scholar


C. Dessenne, P.-Y. Madec and G. Rousset, “Sky implementation of modal predictive control in adaptive optics,” Opt. Lett., 24 339 –341 (1999). OPLEDP 0146-9592 Google Scholar


O. Guyon and J. R. Males, “Adaptive optics predictive control with empirical orthogonal functions (EOFs),” Astron. J., (2017). Google Scholar


L. C. Johnson, D. T. Gavel and D. M. Wiberg, “Bulk wind estimation and prediction for adaptive optics control systems,” J. Opt. Soc. Am. A, 28 1566 –1577 (2011). JOAOD6 0740-3232 Google Scholar


K. Morzinski et al., “Performance of MEMS-based visible-light adaptive optics at lick observatory: closed- and open-loop control,” Proc. SPIE, 7736 77361O (2010). PSISDG 0277-786X Google Scholar


C. Petit et al., “First laboratory validation of vibration filtering with LQG control law for adaptive optics,” Opt. Express, 16 87 –97 (2008). OPEXFF 1094-4087 Google Scholar


C. Correia, J.-P. Véran and G. Herriot, “Advanced vibration suppression algorithms in adaptive optics systems,” J. Opt. Soc. Am. A, 29 185 –194 (2012). JOAOD6 0740-3232 Google Scholar


C. Kulcsár et al., “Vibration mitigation in adaptive optics control,” Proc. SPIE, 8447 84470Z (2012). PSISDG 0277-786X Google Scholar


I. Snellen et al., “Combining high-dispersion spectroscopy with high contrast imaging: probing rocky planets around our nearest neighbors,” Astron. Astrophys., 576 A59 (2015). AAEJAF 0004-6361 Google Scholar


J. Wang et al., “Observing exoplanets with high dispersion coronagraphy. I. The scientific potential of current and next-generation large ground and space telescopes,” Astron. J., 153 183 (2017). ANJOAA 0004-6256 Google Scholar


D. Mawet et al., “Observing exoplanets with high-dispersion coronagraphy. II. Demonstration of an active single-mode fiber injection unit,” Astrophys. J., 838 92 (2017). ASJOAB 0004-637X Google Scholar


I. S. Gradshteyn et al., Table of Integrals, Series, and Products, Academic Press, Oxford (2007). Google Scholar


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.

© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Jared R. Males and Olivier Guyon "Ground-based adaptive optics coronagraphic performance under closed-loop predictive control," Journal of Astronomical Telescopes, Instruments, and Systems 4(1), 019001 (6 February 2018).
Received: 1 June 2017; Accepted: 21 December 2017; Published: 6 February 2018


Back to Top