Characterization of lemniscate atmospheric aberrations in Gemini Planet Imager data
A semi analytic framework for simulating the effects of atmospheric seeing in Adaptive Optics systems on an 8-m telescope is developed with the intention of understanding the origin of the wind-butterfly, a characteristic two-lobed halo in the PSF of AO imaging. Simulations show that errors in the compensated phase on the aperture due to servo-lag have preferential direction orthogonal to the direction of wind propagation which, when Fourier Transformed into the image plane, appear with their characteristic lemniscate shape along the wind direction. We develop a metric to quantify the effect of this aberration with the fractional standard deviation in an annulus centered around the PSF, and use telescope pointing to correlate this effect with data from an atmospheric models, the NOAA GFS. Our results show that the jet stream at altitudes of 100-200 hPa (equivalently 10-15 km above sea level) is highly correlated (13.2) with the strong butterfly, while the ground wind and other layers are more or less uncorrelated.
Characterization of lemniscate atmospheric aberrations in Gemini Planet Imager data
The Gemini Planet Imager (GPI) is an instrument installed on the Gemini South Telescope in Cerro Pachon, Chile, designed to search for thermal emission from young hot extrasolar planets at wide angular separation[Macintosh2014]. The GPI Exoplanet Survey (GPIES) is well under way, with many successful detections[Macintosh2015], but contrast remains limited by residual atmospheric aberrations[Ruffio2017][Poyneer2016] and imperfections in the Adaptive Optics (AO) system[MalesGuyon17]. Notably, the presence of the wind butterfly (so-called because of its figure-8 or lemniscate shape) scatters significant light into the coronagraphic dark hole and causes the residual point-spread function (PSF) to break azimuthal symmetry in the image plane, which reduces the final contrast ratio. This PSF pattern is consistent with wavefront errors from servo lag[Rigaut1998]. A ms delay in positioning of the deformable mirror in response to the wavefront sensor causes a displacement between the actual phase errors on the aperture from atmospheric turbulence when compared to the phase of the applied correction. The Fourier transform of this particular pattern has preferential direction, and appears in the image plane as the wind butterfly.
In this paper, we will demonstrate the servo-lag error as the origin of the wind butterfly in simulations of AO systems based on a theoretical turbulence model that is well-verified with a Kolmogorov power spectrum, as well as presenting our results in analyzing the appearances of the wind butterfly in the GPI Exoplanet Survey data. We have developed a metric to identify image subsets where this effect is highly apparent. The fractional standard deviation in an annulus centered about the PSF is zero in the case in azimuthal symmetry, and approaches one for cases of extreme azimuthal asymmetry. Using a chi-squared minimization routine with a model of the image’s azimuthal dependence, we can extract the wind direction in the image plane, which, along with telescope pointing information, identifies the wind vector in three dimensions. This allows us explore correlations between the wind vector as pointed by the butterfly, and the direction of high-altitude winds from the NOAA GFS database. Evaluating the azimuthal variations in planet detectability caused by the butterfly will help us to evaluate survey sensitivity and potential improvements from faster AO correction.
2 Butterflies in Simulations
In order to motivate understanding the appearance of the wind-butterfly in the GPIES, we would like to first develop a semi-analytic framework that describes turbulence in the atmosphere and how adaptive optics respond to provide a clear picture of the mechanism that brings the wind butterfly into existence. We begin with a short review on standard models of turbulence and their effect on optical imaging systems.
2.1 Theoretical Turbulence Models and Adaptive Optics Simulations
Tartarski[Tatarski1961] has shown that the fluctuations in the optical index of refraction in three dimensions for a Kolmogorov turbulence spectrum follow the form
Where is the index of refraction structure constant and and is the spatial wavevector for an eddy of size . 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[Johansson1994] with
where 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 is the unnormalized inverse Discrete Fourier Transform (DFT) given by
for a discrete array of size with . The discrete indices and 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 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 times the optical path length.
Here is the coordinate system in the aperture at , is the range of altitudes relevant to the turbulent layer at altitude , and the baseline index of refraction of the atmosphere can be approximated[Hardy1998] with
where and are the pressure (in millibars or equivalently hPa) and temperature (in Kelvin) of the atmosphere for a particular altitude. We will obtain a model of the index of refraction for the baseline atmosphere in a later section, but once that is acquired we simply add the fluctuations on top to simulate a particular turbulent instance. For non-zenith observations an additional term of where is the zenith angle should be included in the integral in (4) 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.[MalesGuyon17]
Furthermore, we assume 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
where is the wind velocity at altitude , which is assumed to lie only in the plane at altitude with no vertical component, 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 whose only flaw is a delayed response to test our hypothesis for the origin of the wind butterfly. In essence, this is an ideal open-loop AO simulation. The expression for the compensated phase in the aperture is then
where we sum over the contributions from turbulent layers at altitudes with a flat step interpolation scheme for the structure constant. For a continuous wind velocity profile, one could transform this summation into an appropriate integral, but we will be working with discretely layered wind data. From the compensated phase on the aperture, we can obtain the final image’s intensity distribution with a Fourier transform by assuming the telescope focus operates in a Fraunhofer diffraction limit, so that electric field distribution in the image plane is the Fourier transform of the aperture function[Hecht2002].
Here are the coordinates in the image plane, inside the aperture and zero outside and the brackets denote time average. We consider a telescope of aperature diameter m and we do not consider apodization at the moment.
2.2 The NOAA GFS as an atmospheric model for Cerro Pachon
The National Oceanic and Atmospheric Organization’s (NOAA) Global Forecast System (GFS) is a weather model developed by the National Centers for Environmental Prediction (NCEP), and for the considerations of the paper will be considered the atmospheric truth. The dataset is distributed as gridded data over the entire globe, will half degree scale spatial resolution, and a temporal resolution of every six hours. The dataset contains a significant number of atmospheric variables, although we will be primarily interested in the U (East) and V (North) components of wind velocities at various altitudes for construction of the profile, as well as the temperature as a function of pressure altitude to estimate the baseline index of refraction.
Because pressure variations are equalized at the speed of sound, it is common for an atmospheric model to assume that various altitudes are isobaric, to the point where pressure is the actual coordinate value in the z-direction. It is possible to formulate a rough model which relates the pressure and altitude by assuming the Earth is a uniform sphere with uniform surface temperature with an atmosphere of average molecular mass . The pressure at altitude is then
where is the atmospheric pressure at sea level, is the local surface gravity, and is the ideal gas constant. Using this conversion between atmospheric pressure coordinates to physical altitude, we can use the data from the NOAA GFS to estimate the index of refraction structure constant . Applying a Hufnagel turbulence model[Hardy1998], we use the following equation to convert wind velocities into turbulence strength.
Considering the Gemini South Observatory at Cerro Pachon’s Coordinates of latitude -30:14:26.700 and longitude -70:44:12.096, we select the gridded bin in the GFS centered on latitude 30 South and 70.5 West, from December 7th 2015, to May 17th 2018, dates for which are relevant to the GPI Exoplanet Survey, and calculate the turbulence profile from the wind velocities over the range of altitudes available. The resulting turbulence profile, along with other atmospheric variables is plotted in Figure 1.
As a method of verification for this turbulence profile, we calculate the value of the Fried parameter at m with[Hardy1998]
and find that it is roughly 12 cm, which is better than average. Seeing conditions are arcsec, which are fair, but slightly worse than typical conditions at Cerro Pachon, which typically sees cm. This is due to our estimate integrating to altitude z=0 at sea level, rather than the altitude of the observatory.
2.3 Simulation Results
|U wind (m/s)||-0.38||-0.37||-0.37||-0.38||-0.38||-0.37||-0.32||2.07||5.55||10.3||17.4||24.3||30.1|
|V wind (m/s)||-3.47||-3.47||-3.47||-3.47||-3.47||-3.47||-3.54||-4.20||-6.65||-8.34||-8.01||-5.86||-3.91|
|U wind (m/s)||35.2||40.1||46.7||53.8||56.3||59.9||52.5||18.8||11.8||13.5||3.60||2.70||18.0|
|V wind (m/s)||-2.50||0.50||8.95||16.6||22.4||25.6||15.4||3.19||0.76||6.61||1.50||0.60||-1.20|
By combining the profile we have found above with a particular set of velocities at the various altitudes encapsulated in the GFS model, we should have enough information to perform a rough model of the time evolution of the atmosphere and the response of an AO telescope to seeing such an atmosphere.
To start, L=26 unique Kolmogorov screens are generated, each with a unique realization of the random noise process , and one for each layer of the atmospheric model described earlier. Then the timestep is incremented in units of ms, which is the time delay for an open-loop control system which roughly corresponds to the closed-loop response of GPI, which has a 3 dB bandwidth of around 20 Hz[Poyneer2014]. This time delay also gives nm RMS residual WFS error, at an atmospheric coherence time of ms, also roughly corresponding to models of GPI AO[Poyneer16]. (As an aside, this atmospheric coherence time is unusually short, and bring into question the frozen-flow hypothesis. A more robust simulation would include boiling in the atmosphere with a clever autoregressive technique[Srinath2014], but for now we will continue to assume the frozen-flow hypothesis.) Each kolmogorov screen is translated according to the U and V components of the velocity given by the GFS model instance in Table 1, with the appropriate conversion between m/s and pixels, and with overflow causing the screen to re-enter on the opposite side, also referred to as periodic boundary conditions. At each step, the phase errors on the aperture are calculated from integrating the index of refraction and the fluctuations over the atmosphere. The compensated phase is calculated by subtracting the phase in the current timestep from that of the previous, to simulate a servo-lag error in an AO system. Then the compensated phase is Fourier transformed with a multiplicative factor for the telescope pupil to result in the final image plane. This process is demonstrated visually in Figure 2. One can see the heuristic behavior of the atmosphere and AO systems which generate the butterfly shaped halo in the PSF, and it indeed points along the direction of the wind.
3 Butterflies in the GPIES
Having established a semi-analytic framework and model for the origin of the wind butterfly, we will now examine its presence in the observational data taken during the GPIES campaign. Although the effect is visibly apparent and easy to spot by eye, the sheer number of images taken renders hand separation a futile task. It is with the intention of being able to distinguish images that contain wind butterflies from those who do not that we develop a metric to distinguish for us.
3.1 Evaluation Metrics
For either simulated images, or real telescope data, there are two primary metrics of interest for evaluating the presence and strength of the wind butterfly. The first is the image angular profile, defined as
Which is simply a radial integral of the image over an annulus between and , to ignore complexities with coronagraphic scattering in the dark hole, and simply to examine the azimuthal symmetry or asymmetry of the halo. For azimuthal symmetry, is constant, but in the presence of wind, in general, an image from adaptive optics will not be azimuthally symmetric. To first order, we can approximate the angular profile with a sinusoid of frequency to fit for the preferred butterfly and wind direction.
Such a model is a function of angle in the image as well as a parameter vector , whose components are amplitude, phase, and a constant offset.
In practice, we will only be interested in , as it tells us of the preferred butterfly direction in the image, which is what we are interested in correlating with wind data, the other two degrees of freedom in this model simply allow the fitting to properly converge. From this model, we can construct a total error
which, once integrated over all angles, is a function simply of the parameter vector. We ignore normalization on at the moment since we only care about the minimum value. To obtain the minimum value of such a function, we use scipy.optimize.fmin along with a null guess for the parameters, and it converges on the optimal value for our parameter vector. An example of this fitting is demonstrated in Figure 3.
Of interest is also the fractional standard deviation in the annulus , defined as
where is the mean value over the annulus
and is simply the number of discrete points or pixels under consideration for normalization. The fractional standard deviation is theoretically bounded on , but for any reasonable and continuous image is bounded by , where we consider 0 to be the case of azimuthal symmetry, and 1 to be the extremely asymmetric case of a square wave radial profile. A sine wave squared that varies from a maximum amplitude of one to zero has , to give a typical estimate.
3.2 Telescope Pointing and 3D orientation of GPI Images
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
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 4 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
So that one can think of as pointing in the direction of increasing Azimuth, and pointing towards increasing Elevation. It is left to the reader to show that , and that 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
However, we would instead like to know . By algebraically substituting in our coordinate vectors , and formulas in x,y,z space, we can arrive at an expression for the wind vector in x,y,z space in terms of , , and . This is
With this done, we can easily project the vector onto the ground plane by simply removing the z-component. If we need to find the direction of this wind vector as an azimuth, we can use the following trick
Where , 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
where is the hour angle, is the local sidereal time in radians, is the local latitude, is right ascension, is declination, and the addition of is due to the strange convention that azimuth starts from the South and opens to the West, but here we use the convention that it 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 points towards increasing elevation. Then, we use the techniques described earlier to find preferential butterfly wind vectors using our chi-squared minimization, and project them onto the ground.
3.3 Observational Results
We investigate a selection of 22244 images taken during the course of GPIES campaign, and using our previously described techniques, calculate the fractional standard deviation in the images as well as determine the butterfly’s wind vector and project it onto the ground coordinates. Figure 5 shows a histogram of all of the fractional standard deviations and their total number of occurrences over 40 bins. The average is around .32 and the standard deviation is around .08, and so we select a subset of very strong butterflies with greater than above average. This subset constitutes 3178 images.
The resulting dataset is matched temporally to the NOAA GFS, and correlations between the directions of the strong butterfly and the various wind layers are shown in Figure 6. Due to the availability of atmospheric data, only 1105 images are able to be compared. For each wind layer, a simple linear model is fit to the relationship between the butterfly direction and the wind layers direction, which both exist numerically as degrees azimuth in the coordinates described earlier. Slopes and their respective uncertainties are estimated to give a sense of the strength of the correlation, which should be around zero for uncorrelated, and approach one for strong correlation. The fundamental observation of this paper is the strong correlation of the butterfly and the wind at the altitudes around 100-200 hPa, or around 12-16 km, colloquially known as the jet stream. The rest of the layers exhibit either mild correlations, or strangely near the ground negative correlations, but this may be attributable to correlations in the atmospheric model from the continuity of fluid flow.
To digest the large amount of scatter plot and fit lines into a more concise form, Figure 7 was created. Both the slopes of the fit lines and the Pearson’s R-coefficient defined as
for each scatter plot was graphed as a function of altitude. One can see both that the value of the slope and the strength of the correlation for the scatter diagrams peaks in the appropriate range of the jet stream, although the strength of the correlation, given the by the Pearson’s R-coefficient, seems to peak at a lower altitude than the correspondence of the correlation, given by the slope of the best-fit line. The reason for this is unclear, but does not significantly affect our results at this time.
In order to quantify the likelihood of these correlations spontaneously forming due to chance, a null hypothesis of randomly distributed butterfly directions was compared to a randomly selected sample from the prior wind distributions. Our sample size of 1105 points over one hundred thousand iterations were generated in order to bootstrap an estimate for possible observed values of R or slope. The bootstrapped uncertainties in the sloped varied heavily over altitude, and the resulting areas are shown as shaded on Figure 7. The Pearson-R formed a normal distribution with a mean of zero and standard deviation of regardless of altitude, which puts the significance of our correlations at for , or equivalently a 1 in probability of happening due to chance.
Although it is now clear that there is a strong relationship between the directionality of the wind butterfly and the strength of the jet stream, it is worthwhile to note another particular effect that was discovered in both the simulations and data regarding the magnitude of the the fractional standard deviation. This effect is to note that there exists a maximum value of and thus the apparent strength of the butterfly asymmetry corresponding to a particular value of the wind speed. Figure 8 demonstrate this visually.
Although there is not excellent agreement between the simulation and and the observations, one can see a similar effect in both the simulated and observed with wind velocities between 35-40 m/s in that the fractional standard deviation reaches a maximum. This makes sense when considering the analytic limits of either taking the wind velocity to infinity, or equivalently, making the AO delay infinitely slow. In this limit, the observed phase on the aperture and the correction applied by the deformable mirror will be two entirely uncorrelated Kolmogorov screens, whose difference would similarly be Kolmogorov. This would cause the resulting image to have azimuthal symmetry, as it would just be like observing through an uncorrected atmosphere. The image produced would be atmospheric speckles resolved at forming an approximately Gaussian blob with FWHM . Between the peak velocity and this limit, the simulations show that the butterfly transforms into a roughly elliptically shaped halo, which should have intermediate values of . In the opposite limit, taking the wind speed to zero, or equivalently taking the AO delay to zero, the correction would be perfect, with zero phase on the aperture resulting in a standard diffraction-limited Airy ring, which is also azimuthally symmetric. The very existence of the butterfly in the intermediate regime necessitates there being a particular combination of the wind velocity and AO delay corresponding to maximum azimuthal asymmetry. For our simulations, we find that this constant centimeters, which corresponds to the physical distance between observed and corrected phase aberrations which would produce the greatest azimuthal asymmetry in the images. It is not known if this constant is a function of telescope diameter or observing wavelength.
The implications of this effect in terms of AO improvements is tenuous. Planet detectability algorithms often rely on the assumption that PSF noise is azimuthally symmetric, which is clearly not the case in all images, and so appropriate modeling of the asymmetry of the butterfly halo could provide improvements to detections. Improvements to AO systems response times, that is, decreasing the delay constant would certainly provide improvements as well, although any finite delay will still result in this image asymmetry, unless clever predictive models are developed which can anticipate atmospheric changes to apply corrections psuedo-instantaneously.
To conclude with a summary, we have demonstrated a semi-analytic mechanism to describe the origin on lemniscate atmospheric aberrations or colloquially, the wind butterfly, in adaptive optic telescopes, using simulations of propagating Kolmogorov turbulence and a delay to account for servo-lag errors in a real system. This demonstration motivated an exploratory data mine into images from the GPIES campaign, as well as atmospheric models. By appropriately manipulating the data, we can correlate the projected butterfly direction onto the ground with the direction of various winds layers from the atmospheric models, and we find that the wind butterfly is strongly correlated () with the high altitude jet-stream layers of wind around 10-15 km above the surface of the Earth. This claim reaffirms our understanding of the models that govern turbulence and adaptive optics systems, and provides an underpinning to understanding how planet detectability algorithms should handle azimuthal asymmetry in images taken during the campaign, as well as highlighting the prospective improvements that could arise from faster corrective algorithms.