The Mysterious Lives of Speckles. I. Residual Atmospheric Speckle Lifetimes in Ground-based Coronagraphs

There have been several definitions of the atmospheric speckle lifetime presented in the literature. Early treatments of speckle dynamics were motivated by speckle interferometry, which relies on exposures short enough to "freeze the atmosphere." The optimum exposure time is therefore dependent on the correlation time of the intensity in the focal plane. A common definition for this correlation time was based on a $1/e$ reduction in the autocorrelation of the intensity, which was experimentally found to be on the order of a few ms (Scaddan & Walker 1978; Dainty et al. 1990). Roddier et al. (1982) derived the following expression for the autocorrelation e-folding time in a multi-layer Kolmogorov atmosphere:

$$\begin{eqnarray}&&{\tau }_{\mathrm{boil}}=0.36\displaystyle \frac{{r}_{0}}{\sqrt{{\sigma }_{{ \mathcal V }}^{2}}}\end{eqnarray} \tag{ 1 }$$

called there the "speckle boiling time," where r₀ is Fried's parameter and ${\sigma }_{{ \mathcal V }}^{2}$ is the variance with respect to height of the layer wind velocities ${ \mathcal V }$ above the telescope. At good astronomical sites, values of ${\tau }_{\mathrm{boil}}$ should range from ∼5 to ∼20 ms. Racine et al. (1999) employed Equation (1) for the speckle lifetime and argued that atmospheric speckles dominate photon noise by factors of 100 or more in high-contrast imaging.

The autocorrelation e-folding time is a somewhat arbitrary choice. If we instead consider the evolution of the variance in the estimate of the mean given uncorrelated trials, we can define the speckle lifetime according to the familiar "root-N" reduction in noise:

$$\begin{eqnarray}&&{\sigma }_{\mathrm{mean}}^{2}=\displaystyle \frac{{\sigma }_{o}^{2}}{t/\tau }\end{eqnarray} \tag{ 2 }$$

where τ is the noise lifetime, ${\sigma }_{o}^{2}$ is the variance of a single measurement, ${\sigma }_{\mathrm{mean}}^{2}$ is the variance in the estimate of the mean after time t. This simply states that the variance in a measurement of the mean improves with the number of uncorrelated trials $t/\tau$. Defined this way, τ is given by the integral of the autocorrelation (Fitzgerald & Graham 2006), and in this expression determines the number of uncorrelated realizations of the speckle pattern included in the mean. Note that this assumes $t\gg \tau$. Aime et al. (1986) employed this definition, and argued for longer speckle lifetimes, replacing 0.36 by  1.14 in Equation (1) based on arguments similar to Roddier et al. (1982). Using AO-corrected images, Fitzgerald & Graham (2006) empirically found significantly longer values of τ, giving an overall median of 175 ms at Lick observatory.

Macintosh et al. (2005) considered this time averaging of residual atmospheric speckles behind a coronagraph, directly analyzing the reduction in variance as a function of exposure time. They argued that the speckle lifetime should not be affected by the properties of the AO system nor by the seeing parameter r₀, and instead should be determined by the wind-crossing time of a telescope with diameter D. Based on a simulation, Macintosh et al. (2005) found that

$$\begin{eqnarray}&&\tau \approx 0.3\displaystyle \frac{D}{{ \mathcal V }}\end{eqnarray} \tag{ 3 }$$

for use in Equation (2). ⁷ This result is quite different from the previous predictions, but coarsely in agreement with the experiments by Fitzgerald & Graham (2006) (it is challenging to make more definitive experimental comparisons due to the uncertainty in atmospheric parameters).

Poyneer & Macintosh (2006) analyzed speckle noise correlations under closed-loop control, in a similar fashion to what follows. Their study found that the WFS noise and input turbulence contributed differently to speckle lifetime. They also proposed that the control-loop would impact the lifetime, through the optimization of the control law gain, and that considering long-exposure sensitivity (i.e., through Equation (2)) may lead to different optimizations than by considering raw contrast alone.

More recent experimental studies have used on-sky data with modern ExAO systems. Using focal plane data from SPHERE, Milli et al. (2016) found evidence for speckle decorrelation on several timescales, but were unable to identify an atmospheric component given the various instrument-induced timescales and the relatively low cadence used. In a study using SCExAO at much higher frame rates, Goebel et al. (2018) analyzed the evolution of the intensity using the variance of difference images (a form of autocorrelation). Here short timescale (a few ms) speckle evolution was attributed to the atmosphere, and the authors argued that the AO system modulated the speckles on longer timescales. Stangalini et al. (2017) present a study of speckle lifetimes at the LBT, where they found short-lived speckles attributed to the atmosphere and show some evidence for a difference in lifetime inside and outside the AO control radius.

Our goal in this study is to develop a self-consistent description of the time-domain behavior of speckle noise in AO-fed coronagraphic images from as close to first principles as possible. Our motivation is to understand the true limits of ground-based high-contrast imaging and so we focus on the statistical lifetimes of atmospheric speckles in coronagraphic focal planes. This requires a traversal of speckle statistics, signal processing, wavefront propagation, and control theory. To accomplish our goals we will work through the following tasks, organized by section:

