An explanation for the origin of asymmetry along the preferential axis of the point spread function (PSF) of an AO system is developed. When phase errors from high-altitude turbulence scintillate due to Fresnel propagation, wavefront amplitude errors may be spatially offset from residual phase errors. These correlated errors appear as asymmetry in the image plane under the Fraunhofer condition. In an analytic model with an open-loop AO system, the strength of the asymmetry is calculated for a single mode of phase aberration, which generalizes to two dimensions under a Fourier decomposition of the complex illumination. Other parameters included are the spatial offset of the AO correction, which is the wind velocity in the frozen flow regime multiplied by the effective AO time delay and propagation distance or altitude of the turbulent layer. In this model, the asymmetry is strongest when the wind is slow and nearest to the coronagraphic mask when the turbulent layer is far away, such as when the telescope is pointing low toward the horizon. A great emphasis is made about the fact that the brighter asymmetric lobe of the PSF points in the opposite direction as the wind, which is consistent analytically with the clarification that the image plane electric field distribution is actually the inverse Fourier transform of the aperture plane. Validation of this understanding is made with observations taken from the Gemini Planet Imager, as well as being reproducible in end-to-end AO simulations.

## 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,^{1} 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.^{2} Astrometric measurements of multibody systems probe dynamical constraints on planetary masses and system lifetimes^{3} and provide a spectacular view of Kepler’s laws in action. Spectral measurements of individual giant planets constrain evolutionary and atmospheric models^{4} 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,^{7} 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.^{8} 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.1.

### Angular Spectrum Fresnel Propagation

To evaluate the effects of Fresnel propagation and scintillation, we derive the method of the angular spectrum. For an arbitrary complex illumination $U$ of the electric field, we can consider its decomposition into its angular spectrum $\tilde{U}$ given by the two-dimensional inverse Fourier transform in the plane $(x,y)$ at some constant $z$:

## Eq. (1)

$$U(x,y,z)=\underset{-\infty}{\overset{+\infty}{\iint}}\tilde{U}({k}_{x},{k}_{y},z){e}^{i({k}_{x}x+{k}_{y}y)}\mathrm{d}{k}_{x}\text{\hspace{0.17em}}\mathrm{d}{k}_{y}.$$Direct application of the Helmholtz equation:

to this decomposition will allow us to derive a formula for the propagation through free space of an arbitrary illumination. For a wave propagating with wave vector $\overrightarrow{k}={k}_{x}\widehat{x}+{k}_{y}\widehat{y}+{k}_{z}\widehat{z}$, implying ${k}^{2}={|\overrightarrow{k}|}^{2}={k}_{x}^{2}+{k}_{y}^{2}+{k}_{z}^{2}$, we find that propagation along the*z*axis is constrained by the second-order ordinary differential equation: where we have implicitly included the dependence of $\tilde{U}$ on the particular mode $({k}_{x},{k}_{y})$. This differential equation permits solutions of the form: where the $\pm $ represents a wave traveling in the positive or negative $z$ direction. For a mode with ${k}_{x}$, ${k}_{y}$, and the magnitude of the wave vector constrained by $|k|=2\pi /\lambda $, where $\lambda $ is the wavelength of the propagating light (assumed monochromatic), this means we can find the angular spectrum at some later plane $z$ from the angular spectrum at the origin $z=0$ by simply multiplying by the free-space propagation transfer function $H(z)$, which takes on the form:

## Eq. (5)

$$H({k}_{x},{k}_{y},z)=\mathrm{exp}[\pm i\sqrt{{\left(\frac{2\pi}{\lambda}\right)}^{2}-({k}_{x}^{2}+{k}_{y}^{2})}z].$$Making the substitution ${k}_{x,y}=2\pi {f}_{x,y}$ to represent the true linear wavenumber as in Goodman^{9} recovers this expression for the free-space propagation transfer function:

## Eq. (6)

$$H({f}_{x},{f}_{y},z)=\mathrm{exp}[\pm i\frac{2\pi}{\lambda}\sqrt{1-{\lambda}^{2}({f}_{x}^{2}+{f}_{y}^{2})}z].$$This allows us to consider the propagation of light through free space as a linear spatial filter applied to each Fourier mode of the complex illumination independently.

