1.

Introduction

The advancement of adaptive optics (AO) as a technology has enabled significant progress in astrophysics. Notably, the Gemini Planet Imager (GPI) is one such instrument,¹ where fast and precise correction is pivotal to optimize instrument performance. The results so far have been spectacular. With detections of multiple planetary mass companions around various stars, as well as strong nondetection limits around many more, one can constrain planetary population distributions.² Astrometric measurements of multibody systems probe dynamical constraints on planetary masses and system lifetimes³ and provide a spectacular view of Kepler’s laws in action. Spectral measurements of individual giant planets constrain evolutionary and atmospheric models⁴ of these objects, paving the way toward characterization of extrasolar terrestrial planets.

For this generation of instruments and the next, understanding the point spread function (PSF) of AO instruments on giant telescopes will be important for the development of algorithms optimized in the search for planets.^5,6 The analysis in this paper expands on our previous work,⁷ which demonstrated the origin of azimuthal asymmetry in the PSF as a consequence of the time lag error, to explore asymmetry along the preferential axis introduced by scintillation. This effect has been demonstrated previously by Cantalloube et al.⁸ We will expand on their discussion using a more general method of analyzing the structure of the AO-corrected PSF analytically, as well as validating our conclusions with observations and atmospheric datasets. More specifically, our formalism demonstrates that the asymmetry grows linearly only for small spatial frequencies, and at higher spatial frequencies becomes nonlinear. We include solutions for the zeros of the log of the asymmetry metric, which are image locations with an observable return to symmetry.

The analysis in this paper is presented as a trident—theory, simulations, and observations. The first section derives the method of angular spectrum Fresnel propagation from the time-independent wave equation. This technique allows us to analytically calculate the PSF formed from a single mode of phase aberration that is both scintillated and time-lag corrected. The PSF for this single mode is computed to second order in a Taylor expansion, which is well matched when compared to a numerical solution involving the discrete Fourier transform. In the second section, these methods are extended broadly to waves propagating through an atmospheric model with Kolmogorov turbulence in the frozen flow regime, which reproduces the behavior in a moderately accurate simulation of an entire telescope employing AO. The third section demonstrates the effects from our analytic model are observable in real data taken from the Gemini Planet Imager. We then explore correlations in the observations when combined with a meteorological dataset containing the real wind velocities and directions in the atmosphere during the observations. Finally, this paper concludes with a brief discussion about the importance of this effect in the context of improving AO systems performance from design to postprocessing.

2.

Theoretical Scintillation Analysis

2.2.

Single-Mode Analysis

Now that we have described how free-space propagation affects the complex illumination of the electric field in the Fresnel regime, we will explore an analytic derivation of how this results in an asymmetric PSF for a single mode of phase aberration in one dimension.

When the inverse of the spatial frequency of the mode is much longer than the wavelength of the propagating light (or ${\lambda }^2{f}^2\ll 1$), a binomial expansion of the free-space propagation transfer function:

Eq. (7)

$$H(z)\approx e^{-i\pi \lambda f^{2}z}$$

allows us to perform the analysis up to a constant phase term. In the example of a single mode of sinusoidal phase aberration only (no amplitude errors), with period $p=1/f$, and a phase variation amplitude $\alpha$, the complex illumination can be written $U=A{e}^{i\varphi }$ with

Eq. (8)

$$A=1,$$

Eq. (9)

$$\varphi =\alpha \sin (2\pi x/p).$$

It can be shown¹⁰ that for small phase errors ($\alpha \ll 1$) the illumination vector at propagation distance $z$ takes the form $U(z)=A(z)e^{i\varphi (z)}$, where the phase and amplitude have appropriately been “mixed” due to the scintillation effects:

Eq. (10)

$$A(z)=1+\alpha \sin (2\pi z/z_T)\sin (2\pi x/p),$$

Eq. (11)

$$\varphi (z)=\alpha \cos (2\pi z/z_T)\sin (2\pi x/p).$$

Here $z_T$ is the Talbot length:

Eq. (12)

$$z_T=\frac{\lambda }{1-\sqrt{1-{\lambda }^2{f}^2}}\approx \frac{2}{\lambda f^2}.$$

With the assumption that an ideal AO system corrects phase aberrations after a short delay $t$ due to the servo-lag error, our AO-corrected phase will be the difference between two propagated modes, one shifted along the direction of the wind with the coordinate transform $x\to x-vt$, if the wind velocity $v$ points along the positive $x$ axis, and the other the initially measured phase:

Eq. (13)

$${\varphi }_{\mathrm{AO}}={\varphi }_z(x-vt)-{\varphi }_z(x),$$

Eq. (14)

$${\varphi }_{\mathrm{AO}}=\alpha \cos (2\pi z/z_T)\sqrt{2}\sqrt{1-\cos (2\pi vt/p)}\sin \{2\pi x/p+\arctan [\frac{-\sin (2\pi vt/p)}{\cos (2\pi vt/p)-1}\left]\right\}.$$

To simplify this expression, we will use the substitutions:

Eq. (15)

$$\mathcal{A}=\alpha \sin (2\pi z/z_T),$$

Eq. (16)

$$P=\alpha \cos (2\pi z/z_T)\sqrt{2}\sqrt{1-\cos (2\pi vt/p)},$$

Eq. (17)

$$\mathrm{\Delta }=-\frac{p}{2\pi }\arctan \left[\frac{-\sin (2\pi vt/p)}{\cos (2\pi vt/p)-1}\right],$$

where $\mathcal{A}$ represents the term that modulates the amplitude, $P$ is the term that modulates the phase, and $\mathrm{\Delta }$ is a term that represents a phase shift between the amplitude term and the phase term. With these substitutions, the expressions for our AO-corrected amplitude and phase become

Eq. (18)

$$A_{\mathrm{AO}}=1+\mathcal{A}\sin [2\pi (x-vt)/p],$$

Eq. (19)

$${\varphi }_{\mathrm{AO}}=P\sin [2\pi (x-\mathrm{\Delta })/p].$$

An example of a single mode of AO-corrected illumination and its corresponding electric field distribution in the image plane is demonstrated in Fig. 1. For this figure, the principal observation is how the delayed AO correction and scintillation produce residual errors which are spatially offset, and how the corresponding speckles in the image plane are asymmetric.

Fig. 1

An example of a single mode of AO-corrected complex illumination and corresponding asymmetric speckles in the PSF. Here $\alpha =.01\mathrm{rad}$, $z=10\mathrm{km}$, $p=1\mathrm{m}$, $v=10\mathrm{m}/\mathrm{s}$, and $t=3.2\mathrm{ms}$. The AO-corrected complex illumination shows small aberrations in both amplitude and phase from the Talbot mixing, with a spatial phase offset due to the delayed correction from the servo lag. The corresponding asymmetry is highlighted with red dots in the image plane, which is really the inverse Fourier transform, taken in the Fraunhofer diffraction limit. The large central speckle is due to the constant amplitude, on the left subtracted off from the displayed amplitude so that the small errors are more visible.

2.3.

Open-Loop Model Validation

To validate the assumption that an open-loop AO model can accurately reproduce the behavior of a true closed-loop AO system, the frequency responses of both open-loop and closed-loop AO systems were modeled. Building on the analysis done for GPI’s AO system (see Sec. 4.D of Ref. 1 for a detailed treatment), the standard AO component models and control parameters were used to generate error transfer functions (ETFs). The GPI AO ETF is shown in red in Fig. 2. The system has a read-compute delay of 1.4 ms and a maximum controller gain of 0.3. An equivalent open-loop model, black curve in this figure, was fit by adjusting the effective delay such the measurement is applied 3.2 ms after it was taken. As shown in this figure, these two models agree very well in terms of both magnitude and phase response in the range of 2 to 200 Hz. Given a maximum spatial frequency in the AO-corrected dark hole of $2.78{\mathrm{m}}^{-1}$, these valid temporal frequencies in our model correspond to wind velocities in the range 8 to $72\mathrm{m}/\mathrm{s}$. These velocities encompass the range of possible wind velocities we might see naturally occurring in the jet stream, which typically ranges from 10 to $60\mathrm{m}/\mathrm{s}$. The open-loop model disagrees at the lowest temporal frequencies, which corresponds to static errors in the system. This discrepancy can be addressed by slightly scaling the phase measurement, e.g., 

Eq. (20)

$${\varphi }_{\mathrm{AO}}(t)=\varphi (t)-0.997\varphi (t-\tau ),$$

where $\tau$ is the delay time. The open-loop model does not agree at higher temporal frequencies, most obviously when temporal frequency is the inverse of $\tau$. This is when an ideal open-loop AO system coincidentally achieves perfect correction of the translating Fourier mode, which is not physically realizable. This region is beyond our wind speeds of interest, so the open-loop approximation is suitable for our investigation.

Fig. 2

Bode plot for comparing the transfer function of various open- and closed-loop AO models. By adding an integrating factor to the open-loop model, it becomes possible to recreate the performance of a closed loop system on static errors, although these are outside of the temporal frequency domain we care about, which is approximately from 2 to 200 Hz. In this range, the standard open-loop model performs roughly as well as the integrator in mimicking the closed loop model, although the behavior diverges when the inverse of the temporal frequency becomes equal to the time delay.

2.5.

Taylor Expansion of the Single-Mode PSF

Although it may be possible to calculate the AO-corrected electric field distribution in the image plane by taking the inverse Fourier transform of the AO-corrected aperture illumination by expanding it as an infinite series of Bessel functions, it is not clear that this expression can be easily squared to get the intensity distribution due to an infinite amount of cross terms. Instead we will follow the conventions of Sivaramakrishnan et al.¹² and Perrin et al.¹³ to arrive at the intensity. However, both Sivaramakrishan et al. and Perrin et al. assume that the image plane PSF is the Fourier transform of the aperture plane and not the inverse. For all of their results, this fact does not matter but for ours we must be cautious, and remember to apply the appropriate coordinate transform to recover the proper result. For small $\varphi$, the intensity distribution or PSF $={|\mathcal{F}[u]|}^2$ can be expanded in a Taylor series:

Eq. (23)

$$\mathrm{PSF}\approx {\mathrm{PSF}}_0+{\mathrm{PSF}}_1+{\mathrm{PSF}}_{2,\text{halo}}+{\mathrm{PSF}}_{2,\text{strehl}},$$

whose first few terms are

Eq. (24)

$${\mathrm{PSF}}_0=a{a}^{*},$$

Eq. (25)

$${\mathrm{PSF}}_1=2\mathrm{Im}[a(a^{*}\star {\mathrm{\Phi }}^{*})],$$

Eq. (26)

$${\mathrm{PSF}}_{2,\text{halo}}=(a\star \mathrm{\Phi })(a^{*}\star {\mathrm{\Phi }}^{*}),$$

Eq. (27)

$${\mathrm{PSF}}_{2,\text{strehl}}=-\frac{1}{2}[a(a^{*}\star {\mathrm{\Phi }}^{*}\star {\mathrm{\Phi }}^{*})+a^{*}(a\star \mathrm{\Phi }\star \mathrm{\Phi })].$$

Here, the case change is used as shorthand for the Fourier transform, so $a=\mathcal{F}(A)$ and $\mathrm{\Phi }=\mathcal{F}(\varphi )$. Additionally, $\star$ is the convolution operation and $*$ is the complex conjugate. For our AO-corrected illumination, these can be calculated as follows:

Eq. (28)

$$a_{\mathrm{AO}}=\delta (f)+\frac{\mathcal{A}}{2i}\delta (f-1/p)e^{-i\frac{2\pi }{p}vt}-\frac{\mathcal{A}}{2i}\delta (f+1/p)e^{i\frac{2\pi }{p}vt},$$

Eq. (29)

$${\mathrm{\Phi }}_{\mathrm{AO}}=\frac{P}{2i}\delta (f-1/p)e^{-i\frac{2\pi }{p}\mathrm{\Delta }}-\frac{P}{2i}\delta (f+1/p)e^{i\frac{2\pi }{p}\mathrm{\Delta }},$$

where $\delta$ is the Dirac delta distribution. It is worthwhile to note that here we implicitly are using the entire real line in the Fourier transform, which can be interpreted as using an infinitely large telescope, or as a telescope with an idealized perfect coronagraph. This leads to a PSF with terms:

Eq. (30)

$${\mathrm{PSF}}_0=\delta (f)+\frac{{\mathcal{A}}^2}{4}[\delta (f-1/p)+\delta (f+1/p)],$$

Eq. (31)

$${\mathrm{PSF}}_1=\frac{-\mathcal{A}P}{2}\sin [\frac{2\pi }{p}(vt-\mathrm{\Delta })][\delta (f-1/p)-\delta (f+1/p)],$$

Eq. (32)

$${\mathrm{PSF}}_{2,\text{halo}}=\frac{P^2}{4}[\delta (f+1/p)+\delta (f-1/p)]+\frac{{\mathcal{A}}^2{P}^2}{16}\left[2\right\{1+\cos [\frac{4\pi }{p}(vt-\mathrm{\Delta })]\}\delta (f)+\delta (f+2/p)+\delta (f-2/p)],$$

Eq. (33)

$${\mathrm{PSF}}_{2,\text{strehl}}=\frac{-P^2}{2}\delta (f)-\frac{{\mathcal{A}}^2{P}^2}{32}\{4+2\cos [\frac{4\pi }{p}(vt-\mathrm{\Delta })\left]\right\}[\delta (f-1/p)+\delta (f+1/p)].$$

Indeed, the even order terms are symmetric and the first odd order term ${\mathrm{PSF}}_1$ is responsible for the observed asymmetry. We will examine the ratio of the amplitude for the speckles on either side of the image, and to do so evaluate the PSF at the location appropriate $f=\pm 1/p$:

Eq. (34)

$${\mathrm{PSF}}_0(f=-1/p)={\mathrm{PSF}}_0(f=1/p)=\frac{{\mathcal{A}}^2}{4},$$

Eq. (35)

$$-{\mathrm{PSF}}_1(f=-1/p)={\mathrm{PSF}}_1(f=1/p)=\frac{-\mathcal{A}P}{2}\sin [\frac{2\pi }{p}(vt-\mathrm{\Delta })],$$

Eq. (36)

$${\mathrm{PSF}}_2(f=-1/p)={\mathrm{PSF}}_2(f=1/p)=\frac{P^2}{4}-\frac{{\mathcal{A}}^2{P}^2}{32}\{4+2\cos [\frac{4\pi }{p}(vt-\mathrm{\Delta })\left]\right\}.$$

Our metric for the ratio of the right to left speckle asymmetry is the Taylor expansion sum of the PSF terms evaluated at these appropriate locations:

Eq. (37)

$$\chi =\frac{\mathrm{PSF}(f=1/p)}{\mathrm{PSF}(f=-1/p)},$$

which is a function of propagation distance $z$, mode period $p$, and velocity times delay $vt$. A comparison of the asymmetry metric for both the numerical single-mode scintillation using the discrete Fourier transform and our analytic second-order Taylor expansion is given in Fig. 3.

Fig. 3

The log of the speckle asymmetry ratio as a function of PSF location for parameters $\alpha =.01\mathrm{rad}$, $t=3.2\mathrm{ms}$, and $z=10\mathrm{km}$. The solid (or thin) lines represents the second-order Taylor expansion and the dashed (or thick) lines represent the numerical solution found using the discrete Fourier Transform. The $x$ axis is a proxy for mode length, transformed in PSF location when imaged in H-band at $\lambda =1.6\mu \mathrm{m}$. The zeros corresponding to the particular mode lengths $p$ for which the propagation distance $z$ is an integer multiple of $z_T/4$ are plotted in black stars (or squares). For these, ${\mathrm{PSF}}_1=0$ because either $\mathcal{A}\sim \sin (2\pi z/z_T)=0$ or $P\sim \cos (2\pi z/z_T)=0$, resulting in a symmetric PSF. Although there are additional zeros when the velocity times time delay is comparable to the mode length, visible in both the 30 and $50\mathrm{m}/\mathrm{s}$ case, though their analytic solution is more complicated [set $\frac{2\pi }{p}(\mathrm{\Delta }-vt)=\pi /2+n\pi$, $n\in \mathbb{Z}$].

Highlighting a few observations from this plot: first, we note that for wind layers at $z=10\mathrm{km}$, the first zero crossing is not until roughly 1 arcsec in the image, and so for most of the relevant observations, the strongest asymmetry will be on the side opposite to the direction of the wind, though for higher spatial frequencies corresponding to the edges of the PSF in the image plane, this will not always be true. Second, the strength of the asymmetry is greater for slower wind velocities in this region of interest. Third, it is worthwhile to note that our Taylor expansion generally gets the behavior of the asymmetry analytically correct, although differs from the numerical solution due to the absence of higher order terms in the expansion.

3.

Atmospheric Turbulence Scintillation Simulation

The response of an ideal open-loop AO system to frozen flow Kolmogorov turbulence was developed in simulations using the method of the angular spectrum free-space propagation and Fourier decomposition. These simulations can place many layers of turbulence at arbitrary altitudes, and we explore the effects of single and multiple layers. Readers interested in the precise details of these simulations are referred to Sec. 6, which discusses our techniques at greater length. However, here, we will discuss the results of these simulations in the context of the previous analysis.

To start, the effect of scintillation in the halo of an AO PSF is most clearly seen for the case of a single layer of Kolmogorov turbulence, which is plotted in Fig. 4 for three characteristic wind velocities. The single-layer simulation uses a turbulent layer altitude of $z=10\mathrm{km}$, which roughly corresponds to the altitude of the jet stream. The layer altitude equals the propagation distance of the light when the telescope is pointed at zenith; however, for low elevation pointing, the propagation distance and therefore strength of the scintillation, will grow larger. The jet stream is the strongest layer of turbulence which has enough relevant altitude to scintillate significantly [recall the strength of the scintillation amplitude errors grow like $\sin (2\pi z/z_T)$].

Fig. 4

For a single layer of atmospheric turbulence at $z=10\mathrm{km}$, the scintillation halo at low wind velocities is highly apparent and still noticeable at rapid velocities. Each image has its own unique colorbar, so that the variation in the halo shape is visible, although from looking at the magnitude of the halo intensity, it is clear that slower wind velocities are corrected better in the metric of total scattered light.

Fig. 5

For simulations with many wind layers, here $L=18$ is the number of layers used, the effect of the scintillation asymmetry is significantly less apparent, although at low wind velocities where the asymmetry is strongest, it is still visible. This is likely due to interference from the other wind layers causing the delicate spatial offset in the illumination to become washed out when combined with the finite sampling resolution of our simulation, and for real observations, the detector.

The wind in the single-layer example is propagating directly to the right, along the positive $x$ axis. For the multilayer case, the wind velocity for the jet stream is generally pointing toward 2 o’clock (Fig. 5). As a consequence, we can see that the brighter asymmetric lobe of the PSF is on the opposite side of wind direction. To remain consistent with observations, it is necessary that the image plane is the inverse Fourier transform of the aperture plane in the Fraunhofer limit. Looking at just the single layer, it is apparent that the scintillation dominates for slower wind velocities. This scintillation halo would be close to the true PSF if the atmosphere was actually only a single layer of turbulence. Examining the edges of the PSF, one can identify the region where the asymmetry metric $\log \chi =0$ that was described analytically, noticeable here as dark bands.

The halos for the single-layer case do not match real observations well. However, when we run this simulation with many wind layers, each with unique velocities pulled from an instance of the NOAA GFS to mimic real observing conditions, the extreme scintillation halo washes out. This is likely an interference between multiple independent Kolmogorov layers masking the delicate spatial offset for a single mode needed to create the asymmetry demonstrated earlier. Although at a quick glance the PSF does not appear to vary with the scaled wind velocities cited in their titles, upon closer inspection the asymmetry remains, and is most obvious for the slowest wind velocity. For most observations, with jet stream velocities upward of $50\mathrm{m}/\mathrm{s}$, there are no apparent deviations from the symmetry of the butterfly shaped halo that is often seen. But in the rare cases when the high-altitude winds are very slow, around $10\mathrm{m}/\mathrm{s}$, the scintillation becomes stronger and introduces a noticeable asymmetry in the PSF along the axis of the wind propagation.

4.

Observational Correlations

The previously discussed asymmetry for AO PSFs is observable in real data taken from observations using the Gemini Planet Imager. Over the course of a few years of observations, over 20,000 images were taken in H-band as a part of the GPI Exoplanet Survey. We exercise selection cuts to find a sample of PSFs exhibiting this characteristic asymmetry using the following methods.

The first metric is the fractional standard deviation (FSD) in an annulus 30 to 70 pixels wide centered around the star’s location. We first select images with FSD greater than one standard deviation above the mean. This selects for images where the relatively fast jet stream turbulence was dominant, producing PSFs with the characteristic butterfly shaped halo. Then we construct an angular profile of the image by integrating along the radial axis. This allows us to fit for the preferential wind direction modulo 180 deg. These techniques follow exactly from our previous work,⁷ yet in this work, we expand by defining the asymmetry metric $\chi$ between the two lobes of the PSF. $\chi$ is defined as the ratio of the summed intensities in the half-annulus perpendicular to the axis of wind propagation. This asymmetry metric is used to break the 180-deg symmetry of the wind axis to find the true asymmetric butterfly vector in the image plane. In addition, our metric $\chi$ is folded around 1 by inversion such that the half-annulus with greater intensity is always in the numerator, meaning $\chi$ always takes on values greater than 1, and larger values correspond to greater asymmetry. Note that this $\chi$ is not identical to the $\chi$ described earlier for the one-dimensional case, although both are intensity ratios, this one lives in two dimensions. A large sample of PSFs selected in this manner as displayed in Fig. 6.

Fig. 6

A selection of images taken with the Gemini Planet Imager, ordered by decreasing asymmetry metric $\chi$. Each image is the most strongly asymmetric image taken from each set of observations of a single target, such that there are no duplicates, which shows that this effect is often recurring and not limited to single cases. In addition, each image has been rotated and flipped such that north is up and east is to the left. The angular size of each GPI image is $2.8”\times 2.8”\mathrm{arc}\sec$. In each image, the direction of the strong asymmetry is plotted with a red arrow, and the corresponding direction in degrees azimuth projected onto the ground is printed in the corner of each image, alongside the ratio of the asymmetry metric $\chi$. Glancing over the entire dataset, it is readily apparent that the strong asymmetry direction is predominately pointing into the west.

These asymmetric butterfly vectors are projected from the image plane onto the surface of the Earth, for comparison with meteorological wind data. Readers interested in the trigonometric problem of relating the two are referenced to our techniques in Sec. 7. The directions and velocities of the wind for various layers in the Earth’s atmosphere are taken from the NOAA Global Forecast System,¹⁴ archives of which are available to the public. The distribution of wind directions for the jet stream and the distribution of directions for the asymmetric butterfly vector projected on the ground are displayed in Fig. 7 for a concise comparison. It is apparent that the jet stream predominately points East in accordance with its origin due to the rotation of the Earth. As a consequence, the resulting asymmetric butterfly vector points entirely toward the West. This observation confirms our prior analysis regarding the strength of the asymmetry being opposed to the direction of the wind.

Fig. 7

Angular histogram of the asymmetric butterfly vector’s projection onto the ground for our sample in degrees from azimuth alongside the direction of the jet stream. The asymmetric butterfly nearly always points west while the jet stream is nearly always points east, which is only possible if the image plane electric field distribution is the inverse Fourier transform of the aperture plane. These distributions are not temporally matched but are rather the entire subset observations with strong asymmetry and the entire distribution of jet stream wind directions over the course of the survey. For temporally matched correlations, see Figs. 8 and 11.

To further verify this correlation between the jet stream and the asymmetric butterfly, we take our sample of data where $\chi \ge 1.2$, which is where the asymmetry begins to be noticeable to the eye, and match the time of observation to time of the wind data for the approximate location of Cerro Pachon, Chile, where the Gemini South Telescope takes it observations from, and plot the resulting correlations for the various directions of the wind layers in Fig. 11. For most of the wind layers in the NOAA GFS, there is little to no correlation between the direction of the asymmetric butterfly and the direction of the wind, with the exception of the wind layers around 100 to 300 hPa, which roughly corresponds to altitudes of 10 to 15 km, which is the approximate altitude of the jet stream. The strength of these correlations, as measured by the Pearson $R$ coefficient, and the slopes of the best fit lines are plotted for a concise summary in Fig. 8. The rise of the correlations well beyond the limits imposed by a null hypothesis bootstrap show that these strong correlations are not spurious and a very real phenomenon.

Fig. 8

Distilling the information from the scatter plots in Fig. 11 into their slopes and Pearson $R$ coefficient as a function of altitude from the GFS. The model contains layers all the way to sea level despite the observatory being around 3 km up because of the uniform grid spacing. The strong correlations are visible for the altitudes relevant to the jet stream in both the slope of the best fit line and the $R$ coefficient. These correlations are much greater than the 2 to 4 sigma chances concerning the null hypothesis, generated from bootstrap sampling at random, shown as the shaded regions and solid lines for the slope and the $R$ coefficient, respectively.

To further explore our observations of the asymmetry, we plot $\chi$ versus the velocity of the jet stream and the elevation of the telescope pointing for our subset of observations which have been temporally matched to the wind data in Fig. 9. Both plots have reasonable correlation strength with $R$ coefficient around negative one half, and both of these correlations have an intuitive sense. Our analysis showed that the asymmetry should be strong for slower wind velocities, and this is indeed verified through the first plot. The second has contributions from two different effects. When the telescope is pointing toward zenith ($\text{elevation}=90\mathrm{deg}$), the turbulent wind layers in the jet stream are as close to the telescope as possible along the line of sight. As a consequence, the propagation distance $z$ will be smaller relative to the Talbot length $z_T$ in comparison to when the telescope is pointed lower toward the horizon. A lower elevation value corresponds to the telescope effectively pushing the layer for a given altitude to a further propagation distance, giving the light more time to scintillate. This pushes the relevant mode lengths for the asymmetry closer to the center of the image, see Fig. 12 for a remake of Fig. 3 with $z=25\mathrm{km}$. In addition, a low telescope elevation also has the effect of changing the apparent wind velocity over the aperture. A telescope pointing away from zenith can only observe relatively slower wind velocities than one pointing directly up as long as the winds are parallel to the surface of the Earth. These two effects work together to build the strong correlation between telescope elevation and asymmetry.

Fig. 9

Asymmetry strength $\chi$ when compared to the velocity of the jet stream and the telescope elevation for our sample. Both exhibit strong correlations which corroborate our analytic understanding of the origin of the asymmetry. Slower jet stream velocities directly cause stronger asymmetry, while decreased telescope elevation has two effects. One to decrease the apparent wind velocity and the other to push the turbulent layers father away, giving more distance to scintillate.

5.

Discussion

High-contrast imaging systems, with large actuator counts, are often limited by time lag errors. In particular, when observing bright stars, they are the dominant source of scattered light within the “dark hole” region.¹⁵ Under mid-latitude Chilean conditions, the velocity of the jet stream is often the dominant source of these errors, even if its contribution to the total $r_0$ is moderate. This implies that different observatory sites with slower winds may have a comparative advantage, as well as the merits in scheduling observations around poor atmospheric conditions. In our previous work,⁷ we demonstrate that the jet stream is highly correlated with these errors, in a very large sample of observations.

Scintillation has previously been identified as a performance limiting factor in high-contrast imaging,¹⁶ and here we have demonstrated the severity of this effect, visible in the form of PSF asymmetry. This has implications for both current and future AO systems, especially but not exclusively those designed for high-contrast imaging. Many analyses often assume that various sources of scattered light are uncorrelated¹⁶ for simplicity. In this paper, we show that correlations between scintillation amplitude errors and time lag phase errors exist and can dominate during ideal conditions. If one simply tries to minimize the time lag error alone, without additionally compensating amplitude errors, the asymmetry will become larger as the effective wind velocity is decreased. This effect is worst when considering the asymmetry in low-order modes, which correspond to small separations in the final images, the region where planets or protoplanetary disks are most likely to be found. Although it is still optimal to have as-fast-as-possible correction in the metric of total scattered light in the halo, as computers get faster and algorithms are optimized the asymmetry will begin to play a larger role relative to other errors in the final PSF formed in AO images, leading many to explore possible routes for correction.

Scintillation errors are uncorrectable for an AO system operating in single-DM phase-conjugation mode, even if it is infinitely fast (see Fig. 13). This in turn will set a performance floor for such systems, and the asymmetry detected here provides a first measurement of the level in which those effects begin to dominate. One could address this with a system that corrects amplitude errors as well, using the Talbot or scintillation mixing effect to one’s own advantage. Having two deformable mirrors at two unique conjugate planes in the optical system enables some phase introduced by one DM to transform into amplitude, allowing one to correct amplitude errors from scintillation. This concept has been proposed for space-based coronagraphs¹⁷ to correct static amplitude errors, and laboratory testing is underway.¹⁸ Similar designs have been expressed for improvements in laser communications.^19,20 But we are particularly interested in the future of high-contrast imaging, particularly in the era of ELTs. Such a system could also correct static amplitude errors in an ELT, such as reflectivity variations between segments which have been recoated at different times.^21,22 Driving such a system would require knowledge of both the phase and amplitude of the science wavefront. Space-based coronagraphic instruments can achieve this using the science camera and making several measurements while modulating the speckle field with the DM.^23,24 Another similar focal plane wavefront sensing approach was recently proposed on ground-based telescopes,²⁵ which could also correct for amplitude aberrations from scintillation. Conventional Shack–Hartmann wavefront sensors measure some intensity information but this is complicated by spots in each subaperture moving out the active pixels of the detector, so it is impractical to separate phase and amplitude. Other methods to estimate the amplitude errors could utilize fast interferometric focal-plane sensing,²⁶ or something as simple as adding a high-speed direct pupil-imaging channel to a traditional AO wavefront sensor.

The scintillation halo observed in AO PSFs is a challenge for effective postprocessing of datasets. Many algorithms are designed to subtract a static speckle field with respect to the detector plane whose origin is from optical phase errors from imperfections, misalignment, noncommon-path errors (NCPE), and other systematic sources. Since these speckles have a unique spectral dependence, moving to farther separations at longer wavelengths, instruments imaging with an integral field spectrograph can measure that spectral dependence and remove those aberrations. However, since the scintillation halo intensity is driven by the effectiveness of the AO system, the halo appears fainter at longer wavelengths. Because the halo originates from a spectrum of turbulent modes which decide their final image plane location, the halo does not scale in image location with wavelength, unlike the static speckles best removed with SDI.²⁷ In addition, the scintillation halo is fixed with respect to the direction of the high-altitude winds, and it does not track with the rotation of the instrumental errors or the parallactic rotation of the astrophysical signal in an observing strategy such as ADI.²⁸

Various different postprocessing algorithms^29,30 often use a high-pass filter (HPF) to attempt to eliminate the diffuse background halo, and this is effective for regions of the image at large separations. However, near the coronagraphic mask, the scintillation can have very sharp features, demonstrated analytically as regions where $\log \chi =0$, and observable in simulations as dark regions perpendicular to the axis of the wind direction. When an HPF that preserves the features of a planet is applied to this halo, residuals which vary on spatial scales comparable to the planet are not removed. Often, a quadrupolar residual artifact near the coronagraphic mask if left which is large compared to the speckle residuals in the smooth halo at large separations. These residuals contribute significant noise to planetary detection attempts at the nearest separations, where the likelihood of detection is highest from their population distributions.² In addition, imaging extended objects such as debris disks cannot utilize an HPF in postprocessing, implying a limit to sensitivity even at wide separations for diffuse unpolarized structure.

Various methods of subtracting the scintillation halo have been suggested. One could take an empirical approach, using a PCA style analysis to model the shape of the halo over an averaged population of observations. Acknowledging this effect is unique and must be treated independently with this sort of approach can be effective at improving the final SNR in your detection algorithm.² Another approach may be to model the PSF end-to-end with complete simulations of the instrument and atmosphere, although this approach is significantly limited by the extent to which your simulated instrument can account for all real sources of error. Not only do errors arise from instrumental effects like DM fitting and NCPE, but also the non-Kolmogorov deviations in the turbulent spectrum from environmental effects,³¹ as well as the finite temporal resolution of available atmospheric information. Another path may attempt to estimate the PSF using a reconstruction from measured AO telemetry. However, with current WFS measurements, estimating the amplitude error from scintillation is rather difficult, as current instruments are not designed to measure wavefront amplitude. It is likely that the optimal approach to handling these errors is at the instrument level itself, as discussed previously, with a method to measure and correct the wavefront amplitude in real time. As high-contrast imaging strives for higher and higher performance levels, identification, measurement, estimation, and mitigation of scintillation errors will become increasingly important.

6.

Appendix A: Propagation through the Atmosphere

Tartarski³² has shown that the fluctuations in the optical index of refraction in three dimensions for a Kolmogorov turbulence spectrum follow the form:

The optical path length of a wavefront traversing a turbulent layer in the atmosphere from zenith can be found to first order by integrating the index of refraction over the thickness of the layer, and the accumulated phase is simply the wave vector of the ray $k=2\pi /\lambda$ times the optical path length:

For nonzenith observations, an additional term of $\sec \zeta$ where $\zeta$ is the zenith angle should be included in the integral in Eq. (41) to account for additional atmospheric depth. When the accumulated phase on the aperture is very large, we can subtract off the average phase, which is equivalent to removing the piston term from a Zernike polynomial.³⁵

In order to account for scintillation, each turbulent layer must be propagated according to the angular spectrum rule derived at the beginning of this paper. In order to make this simulation numerically tractable, the propagation through the turbulent phase screen is discretized into two steps, one where the phase is first accumulated according to the entire thickness of the layer at the start, and then second where the wave free-space propagates the entire distance of the layer. For an infinite number of layers, this assumption should recover the true propagation and indeed we are in a regime where the layer thickness is relatively small compared to the total propagation distance. A shorthand summary rule for the angular spectrum propagation is that the complex illumination at propagation distance $z$ is related to the complex illumination at the origin with

Here $\alpha$ and $\beta$ are the coordinates in the image plane, $U_{\mathrm{AO}}$ is our AO corrected complex illumination which includes amplitude errors from scintillation, $\mathcal{A}$ is the aperture function with the Blackman window apodization,³⁶ parameterized radially from the center with $r^2=x^2+y^2$:

7.

Appendix B: Relationship between the Sky and the Ground

The back of the telescope (BT) plane is the simplest way to imagine the relationship between an image on the sky and its orientation relative to the ground. Suppose you have a DSLR on a tripod, or a multimillion-dollar telescope with an Alt-Az tracking system. Either way, your imaging device is pointed at the celestial sphere along the line of sight vector:

Fig. 10

Diagram of coordinates used to orient images on the sky demonstrating the relationship between the ground plane and the BT plane.

Fig. 11

Correlations between the directions of the strong asymmetry of the image PSF and the wind direction for various wind layers in the NOAA GFS. Most wind layers do not exhibit significant correlation, with the exception of the layers around 100 to 250 hPa, which are the pressures corresponding the jet stream, at around 10 to 15 km of altitude.

Fig. 12

The log of the speckle asymmetry ratio for a single mode with a propagation distance of $z=25\mathrm{km}$. Such layers in the atmosphere viewed at zenith are very sparse due to the exponential decline in pressure of the Earth’s atmosphere, and so do not contribute significantly into observations done at zenith. For a jet stream layer at 15 km, a propagation distance of 25 km corresponds to telescope pointing elevation of 40 deg. When comparing this figure to Fig. 3, it becomes clear that modifying the propagation distance $z$ effectively changes the particular mode lengths the asymmetry occurs at to larger mode lengths, pushing the asymmetry to smaller separations, and nearer to the coronagraphic mask.

Fig. 13

Horizontal slice through the center of the PSFs generated in the single-layer atmospheric simulation, overlaid on top of each other for direct comparison. Here it is more explicitly visible as the wind speed decreases how both the total intensity of the halo decreases, as expected for better corrections, but in addition, the asymmetry becomes stronger at lower wind velocities, as expected from our model. For an infinitely fast AO system, the effective wind velocity is $v=0\mathrm{m}/\mathrm{s}$, which demonstrates the residual uncorrected amplitude error.

With such conventions laid out, it becomes easy to identify the location of the image plane on the sky, as it must be perpendicular to the line of sight. Since there are infinitely many such planes, we will use the convention

So that one can think of $\widehat{a}$ as pointing in the direction of increasing azimuth and $\widehat{b}$ pointing toward increasing elevation. It is left to the reader to show that $\widehat{a}·\widehat{b}=0$, and that $\widehat{a}\times \widehat{r}=\widehat{b}$ to verify the orthogonality of these unit vectors as a coordinate system.

With this elaborate set up, it becomes easy to convert vectors in the image plane into vectors in full three-dimensional space, and then project them onto the ground plane. Suppose we have a wind vector which appears in the image plane rotated $\psi$ from $\widehat{a}$ counterclockwise. Such a wind vector is

However, we would instead like to know $\widehat{w}(\widehat{x},\widehat{y},\widehat{z})$. By algebraically substituting in our coordinate vectors $\widehat{a}$ and $\widehat{b}$ formulas in $x$, $y$, $z$ space, we can arrive at an expression for the wind vector in $x$, $y$, $z$ space in terms of $\psi$, az, and el. This is

With this done, we can easily project the vector onto the ground plane by simply removing the $z$ component. If we need to find the direction of this wind vector as an azimuth, we can use the following trick:

However, images in the GPIES are not simply oriented as in the BT plane, but rather can be arbitrarily arranged due to the complexities of postprocessing. Fortunately for us, the orientation of each of the image has been previously calculated in celestial coordinates. These are represented as a CD matrix, which describe how $x$ and $y$ in pixels for the image correspond to right ascension and declination. Using the local sidereal time of the image during the exposure, it is possible to convert coordinates in right ascension and declination to coordinates in azimuth and elevation, using