The impact of high spatial frequency atmospheric distortions on weak lensing measurements
Catherine Heymans, Barnaby Rowe, Henk Hoekstra, Lance Miller, Thomas Erben, Thomas Kitching, Ludovic Van Waerbeke
aa r X i v : . [ a s t r o - ph . C O ] D ec Mon. Not. R. Astron. Soc. , 000–000 (0000) Printed 30 October 2018 (MN L A TEX style file v2.2)
The impact of high spatial frequency atmospheric distortions onweak lensing measurements
Catherine Heymans, ⋆ Barnaby Rowe, , ⋆ Henk Hoekstra, Lance Miller, Thomas Erben, Thomas Kitching, & Ludovic Van Waerbeke Scottish Universities Physics Alliance, Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ, UK. Department of Physics & Astronomy, University College London, Gower Street, London, London, WC1E 6BT, UK. California Institute of Technology, 1200 E California Boulevard, Pasadena, CA, 91106, USA. Leiden Observatory, Leiden University, Niels Bohrweg2, NL-2333 CA Leiden, The Netherlands. Department of Physics, Oxford University, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK. Argelander-Institut f¨ur Astronomie, Universit¨at Bonn, Auf dem H¨ugel 71, 53121 Bonn, Germany. University of British Columbia, 6224 Agricultural Rd., Vancouver, BC, V6T 1Z1, Canada.
30 October 2018
ABSTRACT
High precision cosmology with weak gravitational lensing requires a precise measure of thePoint Spread Function across the imaging data where the accuracy to which high spatial fre-quency variation can be modeled is limited by the stellar number density across the field. Weanalyse dense stellar fields imaged at the Canada-France-Hawaii Telescope to quantify thedegree of high spatial frequency variation in ground-based imaging Point Spread Functionsand compare our results to models of atmospheric turbulence. The data shows an anisotropicturbulence pattern with an orientation independent of the wind direction and wind speed. Wefind the amplitude of the high spatial frequencies to decrease with increasing exposure timeas t − / , and find a negligibly small atmospheric contribution to the Point Spread Functionellipticity variation for exposure times t > seconds. For future surveys analysing shorterexposure data, this anisotropic turbulence will need to be taken into account as the ampli-tude of the correlated atmospheric distortions becomes comparable to a cosmological lensingsignal on scales less than ∼ arcminutes. This effect could be mitigated, however, by cor-relating galaxy shear measured on exposures imaged with a time separation greater than 50seconds, for which we find the spatial turbulence patterns to be uncorrelated. Key words: cosmology: observations - gravitational lensing
Over the next decade we will see all sky surveys being undertakento address key science questions ranging from cataloging all nearearth objects to mapping the Dark Universe. These surveys includethe Large Synoptic Survey Telescope (LSST Science Collaboration2009, LSST), Euclid (Laureijs et al. 2011), and the Panoramic Sur-vey Telescope and Rapid Response System (Kaiser et al. 2010,PanSTARRS). By observing the growth of dark matter structuresover cosmic time, weak lensing may provide the most interest-ing constraints on the dynamical properties of dark energy. Along-side modifications to General Relativity, this dark energy is hy-pothesized as an explanation for the accelerated expansion of theUniverse (see for example Huterer 2010, and references therein).The observational measurement of weak gravitational lensing ishowever far from trivial. In fact it is the technique that sets the ⋆ email: [email protected] & [email protected] most stringent requirements on the optical design of future instru-mentation as we require accurate knowledge of the Point SpreadFunction (PSF). The PSF also needs to be minimal in terms ofits size and ellipticity. Dark matter typically induces a one per-cent change in the ellipticity of lensed distant galaxies. It isthis distortion that we measure to constrain cosmology. In con-trast the telescope, camera and atmosphere can introduce imagedistortions up to the ten percent level. The main challenge inweak lensing is therefore to distinguish the cosmological distortionfrom instrumental and atmospheric distortions in order to recoverthe underlying dark matter signal. A significant effort is under-way to improve the accuracy of weak lensing measurements bothfrom a software (Kitching et al. 2010) and hardware perspective(LSST Science Collaboration 2009). Furthermore ambitious space-based surveys are proposed in order to remove atmospheric seeingand obtain diffraction limited observations (Laureijs et al. 2011).For ground-based observatories, the varying refractive indexin the turbulent layers of the atmosphere produces two effects. It c (cid:13) C. Heymans & B. Rowe et al. blurs the images of point sources and also introduces an instanta-neous ellipticity at the 10% level (de Vries et al. 2007). This canbe viewed a superposition of speckles. As the exposure time in-creases, these speckles de-correlate producing a circular distortionpattern. Telescope optics and the camera distortion then induce anadditional smoothing term and an anisotropic distortion across thecamera field of view. The combined observed PSF size is roughlygiven by adding the size of the atmospheric PSF in quadraturewith the size of the instrument PSF. Optimal instrumentation de-sign therefore aims to minimize the instrument contribution so thatthe resulting PSF is atmosphere dominated. Image distortions areencompassed by the PSF as imaged by stellar point sources. Usingcolour and size information, stars can be separated from galaxiesand used to model the PSF. For weak lensing measurements to date,the PSF model is typically an interpolation between stars, using alow order polynomial model to obtain a PSF at each galaxy loca-tion. This assumes however that the PSF does not vary rapidly onsmall scales (Paulin-Henriksson et al. 2008).In this paper we measure and analyse the image distortion in-duced in short exposure observations by atmospheric aberrations.This is of particular relevance to the PanSTARRS and upcomingLSST ground-based surveys. Both these surveys will image morethan 20,000 sq. degrees with LSST achieving a depth of r ∼ . over 10 years. As transient objects are a key science goal for boththese ambitious surveys repeat short exposure observations are re-quired. For LSST each region of sky will be imaged close to 1000times with 15 second exposures. The current Pan-STARRS surveyimages each region of sky 5 times with 30 second exposures. Asboth surveys progress a deep data set will be collected for opti-mal weak lensing analysis, but with such short exposure times thequestion is whether the rapidly varying atmospheric ellipticity willdominate the PSF degrading the ability to model the PSF in indi-vidual exposures using standard methods.Wittman (2005) first looked at this issue by analysing a set of10 second and 30 second exposure images of a dense stellar fieldfrom the Subaru Telescope taken in one night. The conclusion fromthis paper was that the spurious small scale power induced by theatmosphere would not be a significant factor in the co-added dataof future massive surveys based on the finding that the atmosphericdistortion in exposures separated by ∼ seconds was found to beuncorrelated. This analysis was extended by Asztalos et al. (2007)with an analysis of 15 second exposure images from the GeminiSouth Telescope where evidence was presented for a correlation be-tween the atmospheric distortion and the ground-based wind speedand direction with a consequence that the co-addition of exposurestaken with the same observing site prevailing wind may not de-correlate as concluded by Wittman (2005).An investigation of the image ellipticity of atmosphericaberrations in simulated short exposure data is presented inde Vries et al. (2007). Using ray-tracing and Fourier methods theysimulate a turbulent atmosphere with a single layer Kolmogorovpower spectrum, finding significant 10% ellipticities in the instanta-neous atmospheric PSF. They also investigate the time dependenceof the atmospherically induced ellipticities. Assuming a constantwind speed they show the amplitude of the atmospheric ellipticitydecreases with time t as t − / . This simulation work has been ex-tended by Jee & Tyson (2010) with simulated LSST short exposureimages with a PSF that includes a six layer Kolmogorov powerspectrum turbulence component and a ray-tracing optical distor-tion model. They recover the de Vries et al. (2007) result that asexposure time increases, the atmospheric PSF becomes rounder.Using a Principal Component Analysis (PCA) they model the simu- lated short-exposure PSF across the LSST field-of-view and, keep-ing only the most significant 20 eigen-modes of that PCA analysis,they simulate what is described as a realistic PSF-convolved LSSTimage. They conclude that a typical stellar density is sufficient tomodel the PCA derived PSF and demonstrate their ability to removethat PSF distortion from galaxy shape measurement. By limitingthe PSF to only 20 eigen-modes however, this analysis suppressesthe random and high spatial frequency variation that we investigatein this study.In this paper we assess the impact of atmospheric distortionson weak lensing surveys by analysing short exposure images ofhigh density stellar regions with the wide-field MegaCam Imageron the Canada-France-Hawaii Telescope (CFHT). It is of partic-ular importance to revisit this issue as a significant improvementin the accuracy of weak lensing shape measurement has recentlybeen found when data is analysed using individual exposures com-pared to the analysis of stacked exposures (Miller et al. in prep).Weak lensing data in the future will therefore not be analysed froma co-addition of exposures where all atmospheric turbulent effectshave de-correlated, as it has been to date. With a wide range ofexposure time, wind speeds and directions we test the findings ofde Vries et al. (2007) that the atmospheric ellipticity decreases withtime t as t − / and also the findings of Asztalos et al. (2007) thatthe atmospheric ellipticity is correlated with the wind speed and di-rection. Our results complement efforts to accurately simulate theseeffects (Chang et al. in prep) and inform the hardware and softwaredesign of future ground-based surveys accordingly. This paper isorganised as follows. In Section 2 we describe the data and anal-ysis methods, presenting the results and analysis in Section 3 andconclusions in Section 4. We analyse the CFHT PSF using a series of observations of densestellar fields imaged and archived by the CFHT and listed in Ta-ble 1. We select exposures taken in sub-arcsecond seeing conditionsin the i -band over a range of exposure times. Our main sample of60 exposures imaged in 74 seconds were taken during three dif-ferent observing runs sampling a wide range of wind speeds anddirections. We also analyse nine other samples with exposure timesranging from 1 second to 450 seconds with each set observed ondifferent nights.The data were processed and stellar object catalogues pro-duced as described in Hoekstra et al. (2006) yielding a typical stel-lar density of 7 stars per square arcminute, and in all datasets asignificantly greater stellar density than available in extragalacticimaging survey data taken out of the galactic plane. We parameter-ize the PSF in terms of stellar ellipticity as measured by weightedquadrupole moments Q ij = R d θ W ( θ ) I ( θ ) θ i θ j R d θ W ( θ ) I ( θ ) , (1)where I ( θ ) is the surface brightness of the star, θ is the angulardistance from the star centre and W is a Gaussian weight functionof scale length r g , where r g is the measured dispersion of the PSF(Kaiser et al. 1995). For a perfect ellipse, the weighted quadrupolemoments are related to the weighted ellipticity parameters ε α by (cid:18) ε ε (cid:19) = 1 R (cid:18) Q − Q Q (cid:19) , (2) c (cid:13) , 000–000 tmospheric Distortions Exptime N Exp Proposal ID Date1.0 28 08AH22 2008-04-14/05-072.0 7 03BF99 2004-01-19/2304BF28 2004-12-05/1203BF09 2003-11-18/12-167.5 3 08AQ98 2008-05-09/06-0310.0 33 06BH48 2006-08-19/2230.0 9 03BF99 2004-01-19/2345.0 5 03BQ97 2003-09-3074.0 60 08AH22 2008-04-14/05-08180.0 9 03BL06 2003-10-03300.0 11 08AH53 2008-02-13/14450.0 12 07BC02 2007-11-14/01-10
Table 1.
Table of CFHT archive data used in this analysis listing the ex-posure time, total number of exposures used, the CFHT Proposal ID whichcan be used to query the CFHT Science Data Archive, and the range ofobservation dates. where R is related to object size and given by R = p Q + Q , (3)If the weight function W ( θ ) = 1 in equation 1, the ellipticity orpolarisation | ε | = (1 − β ) / (1 + β ) , where β is the axial ra-tio of the ellipse (see Bartelmann & Schneider 2001). We note thatfor the purposes of this study the quadrupole moment descriptionof the PSF that we adopt is sufficient. For a detailed weak lensinganalysis however, many techniques now determine either a set oforthonormal 2D basis functions or a 2D pixel-based model to de-scribe the PSF (Kitching et al. 2010). Both these methods wouldbe more challenging to model in the presence of high order spatialfrequency variation.In each exposure the variation of the PSF ellipticity ε PSF i andsize R across the field of view is considered to have two compo-nents; a smoothly varying second order polynomial over position, ε model i ( x, y ) and R ( x, y ) for each chip and higher spatial varyingresiduals with δε i = ε PSF i − ε model i , (4) δR = R − R . (5)With a typical number density of stellar objects imaged in an high-Galactic latitude field (40 per MegaCam chip, compared with over700 in this analysis), the chip-wise second order polynomial modelwould be the most complex model that could be accurately fit tolensing survey data (Rowe 2010). Figure 1 shows a typical PSF pattern for a 74 second exposure asimaged by MegaCam at CFHT. The field of view is a square degreewith a pixel scale of . arcseconds and an × CCD patternwhich is visible. The left and right panels show the variation in theamplitude of the average PSF ε and ε respectively across the fieldof view. Each grid point in the greyscale map contains on average 5 Figure 1.
A typical PSF pattern for a 74 second exposure. The left and rightpanels show the variation in the amplitude of the average PSF ε and ε respectively across the field of view. The upper panels show the observedellipticity variation, the middle panels show the second order polynomialmodel fit to the data. The lower panels reveal the residual ellipticities δε and δε components. The residuals show high spatial frequencies and apreferred direction. The ellipticity amplitude in each panel is indicated bythe greyscale shown in the vertical colour bar. stars and spans . × . arcmins. In an image of a typical field, outof the Galactic plane, we would find ∼ useable star in every ∼ grid points ( ∼ . stars per sq arcmin). The upper panels show theobserved ellipticity variation. The middle panels show the secondorder chip-wise polynomial model fit to the data. The lower panelsreveal the residual ellipticities δε and δε components. Figure 2shows the size variation for the same exposure with the size resid-uals (right hand panel) showing the same high spatial frequenciesand a preferred direction as the ellipticity residuals. For the restof the paper we focus mainly on the ellipticity variation, but comeback to the size variation in Section 4 where we show the high spa-tial frequency variation of PSF size is as detrimental to the weaklensing shape measurement as the variation in ellipticity.In order to investigate the high spatial frequency variation ofthe PSF we measure the two-point correlation function of the resid-ual PSF ellipticities in Figure 3, showing the average systematicresidual PSF correlation functions ξ + (upper two panels) and ξ − (lower panel) where ξ sys ± ( θ ) = h ε t ( θ ′ ) ε t ( θ ′ + θ ) i ± h ε r ( θ ′ ) ε r ( θ ′ + θ ) i (6) c (cid:13) , 000–000 C. Heymans & B. Rowe et al.
Figure 2.
A typical PSF size pattern for the same 74 second exposure shownin Figure 1. The left hand panel shows the variation in PSF size R acrossthe field with the greyscale in arcseconds. The right hand panel shows theresidual variation in size after a second order polynomial model has beenremoved. The ripple pattern in size follows the same structure as the ellip-ticity residuals. The right hand greyscale is in units of − arcsec . and ε t,r are the tangential and rotated ellipticity parameters rotatedinto the reference frame joining each pair of correlated objects. Theaverage is taken over all exposures in our sample split by exposurelength. Figure 3 shows the amplitude of the residuals for four setsof decreasing exposure time revealing a characteristic shape to thecorrelation function that we find for all exposure times. ξ + is pos-itively correlated θ < ′ , anti-correlated ′ < θ < ′ , then pos-itively correlated until it becomes consistent with zero θ > ′ .This characteristic shape is also seen in the results of Wittman(2005). ξ − is measured with significance only for the shortest expo-sures (shown lower panel) and is positively correlated for all scales θ < ′ .The residual PSF correlation functions measured betweenconsecutive s , s and s exposures were found to be consistentwith zero. We can therefore conclude that the atmospheric turbu-lence distortion de-correlates in < seconds where the timescaleis set by the CFHT MegaCam read-out, overheads, slew and ac-quisition time. This is in agreement with Wittman (2005) who seta de-correlation time < seconds, limited by the Subaru readtime. The von K´arm´an model for isotropic atmospheric turbulence pre-dicts a projected, two-dimensional power spectrum P ( l ) ∝ (cid:18) l + 1 θ (cid:19) − / , (7)where the angle θ defines the outer scale of the spectrum (e.g.Sasiela 1994). In this model, a physically-motivated modificationof the scale-free Kolmogorov spectrum, turbulent structures in im-ages of the sky become uncorrelated at separations greater than θ .As shown in Appendix A, the isotropic ξ + correlation function forsuch a spectrum may then be written as ξ + ( θ ) ∝ θ / K − / (2 πθ/θ ) , (8)where K ν ( x ) is the modified Bessel function of the second kind(Arfken & Weber 2005). Numerical approximation of these func-tions is simple via the useful series expansion of Kostroun (1980).The functional form of ξ − is given in equation (A4). This expres-sion is significantly more complicated than the ξ + case and wasfound to be prone to numerical instabilities.We now investigate whether the model given in equation (8)
1 10-1.0e-04-5.0e-050.0e+005.0e-051.0e-041.5e-042.0e-042.5e-043.0e-043.5e-04 ξ + t exp = 1 st exp = 7.5 s
1 10-2.0e-060.0e+002.0e-064.0e-066.0e-068.0e-061.0e-051.2e-051.4e-05 ξ + t exp = 74 st exp = 450 s
1 10 θ [arcmin] -5.0e-060.0e+005.0e-061.0e-051.5e-052.0e-052.5e-053.0e-053.5e-054.0e-05 ξ − t exp = 1 st exp = 7.5 s Figure 3.
The ξ + ( θ ) and ξ − ( θ ) correlation function estimates for the PSFellipticity data, with best-fitting von K´arm´an models overlaid on the ξ + re-sults and the corresponding ξ − prediction plotted over the measured ξ − results (see Section 3.1 and Appendix A). Upper panel: the t = 1 s (solidline: best-fitting model) and t = 7 . s (dashed line: best-fitting model) ξ + ( θ ) results. Middle panel: the t = 74 s (solid line: best-fitting model)and t = 450 s (dashed line: best-fitting model) ξ + ( θ ) results. Lower panel:the t = 1 s, and t = 7 . s , ξ − ( θ ) results (the longer exposure PSF patternsshowed ξ − consistent with zero on all scales). is able to reproduce the correlation function results seen in Fig-ure 3. We restrict ourselves to fitting only the ξ + model, whichcan be simply calculated and carries almost all the detectable sig-nal. We go on to compare predictions for ξ − based on these fittingresults to measurements of ξ − in the data. In Figure 3 we haveoverlaid the best-fitting von K´arm´an models to the ξ + correlationfunction estimates for each of the example t = 1 s, 7.5s, 74s &450s exposure data sets. The Levenberg-Marquardt algorithm (e.g.Press et al. 1986) was used to determine the maximum-likelihoodparameter fits to the overall amplitude of the expression in equation(8), and for the outer scale length θ . In the lower panel of Figure 3we use the best-fitting values of θ from the ξ + models to gener-ate a von K´arm´an prediction for the ξ − signal in the t=1s and 7.5sexposure data.It can be seen that the von K´arm´an model gives a reasonablefit to the inner slope of ξ + for the three t = 1 s, 7.5s & 74s ex-posure data sets in Figure 3. For these three sets we also foundbroadly consistent best-fitting values of the outer scale in the range θ = 2 . – . ′ . These results suggest a turbulent origin for thesmall-scale correlated ellipticity measurements in these image, butin all three cases the von K´arm´an model fails to capture both thetrough on scales θ ≃ ′ , and the second peak in correlation near ′ . The discrepancy is also marked in the comparison between the c (cid:13) , 000–000 tmospheric Distortions ξ − model and the data. Here a very distinct excess ‘bump’ is seen,again peaking around θ ≃ ′ in the 1s exposure time data. How-ever, as the amplitude of the ξ − correlation is consistently smallerthan that of ξ + , these results are noisier. We discuss potential ori-gins for these features in Section 4. For the t = 450 s data we founda model amplitude consistent with zero, indicative of the weaknessof the signal for this long exposure data. A key result from the simulations of de Vries et al. (2007) showedthat the ellipticity introduced by atmospheric aberrations decreaseswith time as t − / assuming a constant wind speed. The two-pointcorrelation function of the atmospheric variations ξ + is thereforeexpected to decrease in amplitude as t − . Figure 4 shows the aver-age two-point correlation function | ξ + | as a function of the expo-sure time. Results are shown for two representative angular scaleswith ξ + measured at θ = 0 . ′ (upper) where the von K´arm´an at-mospheric turbulence model provides a good fit to the data, and at θ = 4 . ′ which corresponds to the first negative dip in the correla-tion function. We find that the de Vries et al. (2007) model (showndot-dashed) is a good fit to the data but note that for this analysis,we were unable to select a significant sized sample of exposureswith a constant wind speed as is assumed in the de Vries et al.(2007) theory. We conclude that the scatter we see in the resultsmust be in part driven by the range of wind speeds within our sam-ple. Figure 4 also allows us to compare the amplitude of the at-mospheric turbulence to the expected cosmological signal in theCFHT Lensing Survey (CFHTLenS ). For both angular scales θ shown we over-plot two horizontal bars that indicate the ampli-tude of the two-point shear correlation function ξ + ( θ ) for two lowredshift tomographic bins with a mean redshift of z = 0 . and z = 0 . assuming a WMAP7 cosmology (Komatsu et al. 2011).We should note here that one cannot directly relate the PSF ellip-ticity correlation with that expected from shear, as the effect ofthe PSF on the shear measurement depends on the relative sizeof galaxy such that it would negligible for large galaxies or po-tentially amplified for galaxies comparable to or smaller than thePSF (Paulin-Henriksson et al. 2008). This comparison does how-ever demonstrate that even for the lowest redshift tomographic bin,the residual systematic error from atmospheric turbulence is wellbelow the cosmological signal for the 600 second CFHTLenS ex-posures. In order to compare the turbulence patterns observed to the winddirection we quantify the direction of the atmospheric turbulenceusing correlation as a function of vector separation and directionon the sky, rather than the commonly-used, azimuthally-averagedfunctions ξ + ( θ ) = ξ + ( | θ | ) . We calculate the ξ + residual correla-tion as a function of vector angular separation θ using the DiscreteFourier Transform to find the power spectrum of the δε and δε images (examples of which are shown in Figure 1), then employingthe Wiener-Khinchin theorem to generate images of ξ + ( θ ) . CFHTLenS analyses data from the CFHT Legacy Survey, using both theWide, Deep and Pre-imaging surveys
Figure 4.
The evolution with exposure time of the atmospheric turbulenceas measured by the two-point correlation function ξ + at two different angu-lar scales; 0.7 arcmin(black upper) and 4.9 arcmin (red lower). The best-fitde Vries et al.(2007) prediction is shown dot-dashed for each angular scale.This can be compared to the two pairs of thick horizontal bars over plot foreach angular scale. These bars indicate the expected CFHTLenS cosmolog-ical signal for two tomographic bins with a mean redshift z = 0 . and z = 0 . . determining the exposure time below which the cosmologicalsignal becomes lower than the residual atmospheric PSF signal. Figure 5.
The ξ + residual correlation as a function of 2D vector separation θ for four example 74 second exposures. The dominant direction of thePSF residuals can be compared to the wind direction shown as an arrowinset. The upper left panel shows an example exposure with a wind speed of ms − . The upper right panel shows the exposure also shown in Figure 1and has a wind speed of ms − . The lower left panel has a wind speed of ms − and the lower right panel has a wind speed of ms − . Four example results are shown in Figure 5 revealing a non-isotropic function. The example exposure shown in Figure 1 is thesame exposure presented in the upper right hand panel of Figure 5.A comparison of these two Figures shows that the elliptical ξ + ( | θ | ) has a central dipole with a major axis along the direction of the rip-ples seen in the δε i maps of Figure 1. For other exposures howeverwe see what appears to be a superposition of ripple patterns (see forexample the upper left and lower right panels of Figure 5), poten- c (cid:13) , 000–000 C. Heymans & B. Rowe et al.
Figure 6.
A comparison of wind direction and PSF residual direction as afunction of the ground wind speed measured at the start of each observationfor three samples grouped by exposure time with t s (stars), s > t s (circles) and t = 74 s (crosses). tially arising from different turbulent layers in the atmosphere. Inthese cases the orientation of the central dipole indicates the aver-age ripple direction. Based on these comparisons over the full dataset we use the orientation of the central dipole to determine an ef-fective ripple orientation and we measure this orientation using thequadrupole moments in equation 2) with r g = 2 . arcmin. Apply-ing this method to all exposures with t s, where the amplitudeof the atmospheric turbulence is sufficiently high, we can then com-pare this ripple direction to the wind direction relative to the image(shown by the arrow in the upper inset in Figure 5), calculated asdescribed in Appendix B.Figure 6 shows a compilation of results for three sets of data; t s (stars), s > t s (circles) and t = 74 s (crosses)comparing wind direction and PSF or ripple orientation. We findthat in general there is little evidence of strong correlation betweenwind direction and PSF orientation for the range of exposure timestested. This result is in disagreement with Asztalos et al. (2007)who found a relationship between wind direction and PSF ellip-ticities for wind speeds ranging from 2-6 m/s over 4 consecutivenights. We note that we could have drawn a similar conclusion hadwe analysed a smaller set of exposures, as the points that clusterin Figure 6 are typically taken on the same night. This clusteringwith observation date, if sets of observations were to be analysedin isolation, might well suggest a preferred wind-PSF orientationwhich we do not find when we analyse a large set of data spanningmany years of CFHT imaging. In this study we have analysed short exposure observations of densestellar fields to quantify the high spatial frequency variation of thePSF which we attribute to atmospheric effects. We have shown thatin short exposures with t s, the atmosphere contributes a sig-nificant correlated anisotropy to the PSF on angular scales θ < ′ that would dominate the cosmological signal measured in the lowerredshift bins of a cosmological tomographic analysis. On these an-gular scales the high spatial frequency of the atmospheric aberra-tion is too rapid to model with a typical stellar density and standardmethods. Figure 7.
The measured residual ellipticity correlation in our t = 10 s dataset compared to the expected WMAP7 cosmological two-point shear cor-relation signal for three tomographic bins with a mean redshift z = 0 . , z = 0 . and z = 0 . (solid lines, the lowest amplitude corresponds tothe lowest redshift bin). The data oscillates between positive and negativecorrelation which we indicate on this log-log plot using filled points wherethe data is negative which is primarily between the two dashed vertical lines. This does not mean, however, that multiple short exposure im-ages cannot not be combined in a way that reduces the cumula-tive impact of this (random) atmospheric anisotropy on estimatesof gravitational shear in the final analysis. It was found in thisstudy that the turbulent patterns in consecutive exposures are sta-tistically uncorrelated when separated by timescales s. Ourobservations motivate the combination of information from mul-tiple short exposure images in a way that takes advantage of thisfact. Examples of such approaches might be cross-correlating shearestimates from different exposures to estimate the shear autocor-relation, or combining images into a stacked, co-added image. Ifstacking, great care must taken to retain control and knowledge ofchanges in the PSF due to interpolation and stacking of dithereddata. Rowe, Hirata & Rhodes 2011 present an example of a linearimage stacking approach that seeks to preserve this information.As well as finding that consecutive exposures are uncorrelated,we have confirmed the predictions of de Vries et al. (2007) in find-ing that the amplitude of atmospheric distortion patterns decreaseswith time as t − / . Our results also show that this effect is not asignificant source of error for the CFHTLenS survey where for thefirst time a lensing survey is analysing the 600 second exposures in-dividually (Miller et al in prep), as compared to previous analysesusing a stack (Fu et al. 2008).Figure 7 compares the measured residual ellipticity correlationin our t = 10 s data set with the expected cosmological signal forthree tomographic bins with a mean redshift z = 0 . , z = 0 . , z = 0 . . This is the closest data set analysed to the proposed t = 15 s exposures for LSST. This figure shows that our findingsagree to some extent with Jee & Tyson (2010), who focus on an-gular scales θ > ′ and find that the PSF can be modeled to highaccuracy on these scales. For smaller angular scales, however atmo-spheric turbulence effects will be a significant source of systematicerror in a weak lensing analysis if not corrected for. Furthermore,as a result of the anisotropic nature of the residual correlation func-tion (shown in Figure 5) the effect of atmospheric turbulence mayleak to larger scales.For a future generation all-sky lensing surveyAmara & R´efr´egier (2007) derive a requirement on the vari- c (cid:13) , 000–000 tmospheric Distortions ance of systematic errors to be below σ < − such thatexperiments are limited by statistical noise rather than systematicerrors (see also van Waerbeke et al. (2006) for a similar conclu-sion). To reach this goal Paulin-Henriksson et al. (2008) show thatthis requires σ [ ε PSF ] . − for each ellipticity component and σ [ R ] . − R for ε PSF ≃ . and a typical galaxy/PSFsize ratio of R gal /R PSF & . . For the t = 10 s data we measure σ [ ε PSF ] = q h δε i i = (9 . ± . × − (9) σ [ R ] R = p h δR ih R i = (9 . ± . × − (10)where the error is calculated from the scatter between the 33 ex-posures in the t = 10 s data set showing the variance within theset, not the error on the mean. These measurements are an orderof magnitude larger in both size and ellipticity compared to the re-quirements quoted in Amara & R´efr´egier (2007).Using equation (15) in Paulin-Henriksson et al. (2008) we cancalculate the variance of the systematic errors we would expectfrom the residual atmospheric size and ellipticity variation that wemeasure from the t = 10 s data set: we find σ ∼ − . Thisvalue can be interpreted as the systematic variance produced inmultiple 10s exposures if: (i) the turbulent PSF cannot be mod-elled due to insufficient stellar density; and (ii) no attempt is madeto combine shear estimates from multiple exposures for which theturbulence patterns are mutually statistically independent. We notethat this calculation depends on the accuracy of both PSF elliptic-ity and PSF size characterisation, and both these are affected byatmospheric turbulence.We note that while the expectation of the turbulent contribu-tion to ε PSF across multiple exposures is zero, assuming any staticcomponent is removed, this is not the case for R . This has im-plications for plans to cross-correlate using only shear estimatesfrom different exposures to suppress PSF model errors. The utilityof such a scheme is that the ellipticity cross terms h δε i δε j i (for i = j ) have zero expectation value, and the variance of this prod-uct is the square of the individual variances. This product variancedecreases ∝ /N exp , where N exp is the total number of exposuresbeing cross-correlated, and therefore σ [ ε PSF ] decreases rapidly as ∝ /N exp in the cross-correlation.However, the same approach does not reduce uncertainty in R at the same rate because this quantity has an unknown ex-pectation value which must be characterized. Instead, σ [ R ] ∝ / p N exp . This differing behaviour in the combination of multipleexposures must be considered when comparing the results in equa-tions (9) & (10) to the required levels for the estimation of cosmicshear. The characterization of turbulent variation in R there-fore assumes added importance for upcoming surveys. Combinedwith N exp , the single-exposure value quoted in equation (10) canbe used to estimate this important effect for surveys using multi-ple exposures of duration t ≃ s, assuming similar atmosphericconditions to those at CFHT.Both the results in this paper and Wittman (2005) show an os-cillatory pattern in the residual ξ + correlation function at θ > . arcminutes, a pattern which is not consistent with the model pre-dicted by a simple, isotropic von K´arm´an spectrum (see Figure 3).Wittman (2005) describes this ringing as an artifact of the inter-polation scheme used to model the PSF on each individual chipwhich for the Suprime-Cam camera used in these observations hasa chip size of . ′ × . ′ . In this analysis we apply a similar,low-order chipwise PSF model, where for MegaCam the chip size is . ′ × . ′ . We thus agree with Wittman (2005) that one pos-sible origin for the ‘ringing’ is over-subtraction of a chipwise PSFmodel for the telescope optics, as the first dip occurs at a close tohalf the shorter dimension of the chip.However, it is less clear why chipwise fitting would cause anexcess at ′ . In addition, looking at the ξ + ( θ ) (i.e. not circularlyaveraged) two-point correlation functions (e.g. Figure 5) we seethat the oscillatory pattern is anisotropic in nature but often notaligned with the x - y chip grid. This alignment would be somethingto be expected if it were an artifact of the chipwise PSF model re-moval. We therefore offer an alternative suggestion that these rip-ples have an atmospheric origin, related to anisotropy in the atmo-spheric power spectrum (c.f. Figure 5), which cannot be generatedby simple isotropic models such as the von K´arm´an. This con-clusion is supported by the finding that the residual patterns areuncorrelated between exposures. However, without further data, oraccess to imaging with chips of differing dimensions, it remainsdifficult to say definitively which interpretation is correct.Although not performed explicitly in this study, it is worth-while mentioning the separation of the correlated turbulence sig-nal into ‘E’ and ‘B’ mode contributions (see, e.g., Crittenden et al.2002). The atmospheric PSF pattern for a turbulent sky may con-tain both E- and B-modes; if this PSF is not properly modelledany shear catalogue will then contain an imprint of this pattern. AnE/B-mode decomposition of the cosmic shear signal on the scaleswhere atmospheric turbulence patterns are correlated might there-fore detect this effect in the form of a non-zero B-mode. However,telescope optical abberrations, which also contribute strongly toPSF anisotropy across most angular scales of interest for lensing,tend to generate a more significant E-mode than B-mode (Hoekstra2004). In addition, cosmological effects can also generate B-modes(e.g., Schneider et al. 2002; Vale et al. 2004). The absence of strongB-modes in cosmic shear measurements therefore remains a prob-lematic indicator of the absence of systematic errors due to PSFanisotropy, and the potential presence of B-modes in turbulent PSFpatterns does not alter this fact.Finally, when considering our full sample of data, we havefound little evidence for a strong and consistent relationship be-tween the preferred direction of anisotropic correlations in atmo-spheric distortion and the ground-based wind velocity. We do notnecessarily expect to find such a relationship, as the all the tur-bulent layers responsible for the atmospheric PSF pattern are notnecessarily correlated with the ground conditions. However, we dofind that consecutive observations seem to yield similar angular off-sets between wind direction and turbulent PSF pattern orientation:such effects might lead to an over-interpretation of the correlationbetween wind direction and PSF pattern orientation in more limiteddata sets. In all, these results do not strongly motivate the use of theground-based wind direction and speed as an input to a PCA-typemodel of the atmospheric part of the PSF on short exposure data.However, it may still be useful in modelling out repeatable PSF dis-tortions due to telescope motion under wind (see, e.g., Jarvis & Jain2004). Once again, this motivates data analysis strategies that takeadvantage of the turbulent pattern being uncorrelated between suc-cessive exposures. We thank Chihway Chang, Phil Marshall and the rest of theCFHTLenS collaboration for convincing us to write this analysisinto a paper and for many useful discussions along the way. We c (cid:13) , 000–000 C. Heymans & B. Rowe et al. also thank the anonymous referee for helpful comments and theCFHT QSO team and PIs of the various data sets used in thisanalysis; Jerome Bouvier, Catherine Dougados, Eugene Magnier,Alan McConnachie and John Johnson. CH & BR acknowledge sup-port from the European Research Council under the EC FP7 grantnumbers 240185 (CH) & 240672 (BR). HH acknowledges supportfrom Marie Curie IRG grant 230924 and the Netherlands Organisa-tion for Scientific Research grant number 639.042.814. TE is sup-ported by the Deutsche Forschungsgemeinschaft through projectER 327/3-1 and the Transregional Collaborative Research CentreTR 33 - ”The Dark Universe. TDK was supported by a RAS 2010Fellowship. This study uses data from the CFHT Science DataArchive hosted and supported by the Canada Astronomy Data Cen-tre operated by the National Research Council of Canada with thesupport of the Canadian Space Agency.
Author contributions: All authors assisted with the development andwriting of this paper. CH and BR co-led the statistical analysis, HH con-ceived the project and led the data mining, catalogue production and dataanalysis, LM motivated the project and led the visualization of the effect.
APPENDIX A: CORRELATION FUNCTIONS FOR THEVON K ´ARM ´AN POWER SPECTRUM
As in the case of cosmological shear, the ellipticity correlationfunctions ξ + ( θ ) and ξ − ( θ ) (see Bartelmann & Schneider 2001)due to atmospheric effects can be related to power spectrum of thevon K´arm´an model as ξ + ( θ ) ∝ Z ∞ l d lJ (2 πlθ ) (cid:20) l + 1 θ (cid:21) − / (A1) ξ − ( θ ) ∝ Z ∞ l d lJ (2 πlθ ) (cid:20) l + 1 θ (cid:21) − / , (A2)where J ν ( x ) is the ν -th order Bessel function, and where we areassuming that the spatial power spectrum of ellipticity distortionsdue to turbulence follows P ( l ) as given in equation (7). Using theseexpressions, we calculate the ξ + correlation function as ξ + ( θ ) ∝ θ / K − / (2 πθ/θ ) , (A3)where K ν ( x ) is the modified Bessel function of the second kind(Arfken & Weber 2005). The expression for ξ − is more compli-cated: ξ − ( θ ) ∝ θ − / (cid:0) θ + π θ (cid:1) I − / (2 πθ/θ ) − π θ θ − / I − / (2 πθ/θ )+ θ − / (cid:18) θ − π θ (cid:19) I / (2 πθ/θ )+ 13 π θ θ − / I / (2 πθ/θ ) − θ / θ − (cid:18) θ o − π θ (cid:19) Γ( − / , (A4)where I ν ( x ) is the modified Bessel function of the first kind and Γ( x ) is the Gamma function (Arfken & Weber 2005). The func-tions I ν ( x ) diverge rapidly for x > , and the expression in equa-tion (A4) shows signs of numerical instability for θ & θ . In thisstudy we therefore make a crude approximation to avoid this be-haviour. We manually set ξ − ( θ ) = 0 where θ > . θ for the vonK´arm´an power spectrum model when investigating or plotting ξ − . APPENDIX B: THE PROJECTED DIRECTION OFGROUND WIND IN IMAGES ON THE CELESTIALSPHERE
We consider two general points with altitude-azimuth coordinates ( a , A ) and ( a , A ) in the horizontal coordinate system, whichuses the observer’s local horizon as the fundamental plane. As forany two points on the coordinate sphere, a single great circle may beconstructed that passes through ( a , A ) and ( a , A ) . We label as σ the angular size of the arc between these two points on the greatcircle. Using the spherical law of cosines, it is straightforward toshow that σ is given by cos ( σ ) = cos ( a ) cos ( a ) cos ( A − A )+ sin ( A ) sin ( A ) . (B1)This result will be used below to derive the projected direction ofground wind on telescope images.The pixel y -axis of the MegaCam instrument mounted onCFHT rotates to align with the celestial North Pole. In the hor-izontal coordinate system defined at the telescope, celestial Northhas the altitude-azimuth coordinates N = ( φ, where φ is the geo-graphic latitude of Mauna Kea. We define the telescope pointing di-rection (i.e. the direction of the image field centre) as T = ( a t , A t ) in this coordinate system. The points N and T on the celestialsphere can be connected by the arc of a great circle of constantRight Ascension. The y -axis of CFHT images will align with thisMeridian.The direction of ground wind is measured at CFHT and sup-plied with the astronomical images in degrees relative to compassNorth along the ground. In the horizontal coordinate system atCFHT this direction can therefore be expressed as W = (0 , A w ) ,where A w is the measured wind direction. We may then define asecond great circle passing through W and T . The angle θ at whichthis great circle intersects the Meridian through T and N is there-fore the projected direction of ground wind on the telescope image,relative to the CCD y -axis.The spherical law of cosines gives the following result for thisprojected wind direction angle θ : cos ( θ ) = cos ( γ ) − cos ( α ) cos ( β )sin ( α ) sin ( β ) , (B2)where from equation (B1) we have cos ( α ) = cos ( a t ) cos ( A w − A t ) , (B3) cos ( β ) = cos ( a t ) cos ( φ ) cos ( A t ) + sin ( a t ) sin ( φ ) , (B4) cos ( γ ) = cos ( A w ) cos ( φ ) . (B5)However, the expressions above do not provide the sign of θ , i.e.whether W lies to the West or East of the Meridian through T and N , the y -axis on the telescope image. This information is neededto correctly calculate the angular separation between the wind di-rection and primary PSF residual direction as described in Section3.3. The positive direction on the CCD x -axis of MegaCam canbe associated with positive hour angle H in the equatorial coordi-nate system. If we define the positive θ direction as that movingclockwise from the y -axis, then it is straightforwardly seen thatsgn ( θ ) = sgn ( H w − H t ) (B6)where sgn ( x ) is the signum function, H w is the hour angle of thewind direction vector W and H t is the hour angle of the telescopepointing direction T at the local sidereal time of observation. The c (cid:13) , 000–000 tmospheric Distortions telescope hour angle H t is provided by CFHT as metadata along-side all astronomical images, so it remains solely to determine H w .Using standard formulae for the conversion of horizontal toequatorial spherical coordinates, the declination δ w of the wind vec-tor W is given by sin ( δ w ) = cos ( φ ) cos ( A w ) . (B7)The hour angle of W is then given by H w = arcsin (cid:20) − sin ( A w )cos ( δ w ) (cid:21) = arcsin " − sin ( A w ) p − cos ( φ ) cos ( A w ) . (B8)When substituted into equation (B6) this allows the sign of θ to becalculated. In combination with the determination of the magnitudeof θ from equation (B2), the projected direction of ground wind ona telescope image pointing in the direction T is fully specified. REFERENCES
Amara A., R´efr´egier A., 2007, MNRAS, 381, 1018Arfken G. B., Weber H. J., 2005, Mathematical methods for physi-cists 6th ed.Asztalos S., de Vries W. H., Rosenberg L. J., Treadway T., BurkeD., Claver C., Saha A., Puxley P., 2007, ApJ, 659, 69Bartelmann M., Schneider P., 2001, Physics Reports, 340, 291Crittenden R., Natarajan R., Pen U., Theuns T., 2002, ApJ, 568,20de Vries W. H., Olivier S. S., Asztalos S. J., Rosenberg L. J., BakerK. L., 2007, ApJ, 662, 744Fu L., Semboloni E., Hoekstra H., Kilbinger M., van Waerbeke L.,Tereno I., Mellier Y., Heymans C., Coupon J., 10 other authors2008, A&A, 479, 9Hoekstra H., 2004, MNRAS, 347, 1337Hoekstra H., Mellier Y., van Waerbeke L., Semboloni E., Fu L.,Hudson M. J., Parker L. C., Tereno I., Benabed K., 2006, ApJ,647, 116Huterer D., 2010, General Relativity and Gravitation, 42, 2177Jarvis M., Jain B., 2004, ApJ submitted, astroph/0412234Jee M. J., Tyson J. A., 2010, ArXiv e-printsKaiser N., Burgett W., Chambers K., Denneau L., Heasley J.,Jedicke R., Magnier E., Morgan J., Onaka P., Tonry J., 2010, inSociety of Photo-Optical Instrumentation Engineers (SPIE) Con-ference Series Vol. 7733 of Society of Photo-Optical Instrumen-tation Engineers (SPIE) Conference Series, The Pan-STARRSwide-field optical/NIR imaging surveyKaiser N., Squires G., Broadhurst T., 1995, ApJ, 449, 460Kitching T., Amara A., Gill M., Harmeling S., Heymans C.,Massey R., Rowe B., Schrabback T., Voigt L., GREAT10 Con-sortium ., 2010, ArXiv e-printsKomatsu E., Smith K. M., Dunkley J., Bennett C. L., WMAP Col-laboration 2011, ApJS, 192, 18Kostroun V. O., 1980, Nuclear Instruments and Methods, 172, 371Laureijs et al. R., 2011, Euclid Definition Study ReportLSST Science Collaboration 2009, ArXiv e-printsPaulin-Henriksson S., Amara A., Voigt L., Refregier A., BridleS. L., 2008, A&A, 484, 67Press W., Teukolsky S., Vetterling W., Flannery B., 1986, Numer-ical recipes. Cambridge University PressRowe B., 2010, MNRAS, 404, 350 Rowe B., Hirata C., Rhodes J., 2011, ApJ, 741, 46Sasiela R. J., 1994, Electromagnetic wave propagation in turbu-lence. Evaluation and application of Mellin transformsSchneider P., Van Waerbeke L., Mellier Y., 2002, A&A, 389, 741Vale C., Hoekstra H., van Waerbeke L., White M., 2004, ApJL,613, L1van Waerbeke L., White M., Hoekstra H., Heymans C., 2006, As-troparticle Physics, 26, 91Wittman D., 2005, ApJL, 632, L5 c (cid:13)000