• 1.  
  Section 2: Review speckle statistics to understand how atmospheric speckles dramatically reduce sensitivity.
• 2.  
  Section 3: Develop an analytic framework for analyzing power spectral densities (PSDs), which allows us to derive an expression for the statistical (mean-averaging) lifetime from a PSD.
• 3.  
  Section 4: Show how to calculate the intensity PSD in the coronagraphic focal plane behind a closed-loop AO system with a semi-analytic model.
• 4.  
  Section 5: Employ an end-to-end AO system simulation to (a) validate the statistical lifetime calculations and (b) assess the performance of the semi-analytic model.
• 5.  
  Section 6: Use the analytic framework to analyze speckle lifetime over a range of telescope diameters and atmospheric parameters.
• 6.  
  Section 7: Develop a method to analyze post-processing in the time-domain and discuss its implications for the limits of ground-based high-contrast imaging sensitivity.
• 7.  
  Section 8: Discuss various implications of these results.

The key conclusions can be gleaned from Sections 2, 6, and 8. Sections 3 through 5 and Section 7 contain the technical details.

Here we very briefly review the derivation of the intensity time-series from Males & Guyon (2018), with the goal to introduce the notation needed for this analysis. See that paper for a detailed treatment.

We use a real-valued Fourier basis

$$\begin{eqnarray}\begin{array}{rcl}{M}_{{mn}}^{p}({\boldsymbol{q}}) & = & \cos (2\pi {{\boldsymbol{k}}}_{{mn}}\cdot {\boldsymbol{q}})\\ & & +\ p\sin (2\pi {{\boldsymbol{k}}}_{{mn}}\cdot {\boldsymbol{q}})\end{array}\end{eqnarray} \tag{ 26 }$$

where $p=\pm 1$, ${\boldsymbol{q}}$ is the position vector in the pupil plane and the vector spatial-frequency is

$$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{k}}}_{{mn}}=\tfrac{m}{D}\hat{u}+\tfrac{n}{D}\hat{v}\end{array}\end{eqnarray} \tag{ 27 }$$

where m and n are integer indices, ⁸ D is the aperture diameter, and $(\hat{u},\hat{v})$ are unit vectors. The modal basis described by Equation (26) is normalized, and has the advantage that on the unobscured circular aperture both p modes at a given mn have the same symmetric spatial and temporal PSD. This is not true for pure sines and cosines. We note that this basis is the one derived for the unobscured circular aperture in a 2D Lomb–Scargle analysis (Seilmayer et al. 2020), and, as such, phase-shifted Fourier bases with the same properties exist for more complicated apertures. In this modal basis, the phase in the pupil plane is

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Phi }}({\boldsymbol{q}},t)=\displaystyle \frac{2\pi }{\lambda }\sum _{{mn}}\left[{h}_{{mn}}^{+}(t){M}_{{mn}}^{+}({\boldsymbol{q}})\right.\\ \quad \left.+\,{h}_{{mn}}^{-}(t){M}_{{mn}}^{-}({\boldsymbol{q}})\right]\end{array}\end{eqnarray} \tag{ 28 }$$

where ${h}_{{mn}}^{+}$ and ${h}_{{mn}}^{-}$ are the real-valued coefficients of the modes in the same units as wavelength, and the phase Φ has units of radians.

Males & Guyon (2018) then showed how the residual intensity due to the phase Φ behind a perfect coronagraph is described by

$$\begin{eqnarray}&&\begin{array}{l}{I}_{{\rm{\Phi }}}({\boldsymbol{r}},t)\approx {\left(\displaystyle \frac{2\pi }{\lambda }\right)}^{2}2\sum _{{mn}}\left\{\left[{\left({h}_{{mn}}^{+}(t)\right)}^{2}\right.\right.\\ \quad \left.+\,{\left({h}_{{mn}}^{-}(t)\right)}^{2}\right]\left[{{Ji}}^{2}(\pi {{Dk}}_{{mn}}^{+})\right.\\ \quad \left.\left.+\,{{Ji}}^{2}(\pi {{Dk}}_{{mn}}^{-})\right]+{\rm{cross-terms}}\right\}\end{array}\end{eqnarray} \tag{ 29 }$$

at location ${\boldsymbol{r}}={r}_{x}\hat{x}+{r}_{y}\hat{y}$ where $\hat{x}$ and $\hat{y}$ are the focal plane coordinate unit vectors. ${Ji}$ is the Jinc function

$$\begin{eqnarray}&&{Ji}(x)=\displaystyle \frac{{J}_{1}(x)}{x}.\end{eqnarray} \tag{ 30 }$$

and the modified spatial frequencies are given by

$$\begin{eqnarray}\begin{array}{rcl}{k}_{{mn}}^{+} & = & \sqrt{\left(\displaystyle \frac{{r}_{x}}{\lambda }+\displaystyle \frac{m}{D}\right)+\left(\displaystyle \frac{{r}_{y}}{\lambda }+\displaystyle \frac{n}{D}\right)}\\ {k}_{{mn}}^{-} & = & \sqrt{\left(\displaystyle \frac{{r}_{x}}{\lambda }-\displaystyle \frac{m}{D}\right)+\left(\displaystyle \frac{{r}_{y}}{\lambda }-\displaystyle \frac{n}{D}\right)}.\end{array}\end{eqnarray} \tag{ 31 }$$