## 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:

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}$ withIt can be shown^{10} 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. (11)

$$\varphi (z)=\alpha \text{\hspace{0.17em}}\mathrm{cos}(2\pi z/{z}_{T})\mathrm{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. (14)

$${\varphi}_{\mathrm{AO}}=\alpha \text{\hspace{0.17em}}\mathrm{cos}(2\pi z/{z}_{T})\sqrt{2}\sqrt{1-\mathrm{cos}(2\pi vt/p)}\mathrm{sin}\{2\pi x/p+\mathrm{arctan}[\frac{-\mathrm{sin}(2\pi vt/p)}{\mathrm{cos}(2\pi vt/p)-1}\left]\right\}.$$To simplify this expression, we will use the substitutions:

## Eq. (16)

$$P=\alpha \text{\hspace{0.17em}}\mathrm{cos}(2\pi z/{z}_{T})\sqrt{2}\sqrt{1-\mathrm{cos}(2\pi vt/p)},$$## Eq. (17)

$$\mathrm{\Delta}=-\frac{p}{2\pi}\text{\hspace{0.17em}}\mathrm{arctan}\left[\frac{-\mathrm{sin}(2\pi vt/p)}{\mathrm{cos}(2\pi vt/p)-1}\right],$$## Eq. (19)

$${\varphi}_{\mathrm{AO}}=P\text{\hspace{0.17em}}\mathrm{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.

## 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\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{-1}$, these valid temporal frequencies in our model correspond to wind velocities in the range 8 to $72\text{\hspace{0.17em}\hspace{0.17em}}\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\text{\hspace{0.17em}\hspace{0.17em}}\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.,

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.## 2.4.

### Fraunhofer Diffraction Limit

According to Hecht,^{11} the relationship between the aperture and image planes taken in the far-field limit or the Fraunhofer diffraction limit is simply the Fourier transform, and their derivation results in the expression:

## Eq. (21)

$$E({k}_{x},{k}_{y})=\underset{-\infty}{\overset{+\infty}{\iint}}\mathcal{A}(x,y)\mathrm{exp}[i({k}_{x}x+{k}_{y}y)]\mathrm{d}x\text{\hspace{0.17em}}\mathrm{d}y.$$However, Goodman^{9} also claimed that the relationship between the two is the Fourier transform, yet they obtain the expression:

## Eq. (22)

$$U(x,y)=\frac{{e}^{ikz}{e}^{ik/2z({x}^{2}+{y}^{2})}}{i\lambda z}{\iint}_{-\infty}^{\infty}U(\xi ,\eta )\mathrm{exp}[-i\frac{2\pi}{\lambda z}(x\xi +y\eta )]\mathrm{d}\xi \text{\hspace{0.17em}}\mathrm{d}\eta .$$Comparing the two different expressions, it is not immediately obvious that they are nearly identical. However, after identifying the electric field $E$ and $U$, identifying the geometric relationship between the coordinates in the aperture $(x,y)\to (\xi ,\eta )$ and in the image plane $({k}_{x},{k}_{y})\to (\frac{2\pi x}{\lambda z},\frac{2\pi y}{\lambda z})$, and ignoring the phase prefactors (which do not affect the final intensity), the two answers are comparable with the exception of the sign in the exponent.

Since the choice of defining which Fourier transform is the forward and which is the inverse is arbitrary, both authors choose opposite sign conventions to arrive at the same conclusion that the relationship is the forward transform. However, the physical relationship between the two planes should not have this sign ambiguity. This difference is traceable to an assumption at the beginning of the derivations, where the choice of the direction of phasor rotation in a spherically converging or diverging wave ${e}^{\pm i(kr-\omega t)}/r$ is used as a Green’s function to solve the Huygens–Fresnel diffraction integral.

In practice, the difference between the forward and inverse transforms is essentially a coordinate transform from $x\to -x$ and $y\to -y$, and so the orientation of the PSF is flipped along both axes. For most cases, where the PSF is symmetric, this does not matter. However, for our purposes, properly orienting the PSF is critical and so taking note of this fact is important. Later, we will show that in order to remain consistent with observations, the Fraunhofer diffraction limit should use the inverse transform, with a positive sign in the exponent, which results in the stronger asymmetric lobe of the PSF on the opposite side as the wind.

## 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.^{12} and Perrin et al.^{13} 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}},$$## Eq. (25)

$${\mathrm{PSF}}_{1}=2\text{\hspace{0.17em}}\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}},$$## 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}\text{\hspace{0.17em}}\mathrm{sin}[\frac{2\pi}{p}(vt-\mathrm{\Delta})]\text{\hspace{0.17em}}[\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+\mathrm{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\text{\hspace{0.17em}}\mathrm{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. (35)

$$-{\mathrm{PSF}}_{1}(f=-1/p)={\mathrm{PSF}}_{1}(f=1/p)=\frac{-\mathcal{A}P}{2}\text{\hspace{0.17em}}\mathrm{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\text{\hspace{0.17em}}\mathrm{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:

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.Highlighting a few observations from this plot: first, we note that for wind layers at $z=10\text{\hspace{0.17em}\hspace{0.17em}}\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\text{\hspace{0.17em}\hspace{0.17em}}\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 $\mathrm{sin}(2\pi z/{z}_{T})$].

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 $\mathrm{log}\text{\hspace{0.17em}}\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\text{\hspace{0.17em}\hspace{0.17em}}\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\text{\hspace{0.17em}\hspace{0.17em}}\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,^{7} 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.

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,^{14} 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.

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.

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\text{\hspace{0.17em}\hspace{0.17em}}\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\text{\hspace{0.17em}\hspace{0.17em}}\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.

## 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.^{15} 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,^{7} 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,^{16} 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^{16} 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^{17} to correct static amplitude errors, and laboratory testing is underway.^{18} 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,^{25} 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,^{26} 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.^{27} 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.^{28}

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 $\mathrm{log}\text{\hspace{0.17em}}\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.^{2} 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.^{2} 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,^{31} 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^{32} has shown that the fluctuations in the optical index of refraction in three dimensions for a Kolmogorov turbulence spectrum follow the form:

^{33}with

## Eq. (39)

$$\delta N(\overrightarrow{x},z)=\mathrm{Re}\left\{{\mathcal{F}}^{-1}\right[\xi (\overrightarrow{\kappa},z)\sqrt{{\mathrm{\Phi}}_{N}(\kappa ,z)}\left]\right\},$$## Eq. (40)

$${\eta}_{ab}={\mathcal{F}}^{-1}({\tilde{\eta}}_{pq})=\sum _{p=0}^{P-1}\sum _{q=0}^{Q-1}{\tilde{\eta}}_{pq}\text{\hspace{0.17em}}\mathrm{exp}\left[2\pi i\right(\frac{pa}{P}+\frac{qb}{Q}\left)\right],$$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:

## Eq. (41)

$${\varphi}_{i}(\overrightarrow{x})=k{\int}_{{z}_{i}}^{{z}_{i}+\mathrm{\Delta}{z}_{i}}n(\overrightarrow{x},z)\mathrm{d}z=kn(\overrightarrow{x},z)\mathrm{\Delta}{z}_{i}.$$^{34}with where $\rho $ and $T$ are the pressure (in millibars or equivalently hPa) and temperature (in Kelvin) of the atmosphere for a particular altitude. To obtain a model for the index of refraction as a function of altitude, one can model the pressure as a decaying exponential with a scale height given by the local surface gravity $g$, the mean molecular mass $M$ of the atmosphere, the ideal gas constant $R$, and an assumed isothermal uniform temperature of the surface $T$ with where ${\rho}_{0}$ is the atmospheric pressure at sea level. An atmospheric temperature profile as a function from altitude can be determined empirically, or the values given in the GFS can be used, but small fluctuations in $T$ hardly affect the end value of the index of refraction, compared to the pressure, which dominates.

For nonzenith observations, an additional term of $\mathrm{sec}\text{\hspace{0.17em}}\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.^{35}

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

where $H$ is the free-space propagation transfer function which is implicitly also a function of the particular modes $k$ being propagated. With the complex illumination given at the aperture by the previous description, we invoke the Taylor frozen-flow hypothesis, which requires that the timescale for turbulence is much greater than the time delay with which the AO system will respond. For our simulation, this means that the fluctuations in the field of view simply propagate by translations due to the wind velocity, which can be expressed by## Eq. (45)

$$\delta N[\overrightarrow{x}+\overrightarrow{v}(z)\tau ,{t}_{0}+\tau ]=\delta N(\overrightarrow{x},{t}_{0}),$$## Eq. (46)

$${\varphi}_{\mathrm{AO}}=\varphi (\overrightarrow{x},{t}_{0})-\varphi (\overrightarrow{x},{t}_{0}-\tau ),$$^{11}

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,^{36} parameterized radially from the center with ${r}^{2}={x}^{2}+{y}^{2}$:

## Eq. (48)

$$\mathcal{A}(r)=\frac{1-\gamma}{2}-\frac{1}{2}\text{\hspace{0.17em}}\mathrm{cos}\left[\frac{2\pi (r-D/2)}{D}\right]+\frac{\gamma}{2}\text{\hspace{0.17em}}\mathrm{cos}\left[\frac{2\pi (r-D/2)}{D}\right],$$## 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:

## Eq. (49)

$$\widehat{r}=\u27e8\mathrm{cos}(\mathrm{el})\mathrm{cos}(\mathrm{az}),-\mathrm{cos}(\mathrm{el})\mathrm{sin}(\mathrm{az}),\mathrm{sin}(\mathrm{el})\u27e9,$$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

## Eq. (51)

$$\widehat{b}=\u27e8-\mathrm{sin}(\mathrm{el})\mathrm{cos}(\mathrm{az}),\mathrm{sin}(\mathrm{el})\mathrm{sin}(\mathrm{az}),\mathrm{cos}(\mathrm{el})\u27e9.$$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}\xb7\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

