GENERATION OF REALISTIC INPUT PARAMETERS FOR SIMULATING ATMOSPHERIC POINT-SPREAD FUNCTIONS AT ASTRONOMICAL OBSERVATORIES

,


INTRODUCTION
Images of distant galaxies are distorted by weak gravitational lensing due to inhomogeneities at large scales in the mass distribution in the Universe between the galaxy and the observer.These distortions are called cosmic shear (see, e.g., Kilbinger 2015;Mandelbaum 2018).Spatial correlations in the distortions are sensitive to the properties and evolution of the matter density on large scales and the geometry of space.As the Dark Energy Survey (DES 1 , Amon et al. 2022;Secco et al. 2022), Hyper Suprime-Cam (HSC 2 , Li et al. 2023; Dalal 1 https://www.darkenergysurvey.org/ 2 https://hsc.mtk.nao.ac.jp/ssp/ et al. 2023), and Kilo-Degree Survey (KiDS3 , Asgari et al. 2021) collaborations complete the analysis of their entire surveys, cosmic shear is becoming one of the most precise probes of cosmology -in particular for measurements of the average matter density and the amplitude of matter fluctuations.When the commissioning of the Vera C. Rubin Observatory4 is complete, the Legacy Survey of Space and Time (LSST) is expected to reach and then surpass the precision of existing surveys for cosmic shear by an order of magnitude.
For accurate (unbiased) measurements of cosmic shear, all other sources of correlated image distortions must be modeled and/or removed.One such source is the correlated blurring of images due to Kolmogorov turbulence in the atmosphere, which dominates the point-spread function (PSF) for ground-based instruments such as DES, HSC, KiDS, and Rubin Observatory (e.g., Heymans et al. 2012;Jarvis et al. 2016;Xin et al. 2018).Algorithms are under active development for more accurately modeling and interpolating the PSF across the focal plane (see, e.g., Bertin 2011;Jarvis et al. 2020), correcting galaxy shapes for the PSF, and calibrating measures of cosmic shear from galaxy shapes (Sheldon & Huff 2017;Huff & Mandelbaum 2017;Gatti et al. 2021;Sheldon et al. 2023).Both real astronomical images and high-fidelity simulated images are important tools for developing and optimizing these algorithms and measuring their performance, with simulations playing a unique role because the input parameters we are trying to measure or infer, including shear, are known.High-fidelity atmospheric PSF simulations have played an important role in astronomy -from the development of adaptive optics (Jolissaint 2010;Srinath et al. 2015;Madurowicz et al. 2018) to the optimization of instrumentation and software for precision cosmology, including Rubin Observatory and the associated survey, LSST (Jee & Tyson 2011;Chang et al. 2012;Peterson et al. 2015; The LSST Dark Energy Science Collaboration et al. 2021).Simulations of the atmosphere often use a "thin-screen, frozen-flow" approximation in which the impact of a layer of atmosphere is modeled as a single two-dimensional planar "screen" in which the relative wavefront phase across the plane encodes the impact of variations in index of refraction in the layer of atmosphere (turbulence).The variation in phase across the screen is "frozen" during an exposure, but the screen moves across the field of view to simulate the impact of wind.The altitude dependence of wind speed and direction and of turbulence strength is modeled by including multiple screens, as illustrated in Figure 1.
In preparation for analysis of data from Rubin Observatory, the LSST Dark Energy Science Collaboration (DESC) produced an extensive set of image simulations, called Data Challenge 2 (DC2, The LSST Dark Energy Science Collaboration et al. 2021, hereafter DESC DC2 Sims), in which turbulence parameters were based on median conditions determined in Ellerbroek (2002) from measurements at Cerro Pachón (the location of Rubin Observatory), and wind speeds and directions were drawn randomly from uniform distributions.In Peterson et al. (2015), historical data are used for wind speeds and directions, and turbulence strengths; however, correlations between meteorological parameters at different altitudes, and the relationship between wind and turbulence parameters, are not taken into account.
As described above an important scientific goal in simulating atmospheric PSFs at Rubin Observatory or other wide-field instruments is to predict how weather in the observatory environment impacts correlations in PSF parameters across the focal plane.Therefore, it is impor-Fig.1.-This schematic illustrates the simplified view of the atmosphere used for PSF simulations based on discrete phase screens.Lines of sight (dot-dashed black) for two stars (whose images are located at extrema of the field of view) pass through two phase screens of von Kármán refractive index variations, each with different values for the outer scale L 0 .The phase offset incurred by light passing through each point on each screen is indicated by the color scale in units of wavelength.The columns (solid black) associated with each line of sight show the path of starlight that will reach the telescope aperture (black), along with the relevant phase screen area.The wind vectors (blue arrows) show speed and direction of the wind in the plane of the screen.The primed coordinate systems are perpendicular to the telescope axis, and are related to the ground coordinate system via the altitude and azimuth of the pointing.
tant to include in the simulations altitude-dependent correlations among wind speed, wind direction, and turbulence, as well as realistic temporal variations.We have produced a public software package (psf-weatherstation) that leverages local environmental telemetry and data products from global weather forecasting models to produce turbulence parameters and wind speeds and directions that are realistically correlated with altitude and in time.The package relies on an empirical model of atmospheric turbulence proposed in Osborn et al. (2018) (hereafter OS18), which parameterizes relative turbulence strength at different altitudes in the atmosphere as a function of the wind shear at that location.
Whereas OS18 focuses mainly on fast predictions of turbulence as a function of altitude for real-time adaptive optics corrections in very large telescopes with narrow fields of view, our goal is to parameterize atmospheric conditions that, in an ensemble sense, are representative of a particular site and can be used to generate atmospheric PSFs across a wide field of view.We achieve this through the inclusion of weather-tower telemetry and site-specific empirical distributions of altitude-specific seeing contributions.We then use psfweather-station to predict and study the expected anisotropies in the PSF, the effects of different weather patterns, etc., at the location of a specific observatory.
In this paper, we present an application of the psf-weather-station package to studies of the expected PSF at Cerro Pachón (Rubin Observatory).However, the code is flexible enough to use for other observatories and includes functionality to download necessary datasets from weather forecasting services.We describe how the atmosphere and PSFs are modeled in Section 2, outline the psf-weather-station package in Section 3, describe three sets of inputs to simulations in Section 4, define PSF parameters and twopoint statistics in Section 5, and compare these PSF metrics for the three types of simulation inputs in Section 6and Section 7. We end with a discussion of implications for cosmic shear analyses and future work.

IMAGING THROUGH A TURBULENT ATMOSPHERE
Since stars and galaxies are effectively light sources at infinity, their light can be treated as plane waves when entering the upper atmosphere.During the journey through the atmosphere, points on the surface of a wavefront accrue relative phase shifts.The atmospheric component of an object's PSF is the result of the spatial variations in phase across the telescope pupil.
The phase shifts are caused by variations in the index of refraction in the atmosphere (δn) due to perturbations in air density driven by turbulent mixing of air at different temperatures (Lawrence & Strohbehn 1970;Clifford 1978).These fluctuations in n (referred to as optical turbulence) vary in space and time; therefore, each photon in general incurs a slightly different cumulative phase shift on its path to the telescope pupil.It is convenient to define an "atmospheric column", with diameter roughly that of the telescope pupil, which delineates the volume of turbulent air sampled by the imaged photons from a single source, as illustrated in Figure 1.The atmospheric columns for each object in the field of view overlap at the pupil but diverge with distance from the telescope, resulting in a spatially varying, spatially correlated PSF over the focal plane.
Optical turbulence exists for a range of spatial scales and amplitudes.The spectral density of this turbulence, as a function of spatial frequency κ, can be described by the 3-dimensional von Kármán power spectrum (von Kármán 1948;Tokovinin et al. 1998), where the subscript n denotes index of refraction: This is a modification of Kolmogorov's κ −11/3 power law (Kolmogorov 1941); the (altitude dependent) outer scale parameter L 0 sets an upper bound on the turbulence strength at low spatial frequencies, which would otherwise diverge as κ → 0. The von Kármán turbulence spectrum has an associated spatial correlation function (related to the Fourier transform of Equation 1).The degree of correlation for optical turbulence at a given altitude is set by the turbulence structure constant C 2 n ; when given as a function of altitude h, this is known as the optical turbulence profile (OTP) C 2 n (h).Although turbulent structures are constantly evolving with time, OTPs are typically assumed to be constant during the course of an exposure (Roddier 1981).
Image quality is related to the OTP via the turbulence integral J: For the case when the integration bounds h 1 and h 2 correspond to the entire vertical extent of the atmosphere, J quantifies the total strength of the turbulence experienced by photons passing through the corresponding atmospheric column.C 2 n and J have units of m −2/3 and m 1/3 , respectively.The Fried parameter r 0 is a characteristic length scale that defines the radius of an aperture in which the wavefront phase variance is approximately 1 rad 2 (Fried 1965).It depends on J as well as zenith angle ζ and wavenumber k (Roddier 1981): and is inversely proportional to linear PSF size (FWHM).
In low turbulence conditions, r 0 is large (and J is small), and thus the PSF size is small.The atmosphere is not a static system; in addition to turbulent mixing, there is large-scale motion driven by wind.The component of wind velocity parallel to the telescope pupil translates the optical turbulence across one or more atmospheric columns, leading to correlated phase shifts for photons with different positions and angles at the pupil plane, and therefore correlated PSF shapes across the focal plane.

Atmospheric PSF simulations
Since it is computationally intractable to simulate atmospheric PSFs by calculating the trajectory of each photon through a turbulent (i.e., chaotic) 3D volume of atmosphere, approximate methods have been developed.In this section, we describe the frozen screen approximation, introduced in Section 1, which has been used in the context of weak lensing studies; see Jee & Tyson (2011); Peterson et al. (2015); DESC DC2 Sims for more details.
Measurements of optical turbulence with SCIntillation Detection And Ranging (SCIDAR) instruments show that the atmosphere is often stratified into regions of stronger turbulence separated in altitude by areas of relative calm (Osborn et al. 2018;Osborn & Sarazin 2018).Typically only ∼ 1 km in vertical extent and variable in number, these layers of stronger turbulence dominate the atmospheric contribution to the PSF.These observations motivate a simplified model of the atmosphere that consists of only 2-dimensional phase screens across which the refractive index varies, with each screen representing a layer of turbulence.
The refractive index variations within each phase screen are a realization of von Kármán turbulence.We assume Taylor's frozen flow hypothesis (Taylor 1938), in which the time scales for changes in turbulence are longer than those for changes due to phase screen drift from wind.Under this assumption, it is not necessary to evolve the turbulence structure during a simulated exposure.Instead, each phase screen is assigned a "wind" speed and direction; for each time step ∆t of the simulation, the phase screens are translated accordingly.A schematic of such an atmospheric simulation, with two phase screens, is depicted in Figure 1.After each time step, the phase screen contributions within the atmospheric column (for each star in the field) are summed vertically.These wavefront phase variations have a von Kármán power spectrum (cf.Equation 1 for E n (κ; r 0 )): The wavefront outer scale L 0 is the spatial scale at which correlations in wavefront phase saturate; it can be expressed as the turbulence-weighted sum of the turbulence outer scale L 0 over phase screens i: (Borgnino 1990;Tokovinin et al. 1998).For each star, the wavefront is then Fourier transformed to focal plane coordinates and, after all time steps, added together to form the simulated image -i.e., the PSF.As a function of image coordinates θ x , θ y , where P (u, v) is the aperture transfer function and W (u, v, ∆t) is the wavefront, with each a function of pupil coordinates u, v.The sum over ∆t represents the sum over all simulation time steps during an exposure.The phase screen is the building block of the simplified atmospheric model described above.Because turbulence integrals J i (Equation 2) add linearly, each phase screen contributes to the total turbulence with weight w i = J i ( i J i ) −1 .Although one could use the turbulence integral J i to generate the phase pattern for screen i, it is more natural to use the Fried parameter r 0 because r 0 determines the turbulence power spectrum amplitude.Given that r 0 ∝ J −3/5 (Equation 3), the contribution of the ith screen is r 0,i = w −3/5 i r 0 .By convention, r 0 is specified at λ = 500 nm and zenith angle ζ = 0.
In summary, the input parameters for a simulation of PSFs across a single exposure are the outer scale L 0 , the atmospheric seeing parameterized by r 0 , the number of phase screens and their altitudes, the wind speed and direction for each screen, and the fractional contribution w i of each screen to the total turbulence.

psf-weather-station
In the psf-weather-station software package, we leverage data products from weather forecasting organizations to produce realistically correlated wind and turbulence parameters.We use these input parameters to simulate correlated atmospheric PSFs across a large field of view, as described in later sections.

Motivation
A potential source of bias in analyses of cosmic shear is uncorrected, spatially correlated noise (Mandelbaum 2018).The atmosphere is correlated via the von Kármán power spectrum described in Equation 1and, as we have seen, these spatial correlations translate into angular correlations in the size and shape of the atmospheric PSF in the associated exposure (Heymans et al. 2012;Jarvis et al. 2020).Wind over the telescope plays an integral role in this process, as it moves correlated patches of turbulence through the atmospheric columns that impact the images of different objects, leading to correlations on angular scales larger than the patches.If wind directions are consistent across altitudes, turbulence at different altitudes will imprint a stronger correlation in the PSF than when wind directions at different altitudes are uncorrelated.
Another relevant factor for PSF correlations is the altitude dependence of the optical turbulence profile (OTP), which describes the contribution of each layer to the total turbulence strength (Roddier 1981).Interestingly, one of the drivers of atmospheric turbulence is wind itself -specifically, wind shear (Masciadri et al. 2017)-so we expect that these two factors that influence spatial correlations in PSF parameters are not independent.

Data inputs to psf-weather-station
In psf-weather-station, we separate the atmosphere into two regions based on the typical turbulence regime for those altitudes.The ground layer (GL) is typically defined as the region between ground level and 500 − 1000 m above the telescope, where complex topography and heat sources generate non-Kolmogorov eddies.The free atmosphere (FA) is defined as the region above the ground layer, where turbulence is generally well-described by Kolmogorov statistics.This separation into GL and FA plays an important role in many design choices for psf-weather-station.
The primary sources of data for psf-weatherstation are data products from global weather forecasting organizations such as the European Centre for Medium-Range Weather Forecasts 5 (ECMWF) and the National Oceanic and Atmospheric Administration National Centers for Environmental Prediction6 (NOAA NCEP).The global circulation models (GCM) used in weather forecasting cover the entire globe on a grid of 0.25-2.5 deg resolution and output predictions for dozens of environmental parameters at a number of different altitudes between 0 and ≈80 km above sea level.
Although any of the available GCM data sets can be used in psf-weather-station, the results in this paper are based on data from the ECMWF Reanalysis v5 (ERA5) catalog; these are hourly estimates of global weather conditions based on meteorological measurements assimilated into the GCM.ERA5 was chosen for its denser sampling both in time -hourly -which is useful for sampling conditions throughout the night, and in altitude -137 levels -which is important for capturing vertical wind gradients in the atmosphere.We provide more details on the ECMWF and NOAA NCEP models in the Appendix.
The GCMs give us robust estimates of wind and temperature throughout the free atmosphere as a function of time.However, interactions of the atmosphere with the ground are not accurately captured because topographical features are not modeled at scales smaller than ∼1 km; therefore, the accuracy of the GCM data (both initial conditions and predictions) is limited near the ground.We overcome this limitation in psf-weatherstation by using measurements from a weather tower on the telescope site, rather than GCM data, for the ground layer.Since weather tower data are typically recorded every few minutes, the sampling times of the telemetry information can be matched to the GCM data used for the free atmosphere.The weather tower measurements are optional inputs to psf-weather-station but highly recommended since the GL is of significant importance for the PSF, contributing between 40 and 60% of the turbulence at many observatories (Tokovinin et al. 2003;Tokovinin & Travouillon 2005;Tokovinin et al. 2005).
3.3.Optical turbulence in psf-weather-station Some existing non-hydrostatic atmospheric models with sub-kilometer horizontal resolution are successful in simulating optical turbulence around observatories (Masciadri et al. 2001(Masciadri et al. , 2017)); however, such models are computationally prohibitive when many realizations need to be simulated.On the other hand, useful parameterizations of optical turbulence as a function of environmental quantities -such as temperature, pressure, wind, and kinetic energy -can be adapted from this literature, as is done in OS18 (Osborn et al. 2018).In particular, in OS18, wind shear is assumed to be the sole contributor to the kinetic energy term (i.e., wind shear is the source of turbulent mixing of air) and the temperature and pressure profiles from mesoscale simulations are replaced with GCM data.In exchange for coarser resolution and more limited accuracy, with minimal computational time the OS18 empirical model produces estimates of C 2 n (h) which, as shown in OS18, are broadly consistent with stereo-SCIDAR measurements.
The OS18 model captures variations in turbulence strength with altitude, but not the absolute strength; it requires calibration for the total turbulence J.In addition, the OS18 model significantly under-predicts turbulence in the GL, which is expected since turbulence in the GL can have significant contributions from sources other than wind shear.psf-weather-station combines the OS18 optical turbulence model with complementary information from the literature to produce correlated turbulence parameters -a value of J for each phase screen, including the ground layer -as described below.
Measurements of the altitude dependence of atmospheric turbulence with multi aperture scintillation sensor and differential image motion monitor (MASS-DIMM) instruments at a variety of sites show that turbulence contributions from the FA and the GL are independent (Tokovinin & Travouillon 2005;Tokovinin et al. 2005;Chun et al. 2009).Motivated by this independence, the total turbulence in the GL and FA layers (J GL and J FA ) are treated separately in psf-weather-station.The relative amount of turbulence contributed by each FA layer is calculated with GCM data and OS18, and the total GL and FA integrals are drawn from log-normal distributions fit to published quantiles of measurements of J GL and J FA .
In Figure 2, we illustrate the steps taken to go from raw GCM data to simulation-ready input parameters in psfweather-station, with Cerro Pachón as an example site.Our data sources are ground layer telemetry from the weather tower at Gemini South (also located at Cerro Pachón, about 1 km from Rubin Observatory), ECMWF ERA5 data products,7 and J GL and J FA quantiles from Tokovinin & Travouillon (2005).The plots in the top two rows of the figure show six months of site telemetry (orange dots and histogram) and ERA5 data (purple curves and histogram).The GCM altitude profiles are sampled at 00h, 02h, and 04h local time; each profile has a corresponding co-temporal ground-layer data point.The frequency distributions on the right illustrate the range of wind speeds and the distinctly non-uniform distribution of wind directions, both near the ground and across most altitudes.In this paper and in psf-weather-station, the direction is defined as the angle from which the wind is blowing, East of North.
The heavy purple curve in each of the top two plots in Figure 2 corresponds to an example time that is representative of the median weather conditions.This representative example will serve to illustrate our process in this and the following sections.The OS18 output for this representative example is shown in the bottom panel of Figure 2; the light purple curve is the uncalibrated data and the dark curve is shifted to match a calibration value chosen randomly from the J FA distribution.
Given the number N of desired phase screens, psfweather-station places one in the ground layer, and the free atmosphere is divided into N − 1 equally spaced layers.Here, we choose to calculate turbulence integrals for N = 6 layers of the atmosphere.In the bottom panel of the figure, the purple vertical dashed lines correspond to the boundaries between the six layers; the orange dashed line corresponds to the altitude of ground, here 2737 m at Gemini-South.The dots correspond to the values of the turbulence integrals (J) divided by the width of the layer, and are placed at the C 2 n -weighted position within each FA layer (purple dots) and at 200 m above the elevation of the observatory (orange dot) to match the location of the ground-layer screen in DESC DC2 Sims and thereby facilitate comparisons described in Section 4. The variable number of layers and the turbulence-weighted altitudes of the corresponding phase screens together offer some ability to model the temporal variability in the turbulence profiles observed in SCIDAR data (OS18).The complete set of parameters returned by psf-weather-station is h, J(h), v(h), and ϕ(h), where v(h) and ϕ(h) are the interpolated wind speed and direction at height h of each phase screen.psfweather-station does not currently provide estimates for the outer scale L 0 due to the lack of physically motivated models linking environmental parameters to the outer scale.
Using random draws from turbulence distributions for ground turbulence integrals and calibration of FA C 2 n is not an optimal solution since these draws are not temporally correlated with other environmental parameters.There is currently some observational evidence from a variety of observatories (Tokovinin et al. 2003(Tokovinin et al. , 2005;;Chun et al. 2009) for correlation existing between ground wind speed and J GL , so while we include in psf-weatherstation an option to correlate the random J GL draws with ground wind speed, by default the correlation is set to zero.There is also possible correlation of J GL with ground wind direction (Tokovinin et al. 2003) at some observatories, but we have not yet implemented such an option in psf-weather-station.Since there is only limited empirical evidence of correlations between J FA and FA wind speeds (Tokovinin et al. 2003) -and in the OS18 model turbulence already depends on wind shear -we do not include an option to introduce correlations.This method of using random draws from empirical distributions does somewhat restrict the predictive capabilities of simulations run with psf-weather-station, as we do not expect to recover the average seeing on individual nights.(Predicting the seeing on individual nights would require access to either MASS-DIMM or SCIDAR measurements, at or near the relevant observatory, that could be temporally matched to weather forecasting data products.)However, we expect to recover overall seeing statistics as well as spatial correlations of the PSFs.

SIMULATIONS OF PSFS AT CERRO PACH ÓN
psf-weather-station uses multiple sources of telemetry and vetted models to generate sets of correlated parameters for input to simulations of PSFs across the field of view.In this section, we describe tests of these generated parameters, which aim to quantify how simulations that use as input psf-weather-station parameters compare to earlier generations of atmospheric PSF simulations with uncorrelated parameters.
All simulations described here are generated with the GalSim8 software library (Rowe et al. 2015).We use the same GalSim implementation as described in DESC DC2 Sims, which follows ray-tracing methods developed in Jee & Tyson (2011) and Peterson et al. (2015). 9e generate atmospheric PSF simulations for Rubin Observatory with three types of input parameters.
1. psfws: In the first case, input parameters are generated for Cerro Pachón using psf-weatherstation with the data summarized in Figure 2. Six phase screens are used; the altitudes are allowed to vary according to the C 2 n scheme described in Section 3.3.

bench:
As a second case, we use as a benchmark the input atmospheric parameters used in the DESC Data Challenge 2 (DC2) image simulations (see DESC DC2 Sims).For each of six phase screens, the wind speed and direction are drawn from uniform distributions between 0-20 m/s and 0-360 • , respectively.Small (∼10%) Gaussian variations around the turbulence integrals from Ellerbroek ( 2002) are introduced, but the associated six altitudes remain fixed between simulations.

match:
As a third case we use the same values of input parameters as in bench -except for the wind directions, which are matched to the correlated wind directions used in psfws.The motivation for this match simulation is to identify whether differences between distributions of PSF parameters for the first two cases are mainly driven by the highly correlated wind directions in psfweather-station.
For each of the three cases (psfws, bench, and match), we simulate one 30-second, 3.5-deg exposure (the expected exposure time and field of view for Rubin LSST) for each of the 531 time points in the six months of ERA5 and site telemetry data.Each triplet of simulations has the same outer scale (drawn from a truncated log normal distribution with median of 25 m and used for all phase screens), the same random seed for realizing the turbulent phase variations in the phase screens, and the same atmospheric seeing (drawn uniformly from 0.6 to 1.6 arcsec).The contribution of seeing from each phase screen varies according to the turbulence integrals used for that case (psf-weather-station outputs for psfws; randomized Ellerbroek (2002) for bench and match).
PSF images are generated at 50k random locations across the field of view with a pixel resolution of 0.2 arcsec.Each PSF is drawn with sufficient photons (10 6 ) such that Poisson fluctuations are not significant and then convolved with a Gaussian of 0.35 arcsec FWHM to account for the PSF contribution from optics and sensors.To avoid issues related to overlapping PSF images, each PSF is generated and measured individually on a 50 × 50 pixel grid.

PSF PARAMETERS AND TWO-POINT STATISTICS
We estimate PSF size and shape from the (weighted) second moments Q ij of each PSF intensity profile I(θ x , θ y ): where θ x and θ y correspond to angular position on the focal plane and W (θ x , θ y ) is a weighting function.We have used the GalSim implementation of the HSM adaptive moments algorithm (Hirata & Seljak 2003) to measure PSF Q ij .As a measure of PSF size, we use For PSF shape, we use a definition of ellipticity e that is commonly used in weak lensing analyses10 : The magnitude of e is given by |e| = 1−q 2 1+q 2 , where q is the ratio of the minor to major axis of the second-moment ellipse.If e 2 = 0, then the orientation of the major axis of the ellipse is parallel to or perpendicular to the θ x direction for positive or negative values of e 1 , respectively.If e 1 = 0, then the major axis of the ellipse lies parallel to or perpendicular to an axis rotated +45 • with respect to the θ x axis, for positive or negative values of e 2 , respectively.Non-zero e 1 and e 2 describe orientations in between.
This (complex) ellipticity parameter e has a welldefined response to lensing shear when averaged across an ensemble of galaxy images if the effect of the PSF on the image has been accurately removed.Errors in the model for the PSF size and shape result in multiplicative and additive shear biases, respectively, and the exact impact on the ensemble weak lensing shear observables also depends on their spatial correlations.11 Spatial two-point correlation functions (2PCFs) for both PSF size and shape are relevant in weak lensing analyses.We define the PSF size two-point correlation function as where δσ is the deviation in PSF size σ from the mean size in that exposure, the indices a and b denote a pair of PSFs, and the angle brackets indicate an average over all pairs at each angular separation θ in the field of view.
For calculations of two-point correlation functions for PSF shape parameters, it is most useful to define the complex ellipticity with respect to the separation vector between each pair of PSFs: where the tangential and cross components e t and e × play the roles of e 1 and e 2 , respectively, in Equation 8.In other words, the y axis in Equation 8 is now oriented along the separation vector and the x axis perpendicular to it.A pair of two-point correlation functions for PSF shape are defined in terms of the tangential and cross components of e: ξ ± (θ) = ⟨e t,a e t,b ± e ×,a e ×,b ⟩(θ). (11)

SIMULATION RESULTS AND COMPARISONS
For each type of simulation described in Section 4psfws, match, or bench -we expect the simulated turbulence to imprint structure in the distribution of PSF parameters across a given exposure.As an example of the output of the psfws simulation, we display in Figure 3 the spatial distribution of the PSF size relative to the mean size (top panel) and the PSF shape (|e| and orientation, bottom panel) of the simulated PSFs for a single exposure generated with one of the psfws parameter sets (the representative example described in Section 3.3 and depicted by the dark curves and dots in Figure 2).The PSF parameters are clearly correlated across the field of view and the correlations are not isotropic.
In this section, we quantify these correlations and their anisotropy for each type of simulation inputs.We compare ensemble statistics across exposures for each simulation type, with a focus on the spatial two-point correlations of PSF size and shape, and study the dependence on particular input weather parameters.

Variance in PSF parameters
In Figure 4, we plot the variance of PSF ellipticity e (defined in Equation 8) across each psfws12 simulated exposure versus the input ground-layer wind speed for that exposure.The variance decreases with increasing ground-layer wind speed v(GL); a similar negative correlation with v(GL) is observed for variance in PSF size.This is expected because wind moves the phase screens across the field of view, washing out variations in PSF size or shape due to turbulence structure at all angular scales; therefore, the higher the wind speed, the more the variations are suppressed.
Correlations between variance in PSF size or shape parameters and free-atmosphere wind speed (not shown) are, in general, found to be weaker.This is expected because, instantaneously, different PSFs across the focal plane are more likely to sample independent regions of the phase screen with increasing altitude (see illustration in Figure 1); the resulting variance in a PSF parameter is less likely to decrease due to motion of the screen, compared to the GL screen where different PSFs are more likely to sample overlapping regions of the phase screen.The range of θ (∼ 1 to 100 arcmin) is limited near the low end by the density of PSFs (10 4 per 4.4×10 4 arcmin 2 ) and near the upper end by the size of the field of view (≈ 210 arcmin).The range of angular separations of interest in cosmic shear analyses is also ∼ 1 to 100 arcmin (see, for example, Asgari et al. 2021;Amon et al. 2022;Li et al. 2023).
The values of ξ + (θ) at separations of ∼ 1 arcmin are of order 10 −4 , which is of the same order (or larger than) the values expected and measured for cosmic shear.Since most galaxies used to probe cosmic shear are of similar size to the PSF, this implies that PSF shapes and their variation across the focal plane must be modeled accurately to avoid significant bias on measures of cosmic shear.Errors in the modeling of PSF size and its variation across the focal plane can also lead to biased shear 2PCFs (Rowe 2010;Jarvis et al. 2016).
For all three types of inputs to the simulations (psfws, bench, and match), the values of ξ − (θ) are approximately two orders of magnitude lower than those for ξ + (θ)), with median values that are slightly positive.Ensemble median values are shown as curves for three different sets of inputs to the simulations, described in Section 4. The psf-weather-station input parameters (wind speed, wind direction, and turbulence contribution) are correlated between phase screens at different altitudes (psfws, solid red); the benchmark input parameters are not correlated between phase screens (bench, long-dashed blue); and the matched input parameters are the same as bench but with phase-screen wind directions matched to those in psf-weather-station (match, short-dashed yellow).Shaded areas depict region between 25th and 75th percentile values for psfws simulation inputs (i.e., central 50 percentile values).

2D two-point correlation functions
In order to describe and quantify anisotropies in the distribution of PSF parameters across the focal plane (such as those that are evident for PSF size and ellipticity in Figure 3, and repeated for ellipticity in the top panel in Figure 6), we introduce the angle α to describe the polar angle of the separation vector ⃗ θ between the location of two PSFs, measured with respect to the θ y axis, as illustrated in the top panel in Figure 6.We will specify α in degrees (from 0 to 180 • ) and will continue to specify angular separation θ in arcmin.
Since the magnitude of the 2PCF ξ + is approximately two orders of magnitude greater than that for ξ − , we focus on quantifying anisotropies in ξ + .In the lower three panels in Figure 6, the color scale is constant and depicts the value of ξ + for the single exposure simulated with the representative example of psf-weather-station inputs used to generate the top panel.The first of the lower three panels shows the dependence of ξ + on ∆θ x , ∆θ y . 14On all scales, ξ + (∆θ x , ∆θ y ) is largest for pairs of PSFs with a separation vector oriented at an angle α ≈ 110 • .In addition, for small separations (θ ∼ 1 arcmin), ξ + is enhanced for pair orientations with α ≈ 75 • .
In the lower two panels in Figure 6, we show the same 2PCF but now in "polar" coordinates -i.e., as a function of PSF pair separation θ and orientation angle α -for two ranges of θ: 20 to 120 arcmin, and 1.5 to 6 arcmin.In each plot we see a dark vertical band corresponding to the maximum values of ξ + (θ, α) at orientations consistent with the above estimates of α ≈ 110 • and ≈ 115 • for large and small scales, respectively.For the same exposure, similar features are observed at the same angles α in plots of the 2PCF for PSF size, C(θ, α), for large and small separations.

PARAMETERS
In this section, we probe the relationship between simulation input wind directions and the directions α along which the anisotropic 2PCFs for the output PSF parameters are maximum.To quantify the orientation of the anisotropy in ξ + (θ) at large and small separations, we identify the value of α for which ξ + (θ, α) is maximum (denoted by α max ) for two ranges of separation: θ between 108 and 120 arcmin and between 3 and 6 arcmin, approximately coinciding with the largest and smallest scales currently used in cosmological analyses (e.g., Amon et al. 2022).These ranges correspond to the region above the white dashed line in each of the lower two panels in Figure 6.We find α max for each separation range for each of the 531 simulated exposures for each of the three types of simulation inputs (psfws, bench, match) described in Section 4.
For wind direction, we use both the orientation ϕ(GL) of the ground-layer wind velocity and the orientation ϕ(FA) of the wind velocity in the free-atmosphere layer with highest turbulence strength.
In the scatter plots in Figure 7, we plot α max for the shape 2PCF ξ + (α) versus ϕ(GL) (left) and versus ϕ(FA) (right), for separations of ∼ 100 arcmin (top) and ∼ 1 arcmin (bottom).The projected histograms in the top row illustrate how the distributions of GL and FA wind orientations are uniform for bench simulations and peaked at the same dominant wind orientations for psfws and match simulations. 15 The Pearson correlation coefficient for each set of sim-
ulations is displayed on each scatter plot.A strong correlation between α max for separations of ∼ 100 arcmin and ϕ(GL) (top left scatter plot) exists for all three types of simulation inputs.In contrast, for separations of ∼ 1 arcmin , there is no significant evidence for correlation between α max and ϕ(GL).
The results are quite different for the FA wind direction; there are no significant correlations at separations of ∼ 100 arcmin, but correlations exist at separations of ∼ 1 arcmin with some variation between the three types of inputs.In particular, the correlation is very strong for simulations with realistic altitudedependent FA wind directions and turbulence strengths (psfws, ρ = 0.81 ± 0.04) but is lower when these parameters are chosen randomly with altitude (bench, ρ = 0.47±0.05).When input FA wind directions (but not turbulence strength) are matched to those of psfws, the correlation coefficient has an intermediate value (match, ρ = 0.63±0.05).Although the correlation coefficients are statistically quite significant at separations of ∼ 1 arcmin when FA wind directions are correlated across altitudes, the distributions of α max have higher variance than those at separations of ∼ 100 arcmin.
To understand the FA results, we describe the expected angular scale of PSF correlations and a physical source of anisotropies.
Angular scales: Atmospheric features of a particular physical size that are located at lower altitudes have a larger angular size in the focal plane than features of the same physical size at higher altitudes.Therefore, at one instant in time an optical turbulence pattern near the ground results in spatial variations in the PSF over larger angular separations than the same turbulence pattern located near the top of the free atmosphere.Consider, for example, turbulence with a scale equal to the median outer scale L 0 .(The power spectrum flattens for spatial frequencies lower than ∼ 1/L 0 , as shown in Equation 1.)For L 0 = 25 m and a ground-layer screen at a height of 200 m, the angular scale of turbulence variations is arctan 25 200 ≈ 400 arcmin; for a phase screen at an altitude of 20 km, the angular scale is arctan 25 20000 ≈ 4 arcmin.This illustrates why the ground layer is more relevant at large angular separations and the free atmosphere is more relevant as small separations.
Sources of PSF anisotropies and correlations with wind: While wind speeds and directions can change throughout the night, they are typically fairly stable during the course of an exposure.This is particularly true at Cerro Pachón, where the wind direction is persistently coming from the sea.During a 30-sec exposure then, we expect air at all altitudes to be moving coherently; different layers may be translating at different speeds, but in mostly the same direction.As a consequence, PSF images along the direction of the wind on the focal plane have "seen" much of the same optical turbulence, although in slightly different combinations due to wind speed variation with altitude.The shapes of these PSFs will thus be more correlated with each other than with those in a direction orthogonal to the wind.-Orientation angle αmax at which the anisotropic two-point correlation function ξ + for the PSF shape is maximum versus input wind direction ϕ for each simulated exposure, for the three types of simulation inputs (bench, match, and psfws).The plots in the middle row correspond to large angular separations (θ between 108 and 120 arcmin) and those in the bottom row to small angular separations (θ between 3 and 6 arcmin).The angles ϕ(GL) and ϕ(FA) correspond to the wind velocity direction in the ground layer (left column) and in the highest-turbulence-strength free-atmosphere layer (middle column), respectively, modulo 180 • .The distributions of the angles are shown in the four projected histograms.Correlation coefficients for each type of simulation are reported on each scatter plot (bench in blue, match in gold, psfws in red).
The realistic inputs used to generate the psfws exposures include these correlated wind speeds and directions.In the case of bench simulations, there is no coherent motion of the turbulence in the FA because the wind speed and wind direction for each layer is chosen randomly; therefore, the direction of the highest-turbulence layer is not related to that of the other FA layers.This results in a correlation coefficient ρ between FA wind direction and α max that is suppressed relative to psfws at separations of ∼ 1 arcmin.The match simulations have wind directions that are correlated across altitudes, but speeds are chosen randomly.The random speeds cause turbulence layers to move with greater differences in speed than the smoothly varying wind profiles in psfws, resulting in a slightly suppressed correlation coefficient for match simulations despite the wind directions being matched with psfws.
These conclusions hold when considering the ensemble of FA layers: we found similar results, albeit with less correlation, when using the sum of the turbulenceweighted velocities of all FA layers to define ϕ(FA), rather than choosing the velocity of the single layer with highest turbulence.
Similar correlations are also observed between wind direction and the orientation of anisotropies for PSF size 2PCF C(α), although the values of the correlation coefficients between α max and ϕ(GL) at separations of ∼ 100 arcmin are lower than observed for PSF shape, potentially due to noisier estimates of α max .

CONCLUSIONS AND FUTURE WORK
As described in the introduction, accurate measures of cosmic shear with future astronomical surveys, such as LSST at Rubin Observatory, require unbiased measures of two-point statistics for galaxy shapes, which in turn require unbiased measures of the size and shape of the PSF across the field of view.In this work, we use realistic, altitude-dependent weather and turbulence input provided by the psf-weather-station package (Figure 2) to simulate the atmospheric PSF across the Rubin Observatory field of view.We summarize our findings here: 1.The variance in PSF size and shape across a single exposure decreases as wind speed increases (see Figure 4).
2. The values of the 2PCFs for PSF shape in a single exposure are of the same order as (or larger than) the expected 2PCFs for cosmic shear over the range of angular separations used in cosmic shear analyses: from a few arcmin to over 100 arcmin (see Figure 5).
3. There exist dominant wind directions at Cerro Pachón (red histograms in top panels in Figure 7), which in turn lead to dominant orientations of anisotropies in the 2PCF ξ + (red histograms in right panels in Figure 7).At scales of ∼ 100 arcmin, these anisotropies are due to strong correlations with ground-layer wind direction (upper left scatter plot in Figure 7); at scales of only a few arcmin, they are due to correlations with free-atmosphere wind direction (lower right scatter plot).As discussed in Section 7, these results can be understood in terms of the different angles subtended at different heights by turbulent structure of the same physical scale.
PSF modeling and interpolation methods must accurately capture the anisotropic two-point correlations for PSF size and shape on different scales.In the future, high fidelity simulations generated with psf-weatherstation can be used to test whether current modeling and interpolation methods (e.g., those implemented in Jarvis et al. 2020) can reach the necessary accuracy.One interpolation technique that merits further exploration is anisotropic Gaussian process interpolation; see, for example, Léget et al. (2021) and Fortino et al. (2021) for applications to PSF astrometry.
Because of the dominant wind direction at Cerro Pachón, we expect to see a dominant orientation of the anisotropy in the 2PCF for PSF size and shape with respect to the ground coordinate system, across LSST exposures at Rubin Observatory.The mapping of wind direction on the ground onto sky coordinates is determined by the pointing of the telescope for each exposure.Therefore, the degree to which the dominant wind direction will vary on the sky for a single field depends on the observing strategy for the survey.Using the observing strategy for the 300-square-degree DESC DC2 simulation, we find that the dominant wind direction from Figure 7 will translate to a persistent on-sky anisotropy; further study is needed to understand the implications for the full LSST survey.PSF simulations produced with psf-weather-station input can be used to study this question for a particular survey strategy.
The psf-weather-station software package for producing correlated weather and turbulence input to simulations, configurable to any observatory, is public at https://github.com/LSSTDESC/psf-weather-station and includes installation instructions, documentation, and tutorial notebooks.2019).We acknowledge ECMWF for access to the ERA5 data through the MARS access system.The computing for this project was performed on the Stanford Sherlock cluster.We would like to thank Stanford University and the Stanford Research Computing Center for providing computational resources and support.

AUTHOR CONTRIBUTIONS
Claire-Alice Hébert developed the psf-weatherstation software package, performed the main analysis, and produced all the plots in the paper.Patricia Burchat contributed to defining the project, and advised throughout.C-AH and PB contributed equally to writing the paper.Joshua Meyers advised on all aspects of the project and reviewed the psf-weather-station package.My H. Do contributed to the analysis early in the project as a participant in the Cal-Bridge program for undergraduate students.

GLOBAL CIRCULATION MODELS
As summarized in Section 3.2, multiple organizations around the world produce high-quality weather models and forecasts; we focus here on those from the European Centre for Medium-Range Weather Forecasts (ECMWF) and the National Oceanic and Atmospheric Administration National Centers for Environmental Prediction (NOAA NCEP).Data products from these global models can be very useful for studies of the atmosphere for astronomical applications.Here we summarize the types of data available and considerations for their use in atmospheric PSF simulations.
Both ECMWF and NCEP make available two types of data: (1) analysis products are the best estimate of the state of the atmosphere, produced by combining a numerical weather prediction model with a variety of observations through a process called data assimilation; (2) forecast products are the numerical predictions (based on initial analysis products) for some time into the future.Analysis and forecast data are available in real-time (of use for weather forecasting) and as reanalysis data products: state-of-the-art data assimilation and numerical modeling methods applied to archival data (highly relevant for long-term climate monitoring). 16 The 5th generation ECMWF reanalaysis (ERA5) catalog covers the time period from 1940 to the present and is extensively documented 17 (Hersbach et al. 2020).ERA5 analyses are available hourly, with forecasts initialized at 00h and 18h UTC.At the time of writing, all ECMWF archival data (including ERA5) and subsets of real-time forecasts are available publicly under creative commons. 18 All NOAA NCEP data are available publicly; real-time forecasts are from the Global Forecast System (GFS), and several reanalysis efforts are available for specific time periods. 19The NCEP analysis products are available every 6h.
Data products are output on a uniform grid over the Earth's surface, with resolution 31 km ≈ 0.28 • for ERA5 and NCEP GFS, and 2.5 • for NCEP reanalysis.
Data products are available at heights corresponding to specific levels of pressure rather than specific altitudes.Since higher spatial resolution in the vertical direction allows for more accurate capture of important wind gradients in the atmosphere, we use the output type with densest vertical coverage -called model levels.ECMWF uses 137 model levels.NCEP uses 127 for the period since February 3, 2021, and 64 prior to that time.Model levels follow terrain at the Earth's surface, and the conversion to altitude uses temperature and specific humidity.This conversion has been implemented in psf-weatherstation following ECWMF documentation and example code found in the Q&A section of the ERA5 wiki.
We choose to use ECMWF data products because the temporal and spatial resolution of available data is higher than for NCEP.In addition, ECMWF documentation is more detailed and accessible.

Fig. 2 .
Fig. 2.-Six months of wind data at and above Cerro Pachón (May 2019 through Oct 2019), processed with psf-weatherstation as described in Section 3. We plot wind speed (top) and meteorological direction (direction of wind origin; middle) as a function of altitude, and as a frequency distribution (right).In the top two rows, weather tower measurements near the ground at Gemini South are shown in orange and ECMWF ERA5 data for the free atmosphere are shown in purple.The heavy purple line in each panel corresponds to data from a representative example time (August 6 2019, 08h UTC).In the bottom panel, the uncalibrated C 2 n (h) profile for the same time is shown in light purple; the calibrated profile (scaled by J FA ) is shown in dark purple.The dashed vertical lines depict the boundaries between the altitude bins used to calculate turbulence integrals J for each FA phase screen.The dots correspond to the value of J divided by the corresponding range of altitudes and are placed at 200 m above the elevation of the observatory (orange dot) and the C 2 n -weighted position within each FA layer (purple dots).

Fig. 3 .
Fig. 3.-Spatial distributions of PSF size relative to mean size (top) and PSF ellipticity (bottom) across a simulated 3.5-deg 2 exposure generated with the psf-weather-station input parameters for the representative example depicted by the dark curves and dots in Figure 2. The orientation of each line in the bottom plot corresponds to the orientation of the major axis of the PSF shape, while the length and color contrast are proportional to the magnitude of the ellipticity e defined in Equation 8.

Fig. 4 .
Fig. 4.-Variance of PSF ellipticity e for each simulated exposure versus the ground layer (GL) wind speed in that simulation, for the psfws simulations described in Section 4. The Pearson correlation coefficient ρ is reported on the scatter plot.

6. 2 .
1D two-point correlation functions In Figure 5, we show the value of the two-point correlation function (2PCF) for PSF size (C(θ), top) and PSF ellipticity (ξ + (θ), middle; ξ − (θ), bottom), as a function of angular separation θ between pairs of PSFs. 13The ensemble median values of the 2PCFs are shown as curves for each type of simulation input: psfws (solid red), bench (long-dashed blue), and match (short-dashed yellow).The shaded areas depict the region between the 25th and 75th percentile values of the 2PCFs for psfws simulation inputs (i.e., the central 50 percentile values).

Fig. 5 .
Fig. 5.-Values of the two-point correlation functions for PSF size (C(θ), top) and PSF ellipticity (ξ + (θ), middle, and ξ − (θ), bottom), as a function of angular separation θ between pairs of PSFs.Ensemble median values are shown as curves for three different sets of inputs to the simulations, described in Section 4. The psf-weather-station input parameters (wind speed, wind direction, and turbulence contribution) are correlated between phase screens at different altitudes (psfws, solid red); the benchmark input parameters are not correlated between phase screens (bench, long-dashed blue); and the matched input parameters are the same as bench but with phase-screen wind directions matched to those in psf-weather-station (match, short-dashed yellow).Shaded areas depict region between 25th and 75th percentile values for psfws simulation inputs (i.e., central 50 percentile values).

Fig. 6 .
Fig. 6.-Illustration of anisotropic two-point correlation function of PSF shape for the psfws representative example.The top panel shows the simulated PSF ellipticity across the field of view overlain with a diagram of an example pair of PSF locations (black dots) and the separation between them, in Cartesian and polar coordinates.The color scale in the second panel shows the dependence of the 2PCF ξ + on the coordinates ∆θx and ∆θy.The color scales in lower two panels show the same 2PCF transformed to polar coordinates (separation θ versus angle α) for two ranges of θ: 20 to 120 arcmin, and 1.5 to 6 arcmin.The range above the white dashed line in each of the two bottom panels is referenced in Section 7.

Fig
Fig.7.-Orientation angle αmax at which the anisotropic two-point correlation function ξ + for the PSF shape is maximum versus input wind direction ϕ for each simulated exposure, for the three types of simulation inputs (bench, match, and psfws).The plots in the middle row correspond to large angular separations (θ between 108 and 120 arcmin) and those in the bottom row to small angular separations (θ between 3 and 6 arcmin).The angles ϕ(GL) and ϕ(FA) correspond to the wind velocity direction in the ground layer (left column) and in the highest-turbulence-strength free-atmosphere layer (middle column), respectively, modulo 180 • .The distributions of the angles are shown in the four projected histograms.Correlation coefficients for each type of simulation are reported on each scatter plot (bench in blue, match in gold, psfws in red).
from the Institut National de Physique Nucléaire et de Physique des Particules in France; the Science & Technology Facilities Council in the United Kingdom; and the Department of Energy, the National Science Foundation, and the LSST Corporation in the United States.DESC uses resources of the IN2P3 Computing Center (CC-IN2P3-Lyon/Villeurbanne -France) funded by the Centre National de la Recherche Scientifique; the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231; STFC DiRAC HPC Facilities, funded by UK BEIS National E-infrastructure capital grants; and the UK particle physics grid, supported by the GridPP Collaboration.This work was performed in part under DOE Contract DE-AC02-76SF00515.Generated using Copernicus Climate Change Service information (