Equation (29) is equivalent to a convolution with the PSF, adding contributions from nearby speckles. The cross-terms are a series of terms including products of ${h}_{{mn}}^{+}{h}_{{mn}}^{-}$ and ${h}_{{mn}}^{p}{h}_{m^{\prime} n^{\prime} }^{p^{\prime} }$ (${mn}\ne m^{\prime} n^{\prime}$) with ${Ji}(x){Ji}(x^{\prime} )$. As long as the aperture truncation is accounted for when the long-exposure statistics of h are calculated, the sum can be dropped and these small cross-terms vanish. With this caveat, the long exposure intensity is

$$\begin{eqnarray}&&\left\langle {I}_{{\rm{\Phi }}}({{\boldsymbol{r}}}_{{mn}})\right\rangle ={\left(\displaystyle \frac{2\pi }{\lambda }\right)}^{2}\left\langle | {h}_{{mn}}{| }^{2}\right\rangle .\end{eqnarray} \tag{ 32 }$$

where the discretized position

$$\begin{eqnarray}&&{{\boldsymbol{r}}}_{{mn}}={{\boldsymbol{k}}}_{{mn}}\lambda \end{eqnarray} \tag{ 33 }$$

describes the relationship between the spatial-frequency in the pupil plane and projected-angular position in the focal plane.

The quantity we seek is the temporal PSD of the intensity time-series described by Equation (29). The PSD ${{ \mathcal T }}_{{I}_{{mn}}}(f)$ is defined as in Equation (12).

The average Fourier transform of a time-series formed from the products of uncorrelated zero-mean processes will be 0, so we can employ the same logic as above and write

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal T }}_{{I}_{{mn}}}(f)\approx {\left(\displaystyle \frac{2\pi }{\lambda }\right)}^{2}\left\langle {\left|{{ \mathcal F }}_{t}\left\{{\left({h}_{{mn}}^{+}\right)}^{2}+{\left({h}_{{mn}}^{-}\right)}^{2}\right\}\right|}^{2}\right\rangle .\end{array}\end{eqnarray} \tag{ 34 }$$

We therefore directly relate the temporal PSD of I_ϕ to the temporal PSD of ${\left({h}_{{mn}}^{+}\right)}^{2}+{\left({h}_{{mn}}^{-}\right)}^{2}$. This can be estimated given the PSDs for ${h}_{{mn}}^{+}$ and ${h}_{{mn}}^{-};$ we review how these can be calculated next.

We express the PSDs of the ${h}_{{mn}}^{p}$ in closed loop (${{ \mathcal T }}_{\mathrm{cl},{mn}}(f;g)$) as functions of the open-loop (OL) PSD (${{ \mathcal T }}_{\mathrm{ol},{mn}}(f)$), the error transfer function (ETF), the measurement noise PSD (${{ \mathcal T }}_{\mathrm{ph},{mn}}(f)$), and the noise transfer function (NTF):

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal T }}_{\mathrm{cl},{mn}}(f;g)={{ \mathcal T }}_{\mathrm{ol},{mn}}(f){\left|{{\rm{ETF}}}_{\mathrm{cl}}(f;g)\right|}^{2}\\ \quad +\,{{ \mathcal T }}_{\mathrm{ph},{mn}}(f){\left|{{\rm{NTF}}}_{\mathrm{cl}}(f;g)\right|}^{2}\end{array}\end{eqnarray} \tag{ 35 }$$

where the control loop gain g is optimized to minimize the residual variance in the closed-loop PSD (Madec 1999; Poyneer et al. 2016). Males & Guyon (2018) derived OL PSDs of Fourier modes in frozen-flow turbulence. In simplified form, the WFS noise is given by

$$\begin{eqnarray}&&{\sigma }_{{ph}}^{2}=\displaystyle \frac{{\beta }_{p,{mn}}^{2}}{{F}_{\gamma }{\tau }_{\mathrm{wfs}}}\end{eqnarray} \tag{ 36 }$$

where ${\beta }_{p,{mn}}$ is the WFS sensitivity to photon noise (Guyon 2005), F_γ is the total photon rate (photons/sec) and ${\tau }_{\mathrm{wfs}}$ is the WFS integration time. This gives a white-noise PSD

$$\begin{eqnarray}&&{{ \mathcal T }}_{\mathrm{ph},{mn}}(f)=2\displaystyle \frac{{\beta }_{p,{mn}}^{2}}{{F}_{\gamma }}.\end{eqnarray} \tag{ 37 }$$

We consider ETFs and NTFs for control laws based on general linear filters of the form

$$\begin{eqnarray}&&\begin{array}{l}\widetilde{h}({t}_{i})=\sum _{j=1}^{J}{a}_{j}\widetilde{h}({t}_{i-j})+g\sum _{l=0}^{L}{b}_{l}{\rm{\Delta }}h({t}_{i-l}).\end{array}\end{eqnarray} \tag{ 38 }$$

