Asymmetries in adaptive optics point spread functions

An explanation for the origin of asymmetry along the preferential axis of the 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 towards 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.


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

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Ũ given by the two-dimensional inverse Fourier transform in the plane ðx; yÞ at some constant z: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 7 3 0 Uðx; y; zÞ ¼ ZZ þ∞ −∞Ũ ðk x ; k y ; zÞe iðk x xþk y yÞ dk x dk y : (1) Direct application of the Helmholtz equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 6 6 2 ð∇ 2 þ k 2 ÞUðx; y; zÞ ¼ 0 (2) 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 vectork ¼ k xx þ k yŷ þ k zẑ , implying k 2 ¼ jkj 2 ¼ k 2 x þ k 2 y þ k 2 z , we find that propagation along the z axis is constrained by the second-order ordinary differential equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 5 5 7 d 2 dz 2Ũ ðzÞ þ k 2 zŨ ðzÞ ¼ 0; where we have implicitly included the dependence ofŨ on the particular mode ðk x ; k y Þ. This differential equation permits solutions of the form: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 4 7 9Ũ ðzÞ ¼Ũðz ¼ 0Þe AEik z z ; (4) where the AE 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 jkj ¼ 2π∕λ, where λ 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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 3 2 6 ; 3 5 8 Hðk x ; k y ; zÞ ¼ exp AEi Making the substitution k x;y ¼ 2πf x;y to represent the true linear wavenumber as in Goodman 9 recovers this expression for the free-space propagation transfer function: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 3 2 6 ; 2 7 4 Hðf x ; f y ; zÞ ¼ exp 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.

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 λ 2 f 2 ≪ 1), a binomial expansion of the free-space propagation transfer function: Journal of Astronomical Telescopes, Instruments, and Systems 049003-2 Oct-Dec 2019 • Vol. 5(4) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 6 3 ; 7 5 2 HðzÞ ≈ e −iπλ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 α, the complex illumination can be written U ¼ Ae iϕ with E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 6 3 ; 6 7 3 A ¼ 1; E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 6 3 ; 6 3 1 ϕ ¼ α sinð2πx∕pÞ: It can be shown 10 that for small phase errors (α ≪ 1) the illumination vector at propagation distance z takes the form UðzÞ ¼ AðzÞe iϕðzÞ , where the phase and amplitude have appropriately been "mixed" due to the scintillation effects: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 6 3 ; 5 7 6 AðzÞ ¼ 1 þ α sinð2πz∕z T Þ sinð2πx∕pÞ; E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 6 3 ; 5 3 4 ϕðzÞ ¼ α cosð2πz∕z T Þ sinð2πx∕pÞ: (11) Here z T is the Talbot length: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 6 3 ; 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 → x − vt, if the wind velocity v points along the positive x axis, and the other the initially measured phase: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 6 3 ; 3 9 0 ϕ AO ¼ ϕ z ðx − vtÞ − ϕ z ðxÞ; E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 3 2 6 ; 7 5 2 To simplify this expression, we will use the substitutions: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 3 2 6 ; 6 8 7 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 3 2 6 ; 6 4 5 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 3 2 6 ; 6 2 2 where A represents the term that modulates the amplitude, P is the term that modulates the phase, and Δ 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 3 2 6 ; 5 4 3 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 3 2 6 ; 5 0 1 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.  An example of a single mode of AO-corrected complex illumination and corresponding asymmetric speckles in the PSF. Here α ¼ :01 rad, z ¼ 10 km, p ¼ 1 m, v ¼ 10 m∕s, and t ¼ 3.2 ms. The AOcorrected 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. 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 m −1 , these valid temporal frequencies in our model correspond to wind velocities in the range 8 to 72 m∕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 m∕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., E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 6 3 ; 2 6 3 ϕ AO ðtÞ ¼ ϕðtÞ − 0.997ϕðt − τÞ; (20) where τ is the delay time. The open-loop model does not agree at higher temporal frequencies, most obviously when temporal frequency is the inverse of τ. 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.

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: Aðx; yÞ exp½iðk x x þ k y yÞdx dy: (21) However, Goodman 9 also claimed that the relationship between the two is the Fourier transform, yet they obtain the expression: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 3 2 6 ; 6 7 1 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Þ → ðξ; ηÞ and in the image plane ðk x ; k y Þ → ð 2πx λz ; 2πy λ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 AEiðkr−ω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 → −x and y → −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.

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 ϕ, the intensity distribution or PSF ¼ jF ½uj 2 can be expanded in a Taylor series: whose first few terms are E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 4 ; 6 3 ; 7 0 8 PSF 0 ¼ aa Ã ; E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 5 ; 6 3 ; 6 7 6 PSF 1 ¼ 2 Im½aða Ã ⋆Φ Ã Þ; E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 6 ; 6 3 ; 6 5 5 PSF 2;halo ¼ ða⋆ΦÞða Ã ⋆Φ Ã Þ; E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 7 ; 6 3 ; 6 3 4 PSF 2;strehl ¼ − Here, the case change is used as shorthand for the Fourier transform, so a ¼ F ðAÞ and Φ ¼ F ðϕÞ. Additionally, ⋆ is the convolution operation and Ã is the complex conjugate. For our AO-corrected illumination, these can be calculated as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 8 ; 6 3 ; 5 6 9 a AO ¼ δðfÞ þ E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 9 ; 6 3 ; 5 0 where δ 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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 0 ; 6 3 ; 4 2 8 PSF 0 ¼ δðfÞ þ E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 1 ; 6 3 ; 3 7 4 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 2 ; 6 3 ; 3 3 4 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 3 ; 6 3 ; 2 4 3 Indeed, the even order terms are symmetric and the first odd order term 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 ¼ AE1∕p: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 4 ; 6 3 ; 1 4 7 PSF 0 ðf ¼ −1∕pÞ E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 5 ; 6 3 ; 9 3 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 6 ; 3 2 6 ; 7 4 1 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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 7 ; 3 2 6 ; 6 5 7 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 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. Fig. 3 The log of the speckle asymmetry ratio as a function of PSF location for parameters α ¼ :01 rad, t ¼ 3.2 ms, and z ¼ 10 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 λ ¼ 1.6 μ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, PSF 1 ¼ 0 because either A ∼ sinð2πz∕z T Þ ¼ 0 or P ∼ cosð2π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 m∕s case, though their analytic solution is more compli- Journal of Astronomical Telescopes, Instruments, and Systems 049003-5 Oct-Dec 2019 • Vol. 5(4)

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 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π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 log χ ¼ 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 Fig. 4 For a single layer of atmospheric turbulence at z ¼ 10 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.
Journal of Astronomical Telescopes, Instruments, and Systems 049003-6 Oct-Dec 2019 • Vol. 5 (4) 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 m∕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 m∕s, the scintillation becomes stronger and introduces a noticeable asymmetry in the PSF along the axis of the wind propagation.

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 χ between the two lobes of the PSF. χ 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 χ is folded around 1 by inversion such that the half-annulus with greater intensity is always in the numerator, meaning χ always takes on values greater than 1, and larger values correspond to greater asymmetry. Note that this χ is not identical to the χ 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 Fig. 6 A selection of images taken with the Gemini Planet Imager, ordered by decreasing asymmetry metric χ. 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" × 2.8" 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 χ. Glancing over the entire dataset, it is readily apparent that the strong asymmetry direction is predominately pointing into the west.
Journal of Astronomical Telescopes, Instruments, and Systems 049003-7 Oct-Dec 2019 • Vol. 5 (4) 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 χ ≥ 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 χ 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 (elevation ¼ 90 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 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.

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   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.
Journal of Astronomical Telescopes, Instruments, and Systems 049003-8 Oct-Dec 2019 • Vol. 5 (4) 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. Spacebased 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 log χ ¼ 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. Fig. 9 Asymmetry strength χ 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.
Journal of Astronomical Telescopes, Instruments, and Systems 049003-9 Oct-Dec 2019 • Vol. 5 (4) 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.
where C 2 N is the index of refraction structure constant and κ ¼ 2π∕l is the spatial wave vector for an eddy of size l. Here, we use the standard Kolmogorov power spectrum, which is fractally self-similar at all length scales, although it is in principle simple to extend this model to a Von-Karman spectrum by attenuating the power above and below the outer and inner scales. From the square root of the power spectrum, we can find the fluctuations from the inverse Fourier transform according to Johansson and Gavel 33 with E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 9 ; 6 3 ; 2 6 6 δNðx; zÞ ¼ Re where δN are the fluctuations of the index of refraction from unity in parts per million, ξ is a zero-mean unit-variance complex hermitian Gaussian noise process, and F −1 is the unnormalized inverse discrete Fourier transform given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 0 ; 6 3 ; 1 7 7 η ab ¼ F −1 ðη pq Þ ¼ for a discrete array of size P × Q with P; Q ∈ N. The discrete indices p; a ∈ 0;1; : : : ; P − 1 and q; b ∈ 0;1; : : : ; Q − 1 exist in Fourier and configuration space, respectively. The corresponding forward Fourier transform simply includes negation in the exponent, and we have to pay careful attention the normalization factor used by a routine such as np.fft.fft2, which includes a normalization of 1 PQ on the inverse transform, but no normalization on the forward transform by default.
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π∕λ times the optical path length: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 1 ; 3 2 6 ; 6 5 1 Herex ¼ ðx; yÞ is the coordinate system in the aperture at z ¼ 0, where ρ 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 sec ζ where ζ 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 4 ; 3 2 6 ; 1 2 6 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 Journal of Astronomical Telescopes, Instruments, and Systems 049003-10 Oct-Dec 2019 • Vol. 5 (4) 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 5 ; 6 3 ; 6 8 6 δN½x þṽðzÞτ; t 0 þ τ ¼ δNðx; t 0 Þ; whereṽðzÞ is the wind velocity at altitude z, which is assumed to lie only in the plane at altitude with no vertical component, t 0 is a particular instant in time, and τ is the total time delay for the AO system to respond to a measurement from the wavefront sensor. We also assume a perfect noiseless wavefront sensor and deformable mirror with only a time lag or servo-lag error for an ideal open-loop AO simulation. The expression for the compensated phase in the aperture is a new complex illumination with the amplitude errors from the current timestep and the phase given by two subtracted phases, one from the current timestep and one from the previous, which is our AO correction: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 6 ; 6 3 ; 5 3 4 ϕ AO ¼ ϕðx; t 0 Þ − ϕðx; t 0 − τÞ; where we implicitly have included the contributions from L turbulent layers at various altitudes z with a flat interpolation scheme for the structure constant. From the compensated phase on the aperture, we can obtain the final image's intensity distribution with an inverse Fourier transform by assuming the telescope focus operates in a Fraunhofer diffraction limit, so that electric field distribution in the image plane is the inverse Fourier transform of the aperture function: 11 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 7 ; 6 3 ; 4 1 4 Iðα; βÞ ¼ hjF −1 ðAU AO Þj 2 i: Here α and β are the coordinates in the image plane, U AO is our AO corrected complex illumination which includes amplitude errors from scintillation, A is the aperture function with the Blackman window apodization, 36 parameterized radially from the center with r 2 ¼ x 2 þ y 2 : E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 8 ; 6 3 ; 3 2 7 with an aperture diameter of D ¼ 8 m, a falloff of γ ¼ :16, and the brackets denote time average over the whole length of the simulation. The Blackman apodization simulates a crude coronagraph and dampens the high-order airy rings, whose final intensity in the image plane can swamp the effect of the scintillation halo.

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: where we have assumed the convention of the positive x axis pointing north, and the positive y axis pointing west. This conveniently sets up the positive z axis to point toward zenith, as it should. azimuth is measured from north opening toward the east, and elevation is measured from the horizon upward. See Fig. 10 for an illustration. It is worth noting that the validity of this analogy, as well as is necessary to implement angular differential imaging, a postprocessing technique for combining multiple exposures while the target star moves through the zenith that GPI operates in a fixed parallactic orientation, with the instrument derotator disabled, so that GPI is fixed with respect to the telescope orientation, which is uncommon.
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â as pointing in the direction of increasing azimuth andb pointing toward increasing elevation. It is left to the reader to show thatâ ·b ¼ 0, and thatâ ×r ¼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 ψ fromâ counterclockwise. Such a wind vector is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 5 2 ; 3 2 6 ; 2 1 9ŵ ¼ cosðψÞâ þ sinðψÞb: (52) However, we would instead like to knowŵðx;ŷ;ẑÞ. By algebraically substituting in our coordinate vectorsâ andb formulas in x, y, z space, we can arrive at an expression for the wind vector in x, y, z space in terms of ψ, az, and el. This is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 5 3 ; 3 2 6 ; 1 3 2ŵ ¼ − cosðψÞ sinðazÞ − sinðψÞ sinðelÞ cosðazÞ; − cosðψÞ cosðazÞ þ sinðψÞ sinðelÞ sinðazÞ; sinðψÞ cosðelÞ : (53) Fig. 10 Diagram of coordinates used to orient images on the sky demonstrating the relationship between the ground plane and the BT plane.
Journal of Astronomical Telescopes, Instruments, and Systems 049003-11 Oct-Dec 2019 • Vol. 5 (4) 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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 5 4 ; 6 3 ; 1 3 8 azimuth ¼ 360 deg − arctan 2 w y w x %360 deg; where w x , w y , are the x and y components of the wind vector, respectively, % is the modulo operator, and it is often convenient to use a smart operator like arctan2 to get the quadrant correct. 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 Journal of Astronomical Telescopes, Instruments, and Systems 049003-12 Oct-Dec 2019 • Vol. 5(4) Fig. 12 The log of the speckle asymmetry ratio for a single mode with a propagation distance of z ¼ 25 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 m∕s, which demonstrates the residual uncorrected amplitude error.