## Eq. (53)

$$\widehat{w}=\u27e8-\mathrm{cos}(\psi )\mathrm{sin}(\mathrm{az})-\mathrm{sin}(\psi )\mathrm{sin}(\mathrm{el})\mathrm{cos}(\mathrm{az}),\phantom{\rule{0ex}{0ex}}-\mathrm{cos}(\psi )\mathrm{cos}(\mathrm{az})+\mathrm{sin}(\psi )\mathrm{sin}(\mathrm{el})\mathrm{sin}(\mathrm{az}),\phantom{\rule{0ex}{0ex}}\mathrm{sin}(\psi )\mathrm{cos}(\mathrm{el})\u27e9.$$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:

## Eq. (54)

$$\text{azimuth}=[360\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}-\mathrm{arctan}2(\frac{{w}_{y}}{{w}_{x}}\left)\right]\%360\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg},$$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

## Eq. (55)

$$\text{azimuth}=\left\{\mathrm{arctan}2\right[\frac{\mathrm{cos}(\delta )\mathrm{sin}(h)}{\mathrm{sin}({\varphi}_{0})\mathrm{cos}(\delta )\mathrm{cos}(h)-\mathrm{cos}({\varphi}_{0})\mathrm{sin}(\delta )}\left]\right\}\%360,$$## Eq. (56)

$$\text{elevation}=\mathrm{arcsin}[\mathrm{sin}({\varphi}_{0})\mathrm{sin}(\delta )+\mathrm{cos}({\varphi}_{0})\mathrm{cos}(\delta )\mathrm{cos}(h)],$$## Acknowledgments

This research was sponsored by grants from NSF AST-1411868, NASA NNX14AJ80G, NNX15AC89G, and NNX15AD95G. Research benefited from the Gemini Observatory, operated by AURA for NSF and the Gemini Consortium. Portions of this work were performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DEAC52-07NA27344. Special thanks are owed to Paul Williams, Alfredo Dubra, Julien Milli, Faustine Cantalloube, and Elena Masciadri for their helpful discussions.