Translator Disclaimer
25 October 2019 Asymmetries in adaptive optics point spread functions
Author Affiliations +

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.



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 lifetimes3 and provide a spectacular view of Kepler’s laws in action. Spectral measurements of individual giant planets constrain evolutionary and atmospheric models4 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.


Theoretical Scintillation Analysis


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 U˜ given by the two-dimensional inverse Fourier transform in the plane (x,y) at some constant z:

Eq. (1)


Direct application of the Helmholtz equation:

Eq. (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 vector k=kxx^+kyy^+kzz^, implying k2=|k|2=kx2+ky2+kz2, we find that propagation along the z axis is constrained by the second-order ordinary differential equation:

Eq. (3)

where we have implicitly included the dependence of U˜ on the particular mode (kx,ky). This differential equation permits solutions of the form:

Eq. (4)

where the ± represents a wave traveling in the positive or negative z direction. For a mode with kx, ky, and the magnitude of the wave vector constrained by |k|=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:

Eq. (5)


Making the substitution kx,y=2πfx,y to represent the true linear wavenumber as in Goodman9 recovers this expression for the free-space propagation transfer function:

Eq. (6)


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 λ2f21), a binomial expansion of the free-space propagation transfer function:

Eq. (7)

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=Aeiϕ with

Eq. (8)


Eq. (9)


It can be shown10 that for small phase errors (α1) the illumination vector at propagation distance z takes the form U(z)=A(z)eiϕ(z), where the phase and amplitude have appropriately been “mixed” due to the scintillation effects:

Eq. (10)


Eq. (11)


Here zT is the Talbot length:

Eq. (12)


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 xxvt, if the wind velocity v points along the positive x axis, and the other the initially measured phase:

Eq. (13)


Eq. (14)


To simplify this expression, we will use the substitutions:

Eq. (15)


Eq. (16)


Eq. (17)

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

Eq. (18)


Eq. (19)


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

Fig. 1

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



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  m1, 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., 

Eq. (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.

Fig. 2

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



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)


However, Goodman9 also claimed that the relationship between the two is the Fourier transform, yet they obtain the expression:

Eq. (22)


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 (kx,ky)(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±i(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 xx and yy, 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 =|F[u]|2 can be expanded in a Taylor series:

Eq. (23)

whose first few terms are

Eq. (24)


Eq. (25)


Eq. (26)


Eq. (27)


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:

Eq. (28)


Eq. (29)

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:

Eq. (30)


Eq. (31)


Eq. (32)


Eq. (33)


Indeed, the even order terms are symmetric and the first odd order term PSF1 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=±1/p:

Eq. (34)


Eq. (35)


Eq. (36)


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

Eq. (37)

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

Fig. 3

The log of the speckle asymmetry ratio as a function of PSF location for parameters α=.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 zT/4 are plotted in black stars (or squares). For these, PSF1=0 because either Asin(2πz/zT)=0 or Pcos(2πz/zT)=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 complicated [set 2πp(Δvt)=π/2+nπ, nZ].


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.


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/zT)].

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.


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

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


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.

Fig. 7

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


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

Fig. 8

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


To further explore our observations of the asymmetry, we plot χ 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 zT 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.

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.




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 r0 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 uncorrelated16 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 coronagraphs17 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 algorithms29,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.

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.


Appendix A: Propagation through the Atmosphere

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

Eq. (38)

where CN2 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 Gavel33 with

Eq. (39)

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 F1 is the unnormalized inverse discrete Fourier transform given by

Eq. (40)

for a discrete array of size P×Q with P,QN. The discrete indices p,a0,1,,P1 and q,b0,1,,Q1 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 1PQ 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:

Eq. (41)

Here x=(x,y) is the coordinate system in the aperture at z=0, Δzi is the range of altitudes relevant to the turbulent layer at altitude zi, and the baseline index of refraction of the atmosphere can be approximated34 with

Eq. (42)

where ρ 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

Eq. (43)

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

Eq. (44)

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)

where v(z) is the wind velocity at altitude z, which is assumed to lie only in the plane at altitude with no vertical component, t0 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:

Eq. (46)

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

Eq. (47)


Here α and β are the coordinates in the image plane, UAO 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 r2=x2+y2:

Eq. (48)

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:

Eq. (49)

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.

Fig. 10

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


Fig. 11

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


Fig. 12

The log of the speckle asymmetry ratio for a single mode with a propagation distance of z=25  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.


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. (50)


Eq. (51)


So that one can think of a^ as pointing in the direction of increasing azimuth and b^ pointing toward increasing elevation. It is left to the reader to show that a^·b^=0, and that a^×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 a^ counterclockwise. Such a wind vector is

Eq. (52)


However, we would instead like to know w^(x^,y^,z^). By algebraically substituting in our coordinate vectors a^ and 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 ψ, az, and el. This is

Eq. (53)


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)

azimuth=[360  degarctan2(wywx)]%360  deg,
where wx, wy, 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

Eq. (55)


Eq. (56)

where h=θLα is the hour angle, θL is the local sidereal time in radians, ϕ0 is the local latitude, α is right ascension, δ is declination, and here we use the convention that azimuth starts at north and opens to the east. The modulo is there to handle overflow and the azimuth and elevation are the coordinates on the sky. Once these are calculated, we can orient images relative to the BT plane because a^ points toward increasing azimuth and b^ points toward increasing elevation.


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.



L. A. Poyneer et al., “Performance of the Gemini planet imager’s adaptive optics system,” Appl. Opt., 55 (2), 323 –340 (2016). APOPAI 0003-6935 Google Scholar


E. L. Nielsen et al., “The Gemini planet imager exoplanet survey: giant planet and brown dwarf demographics from 10–100 au,” Astron. J., 158 (1), 13 (2019). ANJOAA 0004-6256 Google Scholar


J. J. Wang et al., “Dynamical constraints on the HR 8799 planets with GPI,” Astron. J., 156 (5), 192 (2018). ANJOAA 0004-6256 Google Scholar


A. Rajan et al., “Characterizing 51 ERI b from 1 to 5 μm: a partly cloudy exoplanet,” Astron. J., 154 (1), 10 (2017). ANJOAA 0004-6256 Google Scholar


J.-B. Ruffio et al., “Improving and assessing planet sensitivity of the GPI exoplanet survey with a forward model matched filter,” Astrophys. J., 842 (1), 14 (2017). ASJOAB 0004-637X Google Scholar


F. Cantalloube et al., “Status of the MEDUSAE post-processing method to detect circumstellar objects in high-contrast multispectral images,” (2018). Google Scholar


A. Madurowicz et al., “Characterization of Lemniscate atmospheric aberrations in Gemini planet imager data,” Proc. SPIE, 10703 107036E (2018). PSISDG 0277-786X Google Scholar


F. Cantalloube et al., “Origin of the asymmetry of the wind driven halo observed in high-contrast images,” Astron. Astrophys., 620 L10 (2018). AAEJAF 0004-6361 Google Scholar


J. W. Goodman, Introduction to Fourier Optics, McGraw-Hill, New York (1996). Google Scholar


P. Zhou and J. H. Burge, “Analysis of wavefront propagation using the Talbot effect,” Appl. Opt., 49 (28), 5351 –5359 (2010). APOPAI 0003-6935 Google Scholar


E. Hecht, Optics, Addison-Wesley, Reading, Massachusetts (2002). Google Scholar


A. Sivaramakrishnan et al., “Speckle decorrelation and dynamic range in speckle noise-limited imaging,” Astrophys. J., 581 (1), L59 –L62 (2002). ASJOAB 0004-637X Google Scholar


M. D. Perrin et al., “The structure of high STREHL ratio point-spread functions,” Astrophys. J., 596 702 –712 (2003). ASJOAB 0004-637X Google Scholar


NOAA NCEP, “Global forecast system analysis dataset,” Google Scholar


V. P. Bailey et al., “Status and performance of the Gemini planet imager adaptive optics system,” Proc. SPIE, 9909 99090V (2016). PSISDG 0277-786X Google Scholar


O. Guyon, “Limits of adaptive optics for high-contrast imaging,” Astrophys. J., 629 (1), 592 –614 (2005). ASJOAB 0004-637X Google Scholar


C. de Jonge et al., “Effect of multiple deformable mirrors in broadband high-contrast coronagraphs,” Proc. SPIE, 10703 107035D (2018). PSISDG 0277-786X Google Scholar


B.-J. Seo et al., “Hybrid lyot coronagraph for WFIRST: high contrast testbed demonstration in flight-like low flux environment,” Proc. SPIE, 10698 106982P (2018). PSISDG 0277-786X Google Scholar


C. Wu et al., “Phase and amplitude beam shaping with two deformable mirrors implementing input plane and Fourier plane phase modifications,” Appl. Opt., 57 (9), 2337 –2345 (2018). APOPAI 0003-6935 Google Scholar


M. C. Roggemann and D. J. Lee, “Two-deformable-mirror concept for correcting scintillation effects in laser beam projection through the turbulent atmosphere,” Appl. Opt., 37 (21), 4577 –4585 (1998). APOPAI 0003-6935 Google Scholar


B. Macintosh et al., “Extreme adaptive optics for the Thirty Meter Telescope,” Proc. SPIE, 6272 62720N (2006). PSISDG 0277-786X Google Scholar


M. Troy et al., “Effects of diffraction and static wavefront errors on high-contrast imaging from the thirty meter telescope,” Proc. SPIE, 6272 62722C (2006). PSISDG 0277-786X Google Scholar


P. J. Borde and W. A. Traub, “High-contrast imaging from space: speckle nulling in a low-aberration regime,” Astrophys. J., 638 (1), 488 –498 (2006). ASJOAB 0004-637X Google Scholar


A. Give’on, “The electric field conjugation—a unified formalism for wavefront correction algorithms,” in Front. Opt. 2009/Laser Sci. XXV/Fall 2009 OSA Opt. & Photonics Tech. Digest, AOWA3 (2009). Google Scholar


B. L. Gerard, C. Marois and R. Galicher, “Fast coherent differential imaging on ground-based telescopes using the self-coherent camera,” Astron. J., 156 106 (2018). ANJOAA 0004-6256 Google Scholar


B. L. Gerard et al., “Fast focal plane wavefront sensing on ground-based telescopes,” Proc. SPIE, 10703 1070351 (2018). PSISDG 0277-786X Google Scholar


C. Marois, D. W. Phillion and B. Macintosh, “Exoplanet detection with simultaneous spectral differential imaging: effects of out-of-pupil-plane optical aberrations,” Proc. SPIE, 6269 62693M (2006). PSISDG 0277-786X Google Scholar


C. Marois et al., “Angular differential imaging: a powerful high-contrast imaging technique,” Astrophys. J., 641 (1), 556 –564 (2006). ASJOAB 0004-637X Google Scholar


J. J. Wang et al., “pyKLIP: PSF subtraction for exoplanets and disks,” (2015). Google Scholar


L. Pueyo et al., “Application of a damped locally optimized combination of images method to the spectral characterization of faint companions using an integral field spectrograph,” Astrophys. J. Suppl. Ser., 199 (1), 6 (2012). APJSA2 0067-0049 Google Scholar


M. Tallis et al., “Air, telescope, and instrument temperature effects on the Gemini planet imager’s image quality,” Proc. SPIE, 10703 1070356 (2018). PSISDG 0277-786X Google Scholar


V. I. Tatarski, Wave Propagation in a Turbulent Medium, Dover, New York (1961). Google Scholar


E. M. Johansson and D. T. Gavel, “Simulation of stellar speckle imaging,” Proc. SPIE, 2200 372 –383 (1994). PSISDG 0277-786X Google Scholar


J. W. Hardy, Adaptive Optics for Astronomical Telescopes, Oxford University Press, New York (1998). Google Scholar


J. R. Males and O. Guyon, “Ground-based adaptive optics coronagraphic performance under closed-loop predictive control,” J. Astron. Telesc. Instrum. Syst., 4 019001 (2018). Google Scholar


R. B. Blackman and J. W. Tukey, “The measurement of power spectra from the point of view of communications engineering—part I,” Bell Syst. Tech. J., 37 (1), 185 –282 (1958). BSTJAN 0005-8580 Google Scholar

Biographies of the authors are not available.

© The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.


Airborne laser advanced concepts testbed
Proceedings of SPIE (August 02 1999)
Single star SCIDAR for Cn2(h) profiling
Proceedings of SPIE (August 17 2005)
Atmospheric effects on active illumination
Proceedings of SPIE (August 22 2005)

Back to Top