Measurements of Speckle Lifetimes in Near-infrared Extreme Adaptive Optics Images for Optimizing Focal Plane Wavefront Control

Unlike speckles subtracted away during post-observation image processing, speckles interferometrically destroyed during observations do not leave behind photon noise. Because a speckle is fundamentally scattered starlight, it is possible to apply equal-amplitude but out-of-phase starlight on top of it and thereby destructively interfere it away. On the other hand, light from a substellar companion is not coherent with the starlight, so it cannot be destructively nulled away.

A number of techniques have been proposed to differentiate between speckles caused by scattered starlight and substellar companions and then mitigate the speckles. The pioneering work in real-time speckle-nulling algorithms and conceptually most straightforward method was proposed by Bordé & Traub (2006). The focal plane, which is where the science detector is located, shows the Fourier transform of the pupil plane, which is where the deformable mirror is located. A sine wave on the DM generates delta functions (PSFs) on the science detector. This is shown in Figure 1. The amplitude, phase, orientation, and frequency of the sine wave on the DM determines the brightness, phase, location, and spacing of the artificial speckles on the science detector. Therefore, one can generate an artificial speckle in the same location as an existing speckle, scan through phase space to find the opposite complex amplitude, and thereby null it away. The upside of this technique is that it is robust and does not require a detailed system model. The primary downside is that several phase measurements need to be made before a speckle can be nulled, so it is not as fast as other algorithms. Martinache et al. (2014) implemented this speckle-nulling algorithm on SCExAO. In practice, it worked well on the SCExAO internal light source for quasi-static speckles, but because the camera used operated at ∼170 Hz and was noisy, and the loop required many phase measurements per iteration, it was not fast enough to affect the dynamic speckles when tested on-sky.

Zoom In Zoom Out Reset image size

A more complex but potentially faster speckle-nulling method known as "electric field conjugation" (EFC) was proposed by Give'on et al. (2007). This method places small displacements on the deformable mirror to derive the electric field and thereby phase of speckles. It has produced the deepest contrasts (on the order of 10⁻⁹) in highly stable high-contrast imaging laboratories, but it requires a detailed system model and high S/N because the actuator pokes are quite small. The theory of this technique forms the basis of several other proposed speckle-nulling ideas. Other techniques for real-time speckle discrimination include coronagraphic focal plane wavefront estimator for exoplanet imaging (COFFEE; Sauvage et al. 2012; Herscovici-Schiller et al. 2018), the self-coherent camera (Baudoz et al. 2006), the phase-shifting interferometer (Serabyn et al. 2011; Bottom et al. 2017), synchronous interferometric speckle subtraction (Guyon 2004), and linear dark field control (Miller et al. 2017). However, apart from the implementation by Martinache et al. (2014) of the technique of Bordé & Traub (2006), none of these methods have been used for closed-loop speckle-nulling on-sky.⁹ This is primarily because they are optomechanically complex, require a higher degree of PSF stability than can be achieved on-sky, require very high Strehl ratios, or require deformable mirror perturbations that would be detrimental to science observations. Jovanovic et al. (2018) provides an informative review of these techniques.

Regardless of the particular algorithm, our speckle lifetime measurements are critical to quantify closed-loop speckle-nulling performance. As loop bandwidth increases, shorter-lived speckles can be destroyed, and therefore the focal plane contrast and detection limits improve. An understanding of the rates at which speckles change is important because it influences the feasibility of implementing the loops described above and gives the bandwidth necessary to reach a given performance level.

Our goal is to measure how speckles evolve as a function of time and brightness. We quantified this using a new technique described below. The single largest detriment to image quality at SCExAO is tip/tilt caused by vibrations from the telescope and instrument (Lozi et al. 2016), which causes translation of the speckles instead of evolution, so it needed to be mitigated. Therefore, we began by aligning all the images. The results of the image alignment are shown in Figure 2. We located the artificial astrometric speckles and aligned images using those, and as a check also cross-correlated each image against the next to detect shifts. Both techniques registered images to ∼0.01 pixel accuracy and produced similar results.

Zoom In Zoom Out Reset image size

Second, we selected all pixels within a radius of 3λ/D <r < 10λ/D (0.12'' < r < 0.40'') of the PSF core. Given that the observed intensity I of a speckle is the square of the modulus of complex amplitude A,

$$\begin{eqnarray}&&I={| {\bf{A}}| }^{2},\end{eqnarray} \tag{ 1 }$$

it follows that

$$\begin{eqnarray}&&\displaystyle \frac{{dI}}{{dt}}=2{\bf{A}}\displaystyle \frac{d{\bf{A}}}{{dt}}.\end{eqnarray} \tag{ 2 }$$

A is a random variable with a characteristic evolutionary timescale (temporal derivative). This is due to linearity between A and turbulence in the pupil, which to first degree, over adequately short timescales (a few ms), has a linear temporal evolution (fixed derivative). On these timescales, the temporal evolution of phase in the pupil is dominated by translation of turbulence layers due to wind. Because we are considering speckles at small angles (3–10λ/D), these correspond to low spatial frequencies on the pupil (period = 2.6–0.8 m). During a few milliseconds, a 10 m/s wind (which is higher than what we observed) translates a turbulence screen by a few cm. This is small relative to the aforementioned spatial frequencies in the pupil. For this reason, temporal evolution of phase can be linearized. There is linearity between phase in the pupil and complex amplitude in the focal plane, as the focal plane shows the Fourier transform of the pupil plane, and Fourier transforms are linear. Therefore, on millisecond timescales, the complex amplitude at the focal plane will evolve linearly.

For this reason, dA/dt is static over short timescales, and brighter (larger I) speckles will change more quickly than dimmer ones. Because of this, we sorted the pixels into four brightness bins (we selected the 20–30 percentile brightness pixels, 40–50 percentile, etc.). For each selection of pixels, we subtracted those pixels from the same pixels time t later. This formed a difference image. Finally, we calculated the standard deviation of these difference images to measure how the pixels changed in brightness over that time interval.¹⁰ If a pixel did not change in brightness, this quantity would be 0. We divided the standard deviation by the peak brightness of the PSF core so that it had units of contrast. We repeated this process for each image compared with 1–260 frames after it, and then averaged all the measurements at a given temporal separation to improve the S/N. We plotted this quantity against time for each brightness bin in Figure 3 (extreme AO, AO188 only, and no AO on the night of 2017 May 31) and 4 (extreme AO on 2017 August 13 and 15).

Zoom In Zoom Out Reset image size

To understand whether the behavior of the standard deviations illustrated in Figures 3 and 4 was in fact due to speckle behavior or merely instrumental effects, we observed the SCExAO internal light source. There were no moving components in this setup. As is illustrated in Figure 5, indeed that speckle behavior is entirely static over the timescales examined. Second, we ran unilluminated frames through the above analysis pipeline to see how noise affected the plot's behavior. Again, the standard deviations were static in time. This means that the behavior in Figures 3 and 4 was indeed caused by the temporal evolution of the on-sky speckles.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In all cases, the brighter speckles change in brightness more rapidly than the dimmer ones. This confirms the behavior we expected from Equation (2). The implication of this is that once a dark hole is present in an image, it is easier for a speckle-nulling loop to keep it dark. Because a dimmer speckle changes more slowly, more speckle-nulling loop iterations can be carried out before it changes. Therefore, the performance of a real-time speckle-nulling loop will be highly affected by the Strehl ratio and/or presence of a coronagraph that creates dark holes in the image.

On sufficiently short timescales, the complex amplitude can be approximated as behaving linearly with time, i.e.,

$$\begin{eqnarray}&&{\bf{A}}(t)\approx {at}+{a}_{0}+{bti}+{b}_{0}i\end{eqnarray} \tag{ 3 }$$

where a and b are the rates at which the real and imaginary parts of A change, respectively, and a₀ and b₀ are their initial values. Inserting this into Equation (1) and computing the measured intensity as a function of time I(t) produces

$$\begin{eqnarray}&&I(t)=({a}^{2}+{b}^{2}){t}^{2}+(2{{aa}}_{0}+2{{bb}}_{0})t+({a}_{0}^{2}+{b}_{0}^{2}).\end{eqnarray} \tag{ 4 }$$

Therefore, the measured intensity of a speckle is a quadratic function of time, and its behavior depends on the complex amplitude's initial value and temporal rate of change. The standard deviation of pixel intensity, which is what we plotted in Figures 3–5, behaves the same way. For example, if a speckle is entirely static, then a = b = 0, and $I(t)={a}_{0}^{2}+{b}_{0}^{2}$, and this produces a flat line in the plot with a y-intercept which depends on the brightness the speckle. This exact case is illustrated in the top panel of Figure 5. On the other hand, the other simple case occurs when the speckles are changing (i.e., $a\ne 0$ and $b\ne 0$) but start from zero intensity (a₀ = b₀ = 0). In this case, the plot would exhibit an upward parabola centered at the origin. We approach this limit (but do not reach it) in the on-sky observations by selecting increasingly dim speckles (i.e., a₀ → 0 and b₀ → 0) for analysis.

Depending on the temporal sampling of the data, the initial brightness of the speckles, and the rates at which they change, this upward parabolic behavior may not be obvious in real-world measurements. The first few data points in the May 31 ExAO plot appear linear, the first few points in the August 13 plot appear downwardly parabolic, and only in the August 15 plot are the first few points plausibly upwardly parabolic. It is possible that our data are not fast enough to recover this predicted effect, and it would become apparent in data with even faster temporal sampling. The upward parabolic behavior assumes that the complex amplitude at a pixel changes linearly with time, and it certainly departs from this after a few milliseconds. Additionally, at sufficiently long timescales (≳0.04 s), the speckles have entirely changed from their initial state, causing the measured standard deviations to plateau. This is why our plots approach asymptotic standard deviations at large temporal separations.

If the behavior described by Equation (4) was evident in the images, one could predict the performance of real-time speckle-nulling loops from it. Because it is unclear whether our data reveal the predicted behavior, we will not do that, but we will explain how it could be done. If our analysis described above could be carried out at even higher temporal sampling or in an extremely high Strehl case, and the upward parabolic behavior revealed itself, one could fit a second-order polynomial to it. Solving this quadratic equation would provide the contrast that could be achieved at a given speckle-nulling loop bandwidth. The y-intercept in this fit would correspond to the read and photon noise (i.e., the noise that would be present if it was possible to image repeatedly with zero time between the images and with zero speckle evolution). A real-time speckle-nulling loop would zero the linear and constant terms with each update by forcing a₀ = b₀ = 0. In this case, I(t) ∝ t². Therefore, the achievable contrast would increase with the square of the speckle-nulling loop bandwidth (where bandwidth is 1/t). This result emphasizes the major impact that speckle-nulling loops will have on detection limits.

Milli et al. (2016) employed a different technique for measuring the timescales over which speckles evolve. Their observations were comparable to ours, so we analyzed our comparatively higher frame rate data using their methods. They collected H-band images with a ∼55 nm bandpass at a frame rate of 1.6 Hz of coronagraphic PSF images from the SPHERE ExAO instrument. They then quantified the change in speckles over the 52 minutes of their on-sky observations by computing Pearson's correlation coefficient (Pearson 1895) for pairs of images. Pearson's correlation coefficient ρ(t_i, t_j) between two frames collected at times t_i and t_j can be defined as follows:

$$\begin{eqnarray}&&\rho ({t}_{i},{t}_{j})=\displaystyle \frac{{\displaystyle \sum }_{x\in S}[I(x,{t}_{i})-\overline{{I}_{{t}_{i}}}][I(x,{t}_{j})-\overline{{I}_{{t}_{j}}}]}{{\sigma }_{{t}_{i}}{\sigma }_{{t}_{j}}{N}_{\mathrm{pix}}}\end{eqnarray} \tag{ 5 }$$

where the sum occurs over all pixels in the region S, I is the instantaneous intensity at a given pixel with position x, $\overline{I}$ is the instantaneous average intensity over the region, σ is the spatial standard deviation over the region, and N_pix is the number of pixels in the region. The normalization is chosen so that ρ = 1 when t_i = t_j (equivalently, if the speckles did not change at all between images at t_i and t_j, ρ = 1). Milli et al. (2016) demonstrated that this quantity is directly related to the contrast obtained after subtracting one image of the pair from the other. The S/N of ρ for a single pair of frames is typically quite poor, so they improved it by averaging the values of ρ computed for many frame pairs separated by equal amounts of time.

After completing their on-sky observations, Milli et al. (2016) observed the instrument's internal light source and repeated the analysis to separate atmospheric effects from instrumental effects. In both cases, they observed two primary temporal regimes of speckle decorrelation. Over the first few seconds of their observations, the speckles decorrelated rapidly with an exponential decay. They fitted the first 30 seconds of their measurements with the equation

$$\begin{eqnarray}&&\rho (t)={\rm{\Lambda }}{e}^{-t/\tau }+{\rho }_{0}\end{eqnarray} \tag{ 6 }$$

and found a best-fitting decay time of τ = 3.5 s for their on-sky observations. They initially posited that these speckles were primarily residuals of uncorrected atmospheric turbulence, but were surprised when then data collected using the internal light source exhibited the same behavior (in that case, they found τ = 6.3 s). They performed additional tests but were not able to identify the fundamental cause of this behavior, and in the end concluded that it was "real [and] related to an internal effect in the instrument independent of the atmospheric conditions or the telescope."

For images separated by greater than a minute, they observed a second regime of speckle change. They observed a linear decorrelation with time, and this was again present in both the on-sky data and also their internal lamp data. They suggested that this change in speckles must be caused by mechanical changes in the instrument (in particular the motion of the image derotator and thermal expansion effects). When they disabled the image rotator, this linear decorrelation of speckles became flat with time, i.e., the speckles were not measurably changing. In summary, they observed two regimes of speckle decorrelation, and these had different timescales for the on-sky and off-sky observations, but both were fundamentally caused by instrumental and not atmospheric effects.

Our data also consisted of H-band images of an ExAO-corrected PSF, but it had three orders of magnitude higher temporal resolution, enabling us to probe far shorter timescales. Additionally, because SCExAO is an entirely different instrument from SPHERE, we hoped our data might shed light on Milli et al. (2016)'s instrumental effects.

Our data had some important differences, however. First, the SPHERE data were collected behind a Lyot coronagraph, whereas our data were noncoronagraphic. Second, while their data were collected over 52 minutes, our data were 7 minutes in length, so we could not examine as long as timescales. We also recorded two 20-minute data sets of the SCExAO internal light source to differentiate between instrumental and atmospheric effects.

We then calculated Pearson's correlation coefficients according to Equation (5) using the annular region centered on the PSF and defined by radial separations 3λ/D <r < 10λ/D (0.12'' < r < 0.40''). Comparing each of the 630 thousand frames for the on-sky data and 1.99 million frames for the internal source data with each of the other frames in that data set would have been computationally prohibitive. Therefore, we compared each frame against 1, 2, 3, ..., 200, 208, 216, ..., 8400, 10080, 11760, ..., 604800 frames later. This enabled us to have maximum (S/N) at each calculated temporal separation, but we sacrificed temporal resolution at higher separations. In other words, we wanted the best time resolution at the shortest possible intervals, but settled for one Pearson's correlation coefficient data point every second at a time intervals greater than six seconds. This simplification reduced the computational load by 99.9%. Even after this, we calculated 852 million correlation coefficients for the on-sky data and 6.04 billion for the internal source data. The data points at a given temporal separation were then averaged to improve signal-to-noise. Because there are fewer frame pairs separated by large amounts of time, signal-to-noise decreased at larger separations. The on-sky data set was seven minutes long, and we calculated correlation coefficients for data points up to six minutes apart to still have reasonable S/N at the larger separations. The computation of the Pearson correlation coefficients for each plot took several CPU-days. Figures 6 and 7 show the results for the on-sky and internal light source data, respectively.

Zoom In Zoom Out Reset image size

We observed three timescales of speckle evolution in the on-sky data. First, and most interestingly, our high frame rate data revealed an effect that Milli et al. (2016) did not observe. Frames separated by ≲2 ms are highly correlated in the on-sky data. After this time, the exponential decay behavior described in the next paragraph dominates. The minimal change in data points separated ≲2 ms is not present in the internal source data, implying that this is an actual atmospheric effect. This timescale is similar to the SCExAO loop bandwidth (2 kHz update rate, 15% gain, and 975 μs latency gives a 3 dB bandwidth of 240 Hz), so it is most likely explained by the changes to the deformable mirror shape. Over timescales of >2 ms, atmospheric speckles become decorrelated because they are corrected by the ExAO loop. Temporal correlations that naturally exist at >2 ms in the atmospheric speckles are erased or greatly attenuated by the ExAO correction; so there is a reduction in speckles with timescales of >2 ms.

We recorded two 20-minute data sets using the SCExAO internal source. The first internal source data set had an unsaturated PSF core, mirroring our 2017 September 11 on-sky observations. However, because this was very nearly 100% Strehl (the slight degradation from 100% Strehl was due to imperfections in the optical elements of the instrument), very few speckles rose above the read noise of the detector. As a result, the calculation of the Pearson's correlation coefficients revealed no real speckle evolution (Figure 7 left plots). This mirrors the findings of Milli et al. (2016), who saw no speckle evolution when they froze the operation of SPHERE's optical elements. To improve signal-to-noise in the speckles, we then saturated the PSF core by a factor of 100 in the second internal source data set (Figure 7 right side). In this case, we observed two timescales of decorrelation as described in the following paragraphs. Although the PSF core of the on-sky data was unsaturated, we could measure speckle evolution there because the Strehl was lower (∼90% instead of ∼100%) and therefore the speckles were brighter.

Zoom In Zoom Out Reset image size

We observed a second timescale on the order of minutes over which the speckles linearly decorrelated in both the on-sky data and the saturated internal source data. We fitted a first-order polynomial to data points between 20 and 350 seconds and found that the speckles decorrelated by 56 parts per million (PPM) on-sky and 26 PPM on the saturated internal source over this period. For comparison, Milli et al. (2016) measured 25 PPM and 9 PPM for these two timescales, respectively, but their fit was performed over a different range of times. Because this decorrelation is present in both the on-sky and internal source observations, it is fundamentally caused by an instrumental effect. Unlike Milli et al. (2016), we did not use the image derotator; we had no moving components between the light source and our detector. This decorrelation must therefore be due to thermal expansions in the mechanics of SCExAO. This data depart from this linear decorrelation between 300 and 360 seconds in the on-sky data, but we think this is due to small number statistics (there are fewer frame pairs with 6 minute separations than, for example, 1 minute separations) than any intrinsic speckle behavior. The internal source data did not exhibit this effect.

The third timescale we observed was an exponential decay in the Pearson's correlation coefficients over the first few seconds for both the saturated internal source and the on-sky data. Milli et al. (2016) also observed this behavior, but in our case it occurred much more quickly. We fitted the points between 0.01 and 5 seconds with Equation (6) and found best-fitting parameters of τ = 0.40 s for the on-sky data and τ = 2.3 s for the saturated internal source. Our measured values for Λ and ρ₀ are provided on the plots, but these have fewer physical implications because they are influenced by the flux level on the detector, as demonstrated by Figure 7. In contrast, Milli et al. (2016) found exponential decay terms of τ = 3.5 s and τ = 6.3 s for their on-sky and internal light source data, respectively, when fitted to timescales between 1 and 30 seconds.¹¹ We observed a much more rapid decay than them, but the exponential decay behavior was present in both observations. Milli et al. (2016) were unable to explain the cause of this effect.

Whereas the analysis presented in Section 3.2 had applications to real-time speckle nulling, the results presented here are primarily applicable to speckle reduction during postprocessing. If the exposure times of individual images in an observing sequence are shorter than the timescales on which speckles evolve, then postprocessing techniques such as KLIP (Soummer et al. 2012) and LOCI (Lafrenière et al. 2007) will be better able to subtract the speckles. For this reason, it is advantageous to use short exposures while conducting observations for which the speckle halo is relevant (this must be balanced against the read noise of the detector, emphasizing the need for low-noise detectors such as SAPHIRA or electron-multiplying CCDs). Although it is preferable to have the speckles as dim as possible, as these have minimal photon noise, it is also ideal for the speckles that still exist to evolve as slowly as possible, because this enables them to be mitigated during data reduction. This emphasizes the need to have the instrument maximally mechanically and thermally stable. Indeed, SPHERE has more thorough environmental/thermal stabilization than SCExAO, which likely explains why the speckle decorrelations measured in this section occur more rapidly than those of Milli et al. (2016).