Here the ${\widetilde{h}}_{i}$ are the estimates of the mode coefficient going back in time, and the ${\rm{\Delta }}{h}_{i}$ are the measurements of the closed-loop residual amplitude of the mode. The filter combines the J previous estimates and L previous measurements to form the optimum estimate for the current time, which is then applied as a DM command.

With J = 1, ${a}_{1}=1$, L = 0, and ${b}_{0}=1$, Equation (38) represents the simple integrator (SI). Males & Guyon (2018) showed how adding more coefficients, using the well-known Linear Prediction (LP) technique, results in significant performance improvements in terms of post-coronagraph contrast. Most, if not all, linear controllers can be cast in the form of Equation (38) (Poyneer et al. 2007; Haffert et al. 2021).

An alternate way to describe the time-evolution of the turbulent wavefront is (see Guyon 2005 and Appendix B of Males & Guyon 2018)

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{{mn}}({\boldsymbol{q}},t)=\displaystyle \frac{2\pi }{\lambda }{h}_{{mn}}^{\dagger }(t)\\ \qquad \times \ \cos \left(2\pi {{\boldsymbol{k}}}_{{mn}}\cdot {\boldsymbol{q}}+{\phi }_{{mn}}(t)\right).\end{array}\end{eqnarray} \tag{ 39 }$$

These relate to the coefficients of the basis defined by Equation (26) as

$$\begin{eqnarray}\begin{array}{rcl}{h}_{{mn}}^{+} & = & \displaystyle \frac{1}{2}{h}_{{mn}}^{\dagger }(t)\left[\cos ({\phi }_{{mn}}(t))+\sin ({\phi }_{{mn}}(t))\right]\\ {h}_{{mn}}^{-} & = & \displaystyle \frac{1}{2}{h}_{{mn}}^{\dagger }(t)\left[\cos ({\phi }_{{mn}}(t))-\sin ({\phi }_{{mn}}(t))\right].\end{array}\end{eqnarray} \tag{ 40 }$$

Under the Taylor frozen flow hypothesis the amplitude ${h}_{{mn}}^{\dagger }$ of each spatial frequency is approximately constant over short periods, which means that ${h}_{{mn}}^{-}({\phi }_{{mn}}(t))$ = ${h}_{{mn}}^{+}({\phi }_{{mn}}(t)+\tfrac{\pi }{2})$. This phase shift yields

$$\begin{eqnarray}&&{{ \mathcal F }}_{t}\left\{{h}_{{mn}}^{-}\right\}={{ \mathcal F }}_{t}\left\{{h}_{{mn}}^{+}\right\}{e}^{i\tfrac{\pi }{2}}.\end{eqnarray} \tag{ 41 }$$

This temporal-phase relationship between the two modes at each spatial frequency must be included in the PSD calculation in order to obtain accurate speckle lifetimes.

In the following recipe for generating PSDs, we will produce time-series of the modal coefficients ${h}_{{mn}}^{p}$. In addition to the temporal phase relationship, there is a modest spatial correlation between the modes, even under good correction. This is illustrated in Figures 2 and 3 of Males & Guyon (2018). This can be accounted for using the covariance matrix ${\boldsymbol{\Sigma }}$ of the ${h}_{{mn}}^{p}$, which is calculated as described in that paper. Given an initially uncorrelated vector ${\tilde{h}}_{{mn}}^{p}$, and the decomposition ${\boldsymbol{\Sigma }}={{\boldsymbol{AA}}}^{T}$, the spatially correlated vector of coefficients is found from

$$\begin{eqnarray}&&{h}_{{mn}}^{p}={\boldsymbol{A}}{\tilde{h}}_{{mn}}^{p}.\end{eqnarray} \tag{ 42 }$$

Zoom In Zoom Out Reset image size

We note that this adds significant computational cost to the following calculations, for relatively minor benefits. In practice we have found that it only modestly improves the match between these calculations and simulations at higher temporal frequencies, and has little to no noticeable impact on speckle lifetime predictions. We include it here for completeness, and consider it optional depending on the goals for employing the model.

We are now ready to calculate ${{ \mathcal T }}_{{I}_{{mn}}}(f)$. Unfortunately, we have found no purely analytic way to calculate $\langle {\left|{{ \mathcal F }}_{t}\left\{{h}^{2}\right\}\right|}^{2}\rangle$ given only $\langle {\left|{{ \mathcal F }}_{t}\left\{h\right\}\right|}^{2}\rangle$. This is further complicated by the spatial correlations and temporal relationship between modes. We therefore employ a Monte Carlo procedure:

• 1.  
  Generate an OL time-series for each ${h}_{{mn}}^{p}$ for both parities. Note that ${h}_{{mn}}^{-}$ is needed for the optional spatial correlation step, even though it is discarded eventually.

  • (a)  
    Generate a Gaussian white-noise time-series of the desired length.
  • (b)  
    Fourier transform the white-noise time-series.
  • (c)  
    Multiply the Fourier transform by the square-root of the OL PSD. We now have ${{ \mathcal F }}_{t}\left\{{h}_{{mn}}^{p}\right\}$ (see Kasdin 1995).
  • (d)  
    Apply the inverse Fourier transform to produce ${h}_{{mn}}^{p}(t)$. If the correlation matrix will be applied in Step 2, then normalize to unit variance.
  • (e)  
    Repeat for each mode.
• 2.  
  Optionally, correlate the modal time-series by multiplying the vector of coefficients at each time-step by the decomposition of the covariance matrix.
• 3.  
  For each spatial frequency mn, replace ${h}_{{mn}}^{-}(t)$ with the phase shifted ${h}_{{mn}}^{+}(t)$.
• 4.  
  Generate a WFS measurement noise time-series for each mode, ${\eta }_{{mn}}^{+}$ and ${\eta }_{{mn}}^{-}$, and Fourier transform them. These are white-noise PSDs, with no phase relationship.
• 5.  
  Generate the closed-loop Fourier transforms for each mode

  $$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal F }}_{t}\{{h}_{{mn},{cl}}^{+}\}={{ \mathcal F }}_{t}\{{h}_{{mn}}^{+}\}\mathrm{ETF}(f)\\ \quad +\ {{ \mathcal F }}_{t}\{{\eta }_{{mn}}^{+}\}\mathrm{NTF}(f)\\ \quad {{ \mathcal F }}_{t}\{{h}_{{mn},{cl}}^{-}\}={{ \mathcal F }}_{t}\{{h}_{{mn}}^{-}\}\mathrm{ETF}(f)\\ \quad +\ {{ \mathcal F }}_{t}\{{\eta }_{{mn}}^{-}\}\mathrm{NTF}(f).\end{array}\end{eqnarray} \tag{ 43 }$$
• 6.  
  Calculate the closed-loop modal coefficient time-series with the inverse Fourier transform. We now have ${h}_{{mn},{cl}}^{+}(t)$ and ${h}_{{mn},{cl}}^{-}(t)$
• 7.  
  Now we can calculate the focal plane intensity as a function of time

  $$\begin{eqnarray}&&{I}_{{mn}}(t)\approx \displaystyle \frac{2\pi }{\lambda }\left[{\left({h}_{{mn},{cl}}^{+}\right)}^{2}+{\left({h}_{{mn},{cl}}^{-}\right)}^{2}\right].\end{eqnarray} \tag{ 44 }$$
• 8.  
  Finally, we calculate $| {{ \mathcal F }}_{t}\{{I}_{{mn}}(t)\}{| }^{2}$. For best results, this is done with a window function applied.
• 9.  
  Now the entire procedure is repeated many times, averaging the result to form the estimate of ${{ \mathcal T }}_{{I}_{{mn}}}(f)$.

In the following sections, we compare the results of these calculations to end-to-end numerical simulations and analyze the results in terms of speckle lifetimes.

We conducted a series of end-to-end numerical simulations of an idealized instrument to validate our calculations, using MagAO-X (Males et al. 2020) as an example. Our goal is not to develop a realistic hardware simulation, rather to ensure that the calculations produce reasonable results.

We simulated a circular unobscured aperture of 6.5 m. An idealized deformable mirror was used, which exactly produced commanded shapes as projection of modes. The DM was capable of making 2400 Fourier modes, corresponding to the Nyquist limit of a 48 × 48 actuator illuminated pupil. The wavefront sensor (WFS) worked by simply projecting modes onto the phase screen. The basis set used for control was the modified Fourier basis (Equation (26)) but orthonormalized with the stabilized Gramm–Schmidt procedure. This was necessary for the loop to be stable under the LP controller. Photon noise was added as Gaussian noise uniformly distributed across the WFS image. The noise was scaled to correspond to the unmodulated pyramid sensor (${\beta }_{p}=\sqrt{2}$, Guyon 2005).

Wavefront propagation was monochromatic using the Fraunhofer approximation, and the science and sensing wavelength were both 800 nm. The WFS was spatially filtered in the Fourier domain to minimize aliasing. The perfect coronagraph was simulated (Cavarroc et al. 2006). The wavefront was spatially filtered before the coronagraph as well, to minimize Gibbs ringing on the dark hole edge due to the digitized pupil.

The simulated system operated at 2000 Hz with a 1.5 frame delay in addition to the 1 frame total sample and hold for WFS integration. To achieve the 1/2 frame, the WFS used the average turbulence between two time-steps. Turbulence was generated in 7 layers as fixed oversized phase screens, translated according to the layer wind speed. Turbulence parameters corresponded to the median LCO model used in Males & Guyon (2018), where ${r}_{0}=0.16$ m, ${L}_{0}=25$ m, and layer averaged wind speed was 18.7 m s⁻¹. The control laws, both SI and LP with 21 coefficients, were optimized on the OL PSDs measured by running the simulation with the control loop off.

Simulations were conducted first for an infinitely bright star, that is without photon noise, to analyze the impact of only the input disturbance dynamics. Next, an 8th magnitude star was simulated to show the effects of photon noise in the WFS. For both cases, simulations with both the SI and LP controllers were conducted. Individual simulation runs covered an elapsed time of 65 s, and the final 60 s were analyzed to be sure that the control system had stabilized. Each of these four cases were run 10 times with a new turbulence phase screen, and the results averaged.

Figure 2 shows post-coronagraph contrast maps, comparing the results of calculations using Males & Guyon (2018) to the simulation outputs. We show results for an infinitely bright ($-\infty$ mag) guide star to show the effects of the input turbulence alone, and for an 8th mag guide star. Figure 3 shows the median radial profile for each of the maps in Figure 2. The agreement is very good, though some small differences are visible in the LP (bottom row). These differences, apparent in both the maps and the profile, are due to the LP controller optimization leading to slightly different controller characteristics. This is due to the difference in the modal basis actually under control in the simulations compared to the model. Another possible source of slight differences is spatial-temporal correlation not fully accounted for in the semi-analytic calculations, which may become apparent at the correction levels under predictive control.

Zoom In Zoom Out Reset image size

Figure 4 evaluates the dynamics of a single spatial frequency component $m,n=10,14$. Note that this required re-fitting the Fourier coefficient after the simulations completed, since the control basis was different. For consistency and brevity we use only this mode for illustration of single-mode dynamics. It was arbitrarily chosen, and is representative of typical results.

Zoom In Zoom Out Reset image size

In the top row, the black lines compare the calculated (dotted–dashed) and simulated (solid) OL PSD. The blue lines compare the output of the calculations to the simulations for the SI, and red lines show the LP. The bottom row shows the ETF for each controller. On the $-\infty$ mag guide star, the OL PSD shows excellent agreement at low frequencies, with a slight divergence in slope beginning at a ∼100 Hz. The SI likewise shows nearly perfect agreement, with the exception of a lower overshoot peak in the simulation. The LP calculations are qualitatively a good match to the simulations, showing the same general shape, but the simulation shows better rejection at lower temporal frequencies. This better low-frequency performance is true on the 8th mag star for both SI and LP. The difference in integrated power at low frequencies is small and does not appreciably impact final contrast, and as we show in the following section leads to good results in intensity.

The disagreements evident in the figure are due to the differences between the idealized continuous-time semi-analytic calculations and the discretized simulations. In particular the empirical optimization of the control laws produces different dynamical behavior in the simulations compared to the model predictions Additionally, the use of an orthogonalized basis in the simulation produces different noise rejection characteristics on the 8th mag star.

Overall the comparison of outputs presented here give us confidence that the semi-analytic model gives reasonable results, and can be used to analyze the closed-loop contrast performance and wavefront dynamics in a coronagraphic AO system.

In Figure 5 we compare the intensity PSDs in the post-coronagraph focal plane. In intensity, the agreement between calculations and simulations is even better than in phase. There is a departure at high frequencies, where the calculations predict stronger peaks than occur in the simulations. The agreements in overall power and the fraction of power contained in low-frequencies, which are most important for determining speckle lifetimes, are very good.

Zoom In Zoom Out Reset image size

In Figure 6 we show the evolution of the variance of the mean with time. The lines show the result of applying Equation (20) to PSDs of the intensity in the simulations. To test the validity of the derivation of speckle lifetime, the points show the direct calculation of the variance in the simulations, as a function of increasing averaging length. Importantly, the solid lines are in nearly perfect agreement with the points, showing the validity of Equation (20). The dotted–dashed lines, which are the output of the semi-analytic calculations, are slightly offset for the LP, which is due to differences in simulated versus predicted contrast (equivalent to the overall variance level) as found in the previous section, but the shape of the predicted evolution with time for both SI and LP controllers is in good agreement with the simulations.

Zoom In Zoom Out Reset image size

Finally in Figure 7 we present maps of speckle lifetime, and in Figure 8 the median radial profile of speckle lifetime. The calculated and simulated lifetimes agree throughout most of the control region. The noteworthy exception occurs in the LP roughly along a diagonal. At wider separations the simulations predict lifetimes ∼50% higher in a wedged shape region. As with the differences in the phase dynamics discussed above, we attribute this to differences between the orthogonalized modal basis used in the simulations and the Fourier modes used in the calcualtions, and the resultant changes in the output of the controller coefficients and gain optimization.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For comparison, the "crossing-time" speckle lifetime of Macintosh et al. (2005) would be $\tau =0.3D/{ \mathcal V }=224$ ms for the turbulence averaged windspeed used in these calculations, with the caveat that the constant of propotionality is expected to vary with coronagraph. We find lifetimes much shorter than this, up to a factor of 10 shorter for the LP controller. Figures 7 and 8 show that speckle lifetime has a strong dependence on control law, and varies with position in the image. We investigate the $D/{ \mathcal V }$ scaling below.

Overall these results show that Equation (25) correctly describes the statistical speckle lifetime for a given intensity PSD, which in turn can be predicted by our semi-analytic model of an ExAO coronagraph.

The results shown here are subject to several caveats. All of the assumptions used in the development of Males & Guyon (2018) apply here. In particular, our analysis is monochromatic. While we assumed that the science and sensing wavelength were the same, neither component of an AO system can be truly monochromatic. Though speckle lifetime formally does not depend on wavelength, which only scales the PSD without changing its shape, AO performance will depend on how light is distributed in the instrument. Furthermore, atmospheric chromaticity will cause differential errors across the bandpass of the WFS, and between the science and WFS bandpasses. These errors will contribute speckles in addition to those explicitly treated here, which will have different dynamics as they are not subject to the same control filters. These speckles can, however, be analyzed with the same tools we developed here. We save this for future work.

A key assumption underlying the mathematics of the PSD analysis is that the speckle intensity is WSS, which means that the mean is constant and that the autocorrelation does not depend on time. As discussed above, this assumption implies that r₀ and the layer ${\boldsymbol{ \mathcal V }}$ are constant such that the statistics of the AO control output process are WSS. In the simulations presented these assumptions are true by construction, and so the comparison to our semi-analytic predictions is valid. On-sky input turbulence parameters will vary, and so the process will not be WSS over long periods. However, these variations occur over finite times. For our results to be useful the timescales of these variations would need to be long enough to estimate a PSD which meaningfully captures the intensity dynamics. What this ultimately means is that we must consider the speckle lifetime as itself a time-variable parameter, but we can expect there to be a representative value over an observation that characterizes long exposure sensitivity.

We have developed a semi-analytic model of a coronagraphic AO system that allows us to estimate long-exposure contrast, the intensity dynamics, and statistical lifetimes. Compared to an idealized end-to-end simulation, the results of the model compare favorably. This powerful technique allows us to analyze the impact of control-law choice, as well as to compare the influence of various telescope and atmosphere parameters on the speckle lifetime. Key benefits of such modeling are that the true long-exposure performance of a system can be determined without the expense of simulations. For instance, the parameter study in Section 6 would have been prohibitively expensive if conducted with simulations. As design efforts continue for the coming GSMTs, such models will be useful for similar parameter studies to optimize high contrast imaging instruments.

We analyzed various parameters with potential impacts on speckle lifetime. Overall, the simplest summary of this analysis is that there is no single thumb-rule-type relationship that applies to all parameters. This modestly contradicts the conclusions of Macintosh et al. (2005). There, the AO system was modeled as a spatial filter, with the result that while the speckle amplitude was reduced, the evolution of correlations depended only on the crossing-time. As a result, neither the control action nor the overall power (seeing) should impact speckle lifetime. Only the wind would change speckle lifetime for a given pupil diameter.

Here we have added the time dimension, fully considering the action of AO control as a temporal filter. Furthermore, we have modeled the impact of high-order predictive control. With this more complete view, we find that the control action does reduce speckle lifetime. There is an influence of the crossing-time as predicted by Macintosh et al. (2005) to some extent, but in addition to the reduction in amplitude, the time-domain influence of the control system is to increase the rate at which one speckle turns into the next. That is, the control reduces the correlation length. This is consistent with the findings of Poyneer & Macintosh (2006), which predicted that control laws could be optimized based on speckle lifetime to maximize sensitivity.

Predictive control, already expected to improve contrast (Poyneer et al. 2007; Correia et al. 2017; Males & Guyon 2018; Haffert et al. 2021), will further improve sensitivity by decreasing the speckle averaging lifetime. On bright stars, predictive control nearly eliminates dependence on the crossing-time as it acts to reduce servo-lag errors.

As ExAO coronagraphs approach the residual atmosphere limit, it will be increasingly important to understand the impact of correlated noise from atmospheric speckles on sensitivity. Our results demonstrate that this is a complicated task, that must consider all relevant parameters of the system under study.

The temporal decorrelation of speckle noise is one of the main benefits of PSF-subtraction post-processing (Marois et al. 2006, 2008), even though algorithms have typically focused on exploiting spatial correlations. Recently, studies have begun to consider the time-domain in post-processing explicitly (Samland et al. 2021), and several algorithms show promise in exploiting short exposure focal-plane data (Walter et al. 2019; Steiger et al. 2021) and WFS telemetry (Rodack et al. 2021, submitted) (Frazin & Rodack 2021, submitted) to estimate speckle intensity. Knowledge of the intensity PSD allows for such algorithms to be analyzed efficiently. We demonstrate this by considering the prospect of telemetry-based reconstruction of the speckle time-series. The algorithm we outlined made use of the main WFS images, but other information sources such as accelerometers, low-order WFS telemetry, and environmental sensors could be used. The key challenge is in the construction of the system propagation model. We used the temporal intensity PSDs to predict the sensitivity gain possible if this can be done. If the first few low-order temporal modes, here described in terms of the Legendre polynomials, can be determined and subtracted from the intensity in the focal plane, large improvements in exposure time are possible. The fact that only a few low-order temporal modes need to be reconstructed has the important consequence that any such post-processing algorithm does not need to accurately determine the intensity at every time interval, and does not need to be perfect. Accuracy in the mean, slope, and other low-order terms of a few percent is all that is needed to lower the speckle lifetime and therefore increase sensitivity.

The Fourier transform of the Legendre polynomials is ⁹

$$\begin{eqnarray}&&{\int }_{-1}^{1}{{\mathscr{P}}}_{\eta }(x){e}^{i2\pi {kx}}{dx}={i}^{\eta }\displaystyle \frac{{J}_{\eta +\tfrac{1}{2}}(2\pi k)}{\sqrt{k}}\end{eqnarray} \tag{ A1 }$$

with J_η denoting the Bessel functions of the first kind. We can then write the Fourier transforms of the orthonormal ${P}_{\eta }(x)$ as

$$\begin{eqnarray}&&{Q}_{\eta }(k)={i}^{\eta }\sqrt{\displaystyle \frac{2\eta +1}{2}}\,\displaystyle \frac{{J}_{\eta +\tfrac{1}{2}}(2\pi k)}{\sqrt{k}}.\end{eqnarray} \tag{ A2 }$$

For analysis of Legendre mode filtering these have the useful property that ¹⁰

$$\begin{eqnarray}&&{\int }_{0}^{\infty }{\left|{Q}_{\eta }(k)\right|}^{2}{dk}=1{\rm{for}}\,{\rm{all}}\,\eta .\end{eqnarray} \tag{ A3 }$$

The expansion of I(t) in terms of the ${P}_{\eta }(\xi )$ is given by Equations (17) and (18). This series expansion in terms of normalized Legendre polynomials is the 1D analog of the expansion in Zernike polynomials used by Noll (1976). Following that analysis we write the covariance between any two ${P}_{\eta }(\xi )$ as

$$\begin{eqnarray}&&\begin{array}{l}\langle {a}_{\eta }^{* }{a}_{\eta ^{\prime} }\rangle =\displaystyle \int \displaystyle \int {P}_{\eta }(\xi )C\left(\displaystyle \frac{T}{2}\xi ,\displaystyle \frac{T}{2}\xi ^{\prime} \right)\\ \quad \times \,{P}_{\eta ^{\prime} }(\xi )d\xi d\xi ^{\prime} \end{array}\end{eqnarray} \tag{ A4 }$$

the autocorrelation is given by

$$\begin{eqnarray}&&C\left(\displaystyle \frac{T}{2}\xi ,\displaystyle \frac{T}{2}\xi ^{\prime} \right)=\left\langle I\left(\displaystyle \frac{T}{2}\xi \right)I\left(\displaystyle \frac{T}{2}\xi ^{\prime} \right)\right\rangle .\end{eqnarray} \tag{ A5 }$$

Equivalently, in the Fourier domain we have

$$\begin{eqnarray}&&\begin{array}{l}\langle {a}_{\eta }^{* }{a}_{\eta ^{\prime} }\rangle =\displaystyle \int \displaystyle \int {Q}_{\eta }^{* }(\kappa ){\rm{\Phi }}\left(\displaystyle \frac{2}{T}\kappa ,\displaystyle \frac{2}{T}\kappa ^{\prime} \right)\\ \quad {Q}_{\eta ^{\prime} }(\kappa )d\kappa d\kappa ^{\prime} \end{array}\end{eqnarray} \tag{ A6 }$$

where for a WSS process

$$\begin{eqnarray}&&{\rm{\Phi }}\left(\displaystyle \frac{2}{T}\kappa ,\displaystyle \frac{2}{T}\kappa ^{\prime} \right)=\displaystyle \frac{2}{T}{ \mathcal T }\,\left(\displaystyle \frac{2}{T}\kappa \right)\delta (\kappa -\kappa ^{\prime} )\end{eqnarray} \tag{ A7 }$$

after the change of variables $f\to (2/T)\kappa$ in the PSD ${ \mathcal T }(f)$.

We now have the covariance of the Legendre coefficients

$$\begin{eqnarray}&&\begin{array}{l}\langle {a}_{\eta }^{* }{a}_{\eta ^{\prime} }\rangle ={i}^{\eta ^{\prime} -\eta }\displaystyle \frac{\sqrt{2\eta +1}\sqrt{2\eta ^{\prime} +1}}{T}\\ \qquad \times \ {\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{J}_{\eta +\tfrac{1}{2}}(2\pi \kappa ){J}_{\eta ^{\prime} +\tfrac{1}{2}}(2\pi \kappa )}{\kappa }{ \mathcal T }\,\left(\displaystyle \frac{2}{T}\kappa \right)d\kappa \end{array}\end{eqnarray} \tag{ A8 }$$

and for $\eta =\eta ^{\prime}$, the variance in a single Legendre mode coefficient over the sample length T is

$$\begin{eqnarray}&&\left\langle {\left|{a}_{\eta }\right|}^{2}\right\rangle =\displaystyle \frac{2\eta +1}{T}{\int }_{0}^{\infty }\displaystyle \frac{{J}_{\eta +\tfrac{1}{2}}^{2}(2\pi \kappa )}{\kappa }{ \mathcal T }\,\left(\displaystyle \frac{2}{T}\kappa \right)d\kappa .\end{eqnarray} \tag{ A9 }$$