Asymmetries in adaptive optics point spread functions
Alexander Madurowicz Bruce Macintosh Vanessa P. Bailey Jeffrey Chilcote Marshall Perrin Lisa Poyneer Laurent Pueyo Jean-Baptiste Ruffio Travis Barman Joanna Bulger Tara Cotten Robert J. De Rosa Rene Doyon Gaspard Duchêne Thomas M. Esposito Michael P. Fitzgerald Katherine B. Follette Benjamin L. Gerard Stephen J. Goodsell James R. Graham Alexandra Z. Greenbaum Pascale Hibon Li-Wei Hung Patrick Ingraham Paul Kalas Quinn Konopacky Jérôme Maire Franck Marchis Mark S. Marley Christian Marois Stanimir Metchev Maxwell A. Millar-Blanchaer Eric L. Nielsen Rebecca Oppenheimer David Palmer Jennifer Patience Abhijith Rajan Julien Rameau Fredrik T. Rantakyrö Dmitry Savransky Anand Sivaramakrishnan Inseok Song Remi Soummer Melissa Tallis Sandrine Thomas Jason J. Wang Kimberly Ward-Duong Schuyler Wolf

Coupled Cluster Green's function formulations based on the effective  Hamiltonians
    Karol Kowalski
    December 30, 2023
==================================================================================




    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 is one such instrument <cit.>, where fast and precise correction are pivotal to optimize instrument performance. The results so far have been spectacular. With detections of multiple planetary mass companions around various stars, as well a strong non-detection limits around many more, one can constrain planetary population distributions.<cit.> Astrometric measurements of multi-body systems probe dynamical constraints on planetary masses and system lifetimes <cit.>, and provide a spectacular view of Kepler's Laws in action. Spectral measurements of individual giant planets constrain evolutionary and atmospheric models <cit.> of these objects, paving the way towards 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 development of algorithms optimized in the search for planets <cit.><cit.>. The analysis in this paper expands on our previous work <cit.>, 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. 2018 <cit.>. We will expand on their discussion by 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 non-linear. 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 the 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 which 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 Adaptive Optics. The third section demonstrates the effects from our analytic model is 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, the paper concludes with a brief discussion about the importance of this effect in the context of improving Adaptive Optics systems performance from design to post-processing.



§ 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

    U(x, y, z) = ∬_-∞^+∞U(k_x, k_y, z) e^i (k_x x + k_y y)dk_xdk_y .

Direct application of the Helmholtz equation

    (∇^2 + k^2)U(x,y,z) = 0

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 wavevector k⃗ = k_x x̂ + k_y ŷ + k_zẑ, implying k^2 = |k⃗|^2 = k_x^2 + k_y^2 + k_z^2, we find that propagation along the z-axis is constrained by the second-order ordinary differential equation

    d^2/dz^2U(z) + k_z^2U(z) = 0 ,

where we have implicitly included the dependence of U on the particular mode (k_x, k_y). This differential equation permits solutions of the form

    U(z) = U(z = 0) e^± i k_z z,

where the ± represents a wave traveling in the positive or negative z direction. For a mode with k_x, k_y, and the magnitude of the wavevector 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

    H(k_x, k_y, z) = exp[ ± i √((2π/λ)^2 - (k_x^2 + k_y^2)) z ] .

Making the substitution k_x,y = 2π f_x,y to represent the true linear wavenumber as in Goodman <cit.> recovers this expression for the free space propagation transfer function

    H(f_x, f_y, z) = exp[ ± i 2π/λ√(1 - λ^2(f_x^2 + f_y^2)) z ] .

This allows us to consider the propagation of light through free space as a linear spatial filter applied to each Fourier mode of the complex illumination independently. 



 §.§ 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 λ^2f^2 << 1), a binomial expansion of the free space propagation transfer function

    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

    A    = 1 
    ϕ   = αsin(2π x/p) .

It can be shown <cit.> that for small phase errors (α << 1) that 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

    A(z)    = 1 + αsin(2π z/z_T)sin(2π x/p) 
    ϕ(z)    = αcos(2π z/z_T)sin(2π x/p) .

Here z_T is the Talbot Length

    z_T = λ/1 - √(1 - λ^2f^2)≈2/λ f^2 .


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

    ϕ_AO   =   ϕ_z(x-vt) - ϕ_z(x) 
       =   αcos(2π z/z_T)√(2)√(1-cos(2π vt/p))sin(2π x/p + arctan[-sin(2π vt/p)/cos(2π vt/p) - 1]) .

To simplify this expression, we will use the substitutions

    𝒜   =   αsin(2π z/z_T) 
    
    P    =   αcos(2π z/z_T)√(2)√(1-cos(2π vt/p))
    Δ   =    -p/2πarctan[-sin(2π vt/p)/cos(2π vt/p) - 1]

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

    A_AO   =    1 + 𝒜sin(2π(x-vt)/p) 
    ϕ_AO   =    Psin(2π(x-Δ)/p).

An example of a single mode of AO-corrected illumination and its corresponding electric field distribution in the image plane is demonstrated in Figure 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.



 §.§ 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 section 4.D of <cit.> 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 Figure 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 figure, was fit by adjusting the effective delay such the measurement is applied 3.2 ms after it was taken.  As shown in the 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 m/s 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 m/s 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. 

    ϕ_AO(t) = ϕ(t) - 0.997ϕ(t - τ),

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

    E(k_x, k_y) = ∬_-∞^+∞𝒜(x,y)exp[i(k_x x + k_y y)]dxdy .

However, Goodman <cit.> also claims that the relationship between the two is the Fourier transform, yet they obtain the expression

    U(x,y) = e^ikze^ik/2z(x^2+y^2)/iλ z∬_-∞^∞ U(ξ,η) exp[ -i2π/λ z(xξ + yη)]dξdη.

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^± 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 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



While 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. <cit.> and Perrin et al. <cit.> 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 = |ℱ[u]|^2  can be expanded in a Taylor series 

    PSF   ≈   PSF_0 + PSF_1 + PSF_2,halo + PSF_2,strehl,

whose first few terms are

    PSF_0    =    aa^* 
    PSF_1    =    2Im[a(a^* ⋆Φ^*)] 
    PSF_2,halo   =    (a ⋆Φ)(a^* ⋆Φ^*) 
    PSF_2,strehl   =    -1/2[a(a^* ⋆Φ^* ⋆Φ^*) + a^*(a ⋆Φ⋆Φ)] .

Here the case change is used as shorthand for the Fourier transform, so a = ℱ[A] and Φ = ℱ[ϕ]. Additionally, ⋆ is the convolution operation and * is the complex conjugate. For our AO-corrected illumination, these can be calculated as follows

    a_AO = δ(f) +𝒜/2iδ(f-1/p)e^-i2π/pvt-𝒜/2iδ(f+1/p)e^i2π/pvt


    Φ_AO = P/2iδ(f-1/p)e^-i2π/pΔ-P/2iδ(f+1/p)e^i2π/pΔ,

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

    PSF_0 = δ(f) + 𝒜^2/4[δ(f-1/p)+δ(f+1/p)]    

    PSF_1 =  -𝒜P/2sin(2π/p(vt-Δ))[δ(f-1/p)-δ(f+1/p) ]    

    PSF_2,halo = P^2/4[δ(f+1/p) + δ(f-1/p)]  + ^2 P^2/16[2(1 + cos(4π/p(vt -Δ)))δ(f)
    + δ(f + 2/p) + δ(f - 2/p)]     

    PSF_2,strehl = -P^2/2δ(f) -𝒜^2P^2/32(4 + 2cos(4π/p(vt-Δ)))[δ(f-1/p)+δ(f+1/p)].   

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=±1/p

    PSF_0(f=-1/p) = PSF_0(f=1/p)    = 𝒜^2/4
    
    - PSF_1(f=-1/p) = PSF_1(f=1/p)    =  -𝒜P/2sin(2π/p(vt-Δ))
    PSF_2(f=-1/p) = PSF_2(f=1/p)    = P^2/4 - 𝒜^2P^2/32[4 + 2cos(4π/p(vt-Δ))] .

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

    χ = PSF(f=1/p)/PSF(f=-1/p) ,

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



§ 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 many of the precise details of these simulations are referred to Appendix A, 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 Figure 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 multi-layer case, the wind velocity for the jet stream is generally pointing towards 2 o'clock. 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. While at 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 upwards 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 a angular profile of the image by integrating along the radial axis. This allows us to fit for the preferential wind direction modulo 180 degrees. These techniques follow exactly from our previous work <cit.>, 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-annuli perpendicular to the axis of wind propagation. This asymmetry metric is used to break the 180 degree 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 Figure 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 Appendix B. The directions and velocities of the wind for various layers in the Earth's atmosphere are taken from the NOAA Global Forecast System <cit.>, 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 Figure 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 towards the West. This observation confirms our prior analysis regarding the strength of the asymmetry being opposed to the direction of the wind.


To further verify this correlation between the jet stream and the asymmetric butterfly, we take our sample of data where χ≥ 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 Figure A2. 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-300 hPa, which roughly corresponds to altitudes of 10-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 Figure 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 Figure 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 towards 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 towards 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 Figure A3 for a remake of Figure 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 <cit.>. 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 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 <cit.>, 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 <cit.>, and here we have demonstrated severity of this effect, visible in the form of PSF asymmetry. This has implications for both current and future adaptive optics systems, especially but not exclusively those designed for high-contrast imaging. Many analyses often assume that various sources of scattered light are uncorrelated <cit.> 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. While 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 adaptive optics 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 Figure A4). This in turn will set a performance floor for such systems, and the asymmetry detected here provides a first measurement of the level at 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 our 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 <cit.> to correct static amplitude errors, and preliminary laboratory testing is underway <cit.>. Similar designs have been expressed for improvements in laser communications <cit.><cit.>. 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 re-coated at different times <cit.><cit.>. 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 <cit.><cit.>. Another similar focal plane wavefront sensing approach was recently proposed on ground-based telescopes <cit.>, 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 <cit.>, or something as simple as adding a high-speed direct pupil-imaging channel to a traditional adaptive optics wavefront sensor. 

The scintillation halo observed in AO PSFs is a challenge for effective post-processing 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, non-common-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 (IFS) 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 <cit.>. 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. <cit.>

Various different post processing algorithms <cit.><cit.> 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. <cit.> In addition, imaging extended objects such as debris disks cannot utilize a HPF in post-processing, 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. <cit.> 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 <cit.>, 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.

plain







§ APPENDIX A: PROPAGATION THROUGH THE ATMOSPHERE

Tartarski<cit.> has shown that the fluctuations in the optical index of refraction in three dimensions for a Kolmogorov turbulence spectrum follow the form

    Φ_N(κ, z) = 0.033 C_N^2(z)κ^-11/3

where C_N^2 is the index of refraction structure constant and and κ = 2π/l is the spatial wavevector 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<cit.> with


    δ N (x⃗,z) = Re[ℱ^-1(ξ(κ⃗,z)√(Φ_N(κ,z)))]

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

    η_ab = ℱ^-1(η̃_pq) = ∑_p=0^P-1∑_q=0^Q-1η̃_pqexp[2π i(pa/P + qb/Q)]

for a discrete array of size P× Q with P,Q ∈ℕ. 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 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 wavevector of the ray k = 2π/λ times the optical path length.

    ϕ_i(x⃗) = k∫_z_i^z_i+Δ z_i n(x⃗,z) dz = k n(x⃗,z)Δ z_i

Here x⃗ = (x, y) is the coordinate system in the aperture at z=0, Δ z_i is the range of altitudes relevant to the turbulent layer at altitude z_i, and the baseline index of refraction of the atmosphere can be approximated<cit.> with 

    N ≡ (n-1)10^6 ≈ 77.6 ρ/T

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


    ρ(z) = ρ_0 e^-Mgz/RT

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 non-zenith observations an additional term of ζ where ζ is the zenith angle should be included in the integral in (40) 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.<cit.>

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


    u(z) = ℱ^-1[H(z) ℱ[u]]

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

    δ N(x⃗+v⃗(z)τ,t_0+τ) = δ N(x⃗, t_0)

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, t_0 is an particular instant in time, and τ is the total time delay for the adaptive optics 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 of 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

    ϕ_AO = ϕ(x⃗,t_0) - ϕ(x⃗,t_0-τ)

where we implicitly have included the 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<cit.>.


    I(α, β) = | ⟨ℱ^-1(𝒜U_AO)⟩ |^2

Here α, β are the coordinates in the image plane, U_AO is our AO corrected complex illumination which includes amplitude errors from scintillation, 𝒜 is the aperture function with the Blackman window apodization <cit.>, parameterized radially from the center with r^2 = x^2 + y^2


    𝒜(r) = 1 - γ/2 - 1/2cos(2π(r-D/2)/D) + γ/2cos(2π(r-D/2/D)

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: THE 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 multi-million dollar telescope with an Alt-Az tracking system. Either way[It is worth noting that the validity of this analogy, as well as is necessary to implement Angular Differential Imaging, a post-processing 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.], your imaging device is pointed at the celestial sphere along the line of sight vector 

    r̂ = ⟨cos(el)cos(az), -cos(el)sin(az), sin(el)⟩

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 towards Zenith, as it should. Azimuth is measured from North opening towards the East, and elevation is measured from the horizon upwards. See Figure A1 for an illustration.



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

    â   =   ⟨ - sin(az), - cos(az), 0 ⟩
    b̂   =   ⟨-sin(el)cos(az), sin(el)sin(az), cos(el)⟩

So that one can think of â as pointing in the direction of increasing Azimuth, and b̂ pointing towards 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

    ŵ = cos(ψ)â + sin(ψ)b̂

However, we would instead like to know ŵ(x̂,ŷ,ẑ). By algebraically substituting in our coordinate vectors â, 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

    ŵ = ⟨ -cos(ψ)sin(az) - sin(ψ)sin(el)cos(az), 
    
    -cos(ψ)cos(az) + sin(ψ)sin(el)sin(az), 
    sin(ψ)cos(el)⟩

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

    azimuth =(360 - arctan 2(w_y/w_x)) % 360

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


    azimuth   =   (arctan2(cos(δ)sin(h)/sin(ϕ_0)cos(δ)cos(h) - cos(ϕ_0)sin(δ)) + 180) % 360 
    elevation   =   arcsin(sin(ϕ_0)sin(δ) + cos(ϕ_0)cos(δ)cos(h)

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 â points towards increasing azimuth and b̂ points towards increasing elevation.



§ APPENDIX C: CALCULATING THE 2ND ORDER PSF

For some electric field U incident upon an aperture, described by

    U = A e^iφ

and defining the change of case to represent the Fourier Transform

    a    ≡    [A] 
    
        Φ    ≡    [φ] 
     
        u    ≡    [U]

The PSF formed by such an electric field can be calculated with (Sivaramakrishnan et al. and Perrin et al.)

    =    |u|^2 = uu^* 
       ≈   _0 + _1 + _2

Where the first few terms of the Taylor series are

    _0    =    aa^* 
    _1    =    2 Im[a(a^* ⋆ Φ^*)] 
    _2    =    (a⋆ Φ)(a^*⋆ Φ^*) - 1/2[a(a^*^*^*) + a^*(a)]

In our one dimensional analytic model, the amplitude and phase that result after a time-delayed AO phase correction of an incoming scintillated wavefront are

    A    =    1 + sin(2π(x-vt)/p) 
    
        φ    =    Psin(2π(x-Δ)/p)

where we have used the shorthand notations

    =    αsin(2πz/z_T) 
    
        P    =    αcos(2πz/z_T)√(2)√(1 - cos(2πvt/p))
    
        Δ    =    -p/2πarctan[-sin(2πvt/p)/cos(2πvt/p) - 1]

to simplify the expressions. In the following sections, we calculate a, Φ, and then use them to obtain PSF_0, PSF_1, PSF_2.



 §.§ Calculating a

Starting with the definition of a and the definition of the Fourier Transform

    a =  [A] = ∫_-∞^∞ [1 + sin(2π/p(x-vt))]e^-i2πfxx̣

and using the property that sin(x) = 1/2i(e^ix - e^-ix)

    a    =   ∫_-∞^∞[1 + /2i(e^i2π/p(x-vt) - e^-i2π/p(x-vt))]e^-i2πfxx̣
    
        a    =   ∫_-∞^∞[e^-2πfx + /2i( e^i2πx(1/p - f)e^-i2π/pvt -  e^-i2πx(1/p + f)e^i2π/pvt)]x̣

along with the identity that δ(x-α) = ∫_-∞^∞ e^i2πf(x-α)f̣

    a = δ(f) + /2iδ(f-1/p)e^-i2π/pvt - /2iδ(f+1/p)e^i2π/pvt




 §.§ Calculating Φ

Likewise

    Φ =  [φ]    =   ∫_-∞^∞ Psin(2π/p(x-Δ))e^-i2πfxx̣
       =   ∫_-∞^∞P/2i(e^i2π/p(x-Δ) - e^-i2π/p(x-Δ))e^-i2πfxx̣
       =   ∫_-∞^∞P/2ie^i2πx(1/p - f)e^-i2π/pΔ -P/2ie^-i2πx(1/p - f)e^i2π/pΔx̣
       =   P/2iδ(f-1/p)e^-i2π/pΔ -P/2iδ(f+1/p)e^i2π/pΔ




 §.§ Calculating PSF_0

Since _0 = aa^*,

    _0 = [δ(f) + /2iδ(f-1/p)e^-i2π/pvt - /2iδ(f+1/p)e^i2π/pvt] ×
    [ δ(f) + /-2iδ(f-1/p)e^i2π/pvt - /-2iδ(f+1/p)e^-i2π/pvt]

And using the mathematically non-rigorous (though conceptually clear) notions that [δ(f-α)]^2 = δ(f-α) and δ(f-α)δ(f-β) = 0 when α≠β, then the cross terms cancel and the complex terms go to their moduli squared and the equation simplifies to

    _0 = δ(f) + ^2/4[δ(f-1/p) + δ(f+1/p)].




 §.§ Calculating PSF_1

Starting with the definition of PSF_1

    _1 = 2Im[a(a^*^*)],

and the definition of the convolution operation

    (a^* ⋆ Φ^*)(f) = ∫_-∞^∞ a^*(f')Φ^*(f-f')f̣' ,


    (a^* ⋆ Φ^*)(f) = ∫_-∞^∞[ δ(f') - /2iδ(f'-1/p)e^i2π/pvt + /2iδ(f'+1/p)e^-i2π/pvt] ×
    [P/-2iδ((f-f')-1/p)e^i2π/pΔ +P/2iδ((f-f')+1/p)e^-i2π/pΔ]

and using the sifting property ∫_-∞^∞ f(t) δ(t-τ)ṭ = f(τ) six times, once for each term,

    ∫_-∞^∞ δ(f') P/-2iδ((f-f')-1/p)e^i2π/pΔf̣'    = P/-2ie^i2π/pΔδ(f-1/p) 
    ∫_-∞^∞ δ(f') P/2iδ((f-f')+1/p)e^-i2π/pΔf̣'    = P/2ie^-i2π/pΔδ(f+1/p) 
    ∫_-∞^∞/-2iδ(f'-1/p)e^i2π/pvtP/-2iδ((f-f')-1/p)e^i2π/pΔf̣'    = - P/4e^i2π/p(vt +Δ)δ(f-2/p) 
    ∫_-∞^∞/-2iδ(f'-1/p)e^i2π/pvtP/2iδ((f-f')+1/p)e^-i2π/pΔf̣'    =  P/4e^i2π/p(vt -Δ) δ(-f) 
    ∫_-∞^∞/2iδ(f'+1/p)e^-i2π/pvtP/-2iδ((f-f')-1/p)e^i2π/pΔf̣'    =  P/4e^-i2π/p(vt -Δ) δ(f) 
    ∫_-∞^∞/2iδ(f'+1/p)e^-i2π/pvtP/2iδ((f-f')+1/p)e^-i2π/pΔf̣'    = - P/4e^-i2π/p(vt +Δ) δ(f + 2/p) .

Here it is worth noting that the δ function is even and so obeys the relationship that δ(-x) = δ(x). Furthermore, factoring out common terms

    (a^* ⋆ Φ^*) =  P/4[e^i2π/p(vt -Δ) δ(f)    -e^i2π/p(vt +Δ)δ(f-2/p) 
    
        e^-i2π/p(vt -Δ) δ(f)    -e^-i2π/p(vt +Δ)δ(f + 2/p) ] 
    
        + P/2i[e^-i2π/pΔδ(f+1/p)    - e^i2π/pΔδ(f-1/p)] .

Multiplying by a, only the cross terms with δ(f), δ(f-1/p), and δ(f+1/p) remain

    (a(a^* ⋆ Φ^*))(f) =  P/4[e^i2π/p(vt -Δ) δ(f) 
         + e^-i2π/p(vt -Δ) δ(f) ] 
    
         +  P/4 i^2[ -δ(f-1/p)e^-i2π/pvte^i2π/pΔ - δ(f+1/p)e^i2π/pvte^-i2π/pΔ]

factoring again,

    (a(a^* ⋆ Φ^*))(f) =  P/4[e^i2π/p(vt -Δ)(δ(f) + d(f+1/p)) + e^-i2π/p(vt -Δ)(δ(f) + δ(f-1/p)) ]

And taking twice the imaginary part, using Im[e^iαx] = sin(αx)

    2Im[(a(a^* ⋆ Φ^*))](f) =  P/2[sin(2π/p(vt -Δ))(δ(f) + δ(f+1/p)) 
    
        + sin(-2π/p(vt -Δ))(δ(f) + δ(f-1/p)) ]

However, since sine is odd and sin(-x) = -sin(x), the term with δ(f) cancels and

    _1 = 2Im[(a(a^* ⋆ Φ^*))](f) =  P/2sin(2π/p(vt -Δ)) [δ(f+1/p) - δ(f-1/p)]

and as we can see from the minus sign it is antisymmetric.



 §.§ Calculating PSF_2

Since

    _2 = (a⋆ Φ)(a^*⋆ Φ^*) - 1/2[a(a^*^*^*) + a^*(a)],

we start first by calculating the first term _2,halo = (a⋆ Φ)(a^*⋆ Φ^*). Since we have already calculated (a^*⋆ Φ^*) in the previous section (see equation <ref>), we only need to find (a ⋆ Φ)  . While, in general (g^*⋆ h^*) ≠ (g⋆ h)^* for generic complex functions g and h, here the complex conjugate does not modify the δ function, which has no imaginary part, and therefore does not change the convolution integral, as it is the only part of the function dependent on the variable of integration f'. Since the complex conjugate only modifies the multiplicative constants with i and the phasors e^iα, we can simply take the complex conjugate of (<ref>) and arrive at the answer.

    (a⋆ Φ) =  P/4[e^-i2π/p(vt -Δ) δ(f)    -e^-i2π/p(vt +Δ)δ(f-2/p) 
    
        e^+i2π/p(vt -Δ) δ(f)    -e^+i2π/p(vt +Δ)δ(f + 2/p) ] 
    
        + P/-2i[e^+i2π/pΔδ(f+1/p)    - e^-i2π/pΔδ(f-1/p)] .

This convenient fact means that when we multiply the two terms together, nearly all of the complex phasors will go to their moduli squared, aka the constant 1, simplifying the expression very nicely. However, since cross terms from both of the δ(f) terms will exist, one with each + and - in the exponent, two of them cancel but two of them do not, and their exponents add up to having 4π terms.

    (a⋆ Φ)(a^*⋆ Φ^*)  = ^2 P^2/16[(2 + e^+i4π/p(vt -Δ) + e^-i4π/p(vt -Δ))δ(f) 
    
        + δ(f + 2/p) + δ(f - 2/p)] + P^2/4[δ(f+1/p) + δ(f-1/p)]

These two exponents can be combined with the identity that e^iαx + e^-iαx = 2 cos(αx)

    _2,halo = ^2 P^2/16[2(1 + cos(4π/p(vt -Δ)))δ(f)
        + δ(f + 2/p) + δ(f - 2/p)] 
     
        + P^2/4[δ(f+1/p) + δ(f-1/p)]


To find PSF_2,strehl, the second second-order term, we must perform an additional convolution, this time, with twelve terms.

    a^*^*^* = ∫_-∞^∞  (a^*^*)(f') Φ^*(f-f') f̣'


    = ∫_-∞^∞{ P/4[e^i2π/p(vt -Δ) δ(f') + e^-i2π/p(vt -Δ)δ(f') - e^i2π/p(vt +Δ)δ(f'-2/p)  - e^-i2π/p(vt +Δ)δ(f' + 2/p) ] 
        + P/2i[e^-i2π/pΔδ(f'+1/p) - e^i2π/pΔδ(f'-1/p)] }×P/-2i[e^i2π/pΔδ((f-f')-1/p) - e^-i2π/pΔδ((f-f')+1/p)] f̣'


    ∫_-∞^∞ P/4e^i2π/p(vt -Δ)δ(f')×P/-2ie^i2π/pΔδ((f-f')-1/p) f̣'    = - P^2/8ie^i2π/pvtδ(f-1/p) 
    ∫_-∞^∞ P/4e^i2π/p(vt -Δ)δ(f')×P/2i e^-i2π/pΔδ((f-f')+1/p) f̣'    =  P^2/8ie^i2π/p(vt -2Δ)δ(f + 1/p) 
    
        ∫_-∞^∞ P/4e^-i2π/p(vt -Δ)δ(f')×P/-2ie^i2π/pΔδ((f-f')-1/p) f̣'    = - P^2/8ie^-i2π/p(vt -2Δ)δ(f-1/p) 
    ∫_-∞^∞ P/4e^-i2π/p(vt -Δ)δ(f')×P/2i e^-i2π/pΔδ((f-f')+1/p) f̣'    =  P^2/8ie^-i2π/pvtδ(f+1/p)
    
        ∫_-∞^∞- P/4e^i2π/p(vt +Δ)δ(f'-2/p)×P/-2ie^i2π/pΔδ((f-f')-1/p) f̣'    =  P^2/8ie^i2π/p(vt +2Δ)δ(f'-3/p) 
    ∫_-∞^∞- P/4e^i2π/p(vt +Δ)δ(f'-2/p)×P/2i e^-i2π/pΔδ((f-f')+1/p) f̣'    = - P^2/8ie^i2π/pvtδ(f'- 1/p)
    
        ∫_-∞^∞- P/4e^-i2π/p(vt +Δ)δ(f' + 2/p)×P/-2ie^i2π/pΔδ((f-f')-1/p) f̣'    =  P^2/8ie^-i2π/pvtδ(f + 1/p) 
    ∫_-∞^∞- P/4e^-i2π/p(vt +Δ)δ(f' + 2/p)×P/2i e^-i2π/pΔδ((f-f')+1/p) f̣'    = - P^2/8ie^-i2π/p(vt +2Δ)δ(f + 3/p) 
    
        ∫_-∞^∞P/2ie^-i2π/pΔδ(f'+1/p)×P/-2ie^i2π/pΔδ((f-f')-1/p) f̣'    = P^2/-4i^2δ(f) 
    ∫_-∞^∞P/2ie^-i2π/pΔδ(f'+1/p)×P/2i e^-i2π/pΔδ((f-f')+1/p) f̣'    = P^2/4i^2e^-i4π/pΔδ(f+2/p) 
    
        ∫_-∞^∞P/-2ie^i2π/pΔδ(f'-1/p)×P/-2ie^i2π/pΔδ((f-f')-1/p) f̣'    = P^2/4i^2e^i4π/pΔδ(f-2/p) 
    ∫_-∞^∞P/-2ie^i2π/pΔδ(f'-1/p)×P/2i e^-i2π/pΔδ((f-f')+1/p) f̣'    = P^2/-4i^2δ(f)

Multiplying by a, only the terms with δ(f), δ(f-1/p), and δ(f+1/p) remain

    a(a^*^*^*)  = [ δ(f) + /2iδ(f-1/p)e^-i2π/pvt - /2iδ(f+1/p)e^i2π/pvt] ×{ P^2/8i[ ... ] + P^2/4[ ... ]}
    
        = P^2/2δ(f) + ^2P^2/16 i^2[( -e^i2π/pvt-e^-i2π/p(vt -2Δ)-e^i2π/pvt)e^-i2π/pvtδ(f-1/p)
    
         - (e^i2π/p(vt -2Δ) + e^-i2π/pvt + e^-i2π/pvt)e^i2π/pvtδ(f+1/p)] 
    
         = P^2/2δ(f) + ^2P^2/16[( 2 + e^-i4π/p(vt -Δ))δ(f-1/p) + (2 + e^i4π/p(vt -Δ))δ(f+1/p)]

where the ... represent the terms in (45) – (57) after they have been properly factored.
It is a good exercise to perform the identical calculations with every single term complex conjugated, to find the corresponding second half of the second second-order term a^*(a), which you may choose to do so on your own if you wish. Here we'll simply use the shortcut that since the δ function is strictly real, a^*(a) = [a(a^*^*^*)]^*, to much more rapidly arrive at the answer.

    a^*(a) = P^2/2δ(f) + ^2P^2/16[( 2 + e^i4π/p(vt -Δ))δ(f-1/p) + (2 + e^-i4π/p(vt -Δ))δ(f+1/p)]

Adding these two together, and using e^iαx + e^-iαx = 2 cos(αx),

    _2,strehl   = -1/2[a(a^*^*^*) + a^*(a)] 
       = -1/2{ P^2δ(f) + ^2P^2/16( 4 + 2cos(4π/p(vt -Δ)) )[δ(f-1/p) + δ(f+1/p)]}

In summary, we have calculated the following terms in the second order expansion of the AO-corrected and scintillated PSF.

    PSF_0 = δ(f) + 𝒜^2/4[δ(f-1/p)+δ(f+1/p)]    

    PSF_1 =  -𝒜P/2sin(2π/p(vt-Δ))[δ(f-1/p)-δ(f+1/p) ]    

    PSF_2,halo = P^2/4[δ(f+1/p) + δ(f-1/p)]  + ^2 P^2/16[2(1 + cos(4π/p(vt -Δ)))δ(f)
    + δ(f + 2/p) + δ(f - 2/p)]     

    PSF_2,strehl = -P^2/2δ(f) -𝒜^2P^2/32(4 + 2cos(4π/p(vt-Δ)))[δ(f-1/p)+δ(f+1/p)].   

In general, obtaining an expression for the n^th-order term in the expansion will require performing an n-fold convolution, so calculations beyond the second order become increasingly tedious to perform.



§ ADDITIONAL FIGURES





