Atmospheric turbulence profiling with multi-aperture scintillation of a Shack-Hartmann sensor
MMNRAS , 1–11 (2020) Preprint 13 January 2021 Compiled using MNRAS L A TEX style file v3.0
Atmospheric turbulence profiling with multi-aperture scintillation of aShack-Hartmann sensor
Hajime Ogane, ★ Masayuki Akiyama, Shin Oya and Yoshito Ono Astronomical Institute, Tohoku University, 6-3 Aramaki, Aoba-ku Sendai, Miyagi 980-8578, Japan National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan Subaru Telescope, National Astronomical Observatory of Japan, 650 North Aohoku Place Hilo, HI 96720, USA
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Adaptive optics (AO) systems using tomographic estimation of three-dimensional structure of atmospheric turbulence requiresvertical atmospheric turbulence profile, which describes turbulence strength as a function of altitude as a prior information. Wepropose a novel method to reconstruct the profile by applying Multi Aperture Scintillation Sensor (MASS) method to scintillationdata obtained by a Shack-Hartmann wavefront sensor (SH-WFS). Compared to the traditional MASS, which uses atmosphericscintillation within 4 concentric annular apertures, the new method utilizes scintillation in several hundreds of spatial patterns,which are created by combinations of SH-WFS subapertures. Accuracy of the turbulence profile reconstruction is evaluated withBayesian inference, and it is confirmed that turbulence profile with more than 10 layers can be reconstructed thanks to the largenumber of constraints. We demonstrate the new method with a SH-WFS attached to the 50 cm telescope at Tohoku universityand confirm that general characteristics of atmospheric turbulence profile is reproduced.
Key words: atmospheric turbulence – scintillation – Shack-Hartmann wavefront sensor – adaptive optics
Refractive index fluctuation in the Earth’s atmosphere distorts thewavefront or equiphase surface of the starlight and causes blur ofthe stellar image. Adaptive optics (AO) systems realize diffraction-limited spatial resolution images with ground-based large aperturetelescopes. In AO systems, with using a natural guide star (NGS) oran artificial laser guide star (LGS) as a reference source, the distortionof wavefront is measured by a wavefront sensor (WFS) and correctedby a deformable mirror (DM) in a timescale of ∼ ¥ ‘ere et al. 2014, Rigaut et al. 2014, Minowaet al. 2017). These systems measure the wavefront distortion in sev-eral lines of sight and reconstruct the distortions optimized in thedirection of science objects using tomographic estimation, aimingAO correction for lower wavefront error (Laser Tomography Adap-tive Optics, LTAO) or wider field of view (Multi Conjugate AdaptiveOptics, MCAO, Beckers 1988, Rigaut & Neichel 2018; Ground LayerAdaptive Optics, GLAO, Rigaut 2002, Tokovinin 2004; Multi ObjectAdaptive Optics, MOAO, Hammer et al. 2004, Vidal et al. 2010).The tomographic turbulence estimation is essential technique fornext-generation giant segmented mirror telescopes (GSMTs), whichhave 30m-class primary mirrors.The tomographic estimation of the three-dimensional turbulence ★ E-mail: [email protected] structure requires prior information of the strength of atmosphericturbulence as a function of altitude, which is called atmosphericturbulence profile. Tomographic reconstruction matrix is computedfrom the positions of guide stars and vertical atmospheric turbulenceprofile. Imperfect prior information of the turbulence profile causestomographic error, which accounts for a large fraction of the totalAO error budget (Gilles et al. 2008). Because atmospheric turbulenceprofile varies with time, tomographic reconstruction matrix shouldbe updated in a timescale of tens of minutes, which corresponds tothe typical time scale of the profile time evolution (Gendron et al.2014, Farley et al. 2020).Altitude resolution is an important parameter for the atmosphericturbulence profiling to reflect the precise turbulence distribution inthe tomographic reconstruction matrix. Fusco & Costille (2010) andCostille & Fusco (2012) studied the effect of the number of lay-ers in turbulence profile on the tomographic error in ELT-scaletomographic AO systems. They created less resolved turbulenceprofile by under-sampling the original highly-resolved (250 layers, Δ ℎ ∼
100 m) turbulence profile obtained by a balloon experimentand investigated the impact of the number of layers for reconstruc-tion on the tomographic error. The conclusion is that at least 10-20layers are needed to achieve the tomographic error comparable tothat obtained using 250-layer profile. The required number of layersdepends on LGS asterism diameter and tolerated tomographic errorand can be reduced by optimizing the altitude combination of theprofile or optimized compression method (Saxenhuber et al. 2017,Farley et al. 2020).A number of methods to obtain real-time atmospheric turbulenceprofile based on optical triangulation have been developed. Scintilla-tion detection and ranging (SCIDAR; Rocca et al. 1974) uses spatial © a r X i v : . [ a s t r o - ph . I M ] J a n H.Ogane et al. correlation of scintillation map on the pupil plane from two brightstars with an angular separation of several arcseconds. Because vari-ance of intensity fluctuation is proportional to apparent altitude withthe power of 5/6, SCIDAR does not have sensitivity to turbulence atlow altitudes. Generalized-SCIDAR (G-SCIDAR, Avila et al. 1997)has overcome this limitation by detecting scintillation on a plane atsome distance away from the pupil plane. Slope detection and rang-ing (SLODAR; Wilson 2002) has the same triangulation principlewith SCIDAR but it uses spatial correlation of phase map on thepupil plane, and the correlation have sensitivity to the ground layer.Because these methods are based on optical triangulation betweentwo stars, these methods can be limited by the availability of doublestars. In addition, these methods do not have any sensitivity to tur-bulence at high altitudes since spatial correlation length created byhigh turbulence layer is larger than the size of the pupil.Multi aperture scintillation sensor and differential image motionmonitor (MASS-DIMM; Kornilov et al. 2007) is one of the mostcommon profilers, which uses a single star and has lower altituderesolution compared to SLODAR. This is a combined method ofMASS (Tokovinin et al. 2003, Kornilov et al. 2003) and DIMM(Sarazin & Roddier 1990). MASS reconstructs 6-layer profile offree-atmosphere based on the scintillation from a single bright stardetected by some spatial patterns on the pupil plane. On the otherhand, DIMM measures total seeing or Fried parameter based on thevariance of differential image motion through two small apertures.Ground-layer seeing is estimated with the difference between theMASS seeing (free-atmosphere seeing) and the DIMM seeing (totalseeing). Then, 7-layer (ground-layer plus 6 MASS layers) turbulenceprofile can be obtained.Recently, in order to improve the altitude resolution of MASS,fine spatial sampling of scintillation using fast and low-noise de-tector has been demonstrated by Full Aperture Scintillation Sensor(FASS; Guesalaga et al. 2016, Guesalaga et al. 2020) method. FASSmeasures angular power spectrum of scintillation on the pupil planeand compare it with simulated spectra to reconstruct a turbulenceprofile with 14 layers from 0.3 km to 25 km above the ground.In this paper, we propose another new turbulence profiling method,which carries out scintillation measurement similar to MASS usingspot brightness fluctuation data of Shack-Hartmann wavefront sen-sor (SH-WFS). Hereafter we call this new method SH-MASS. Com-pared to the traditional MASS instrument, which measures scintilla-tion with 4 concentric annuli, SH-WFS can measure the scintillationwith many spacial patterns created from combinations of subaper-tures of SH-WFS. Thanks to the larger number of measurements,SH-MASS approach utilizes more constraints on the turbulence pro-file than MASS. These constraints make it possible to estimate theatmospheric turbulence profile with high altitude resolution by theobservation of scintillation of a single star. As an additional meritof using SH-WFS, it is worth noting that SH-WFS can be combinedwith above-mentioned slope-based profiling methods. For example,DIMM can be conducted using differential image motions of twoSH-WFS spots and SLODAR can also be executed simultaneously iftwo SH-WFSs are available.This paper is organized as follows. In section 2, the principle ofMASS is reviewed and application of the method to SH-WFS data isexplained. In section 3, the response function of SH-MASS is evalu-ated through simulations with a single layer turbulence. In section 4,on-sky experiment conducted with 50cm telescope at Tohoku univer-sity is explained and the results of atmospheric profiling are shown.In section 5, remaining problems and limitations of SH-MASS arediscussed. Finally, the whole contents are summarized in section 6. Spatial frequency [ m ]10 S c i n t ill a t i o n p o w e r [ m ]
20 km propagation10 km propagation1 km propagation
Figure 1.
One-dimensional spatial power spectrums of scintillationare shown. Vertical axis means Φ 𝐼 ( 𝑓 𝑥 , 𝑓 𝑦 ) integrated at constant 𝑓 ;2 𝜋 𝑓 Φ 𝐼 ( 𝑓 𝑥 , 𝑓 𝑦 ) . Monochromatic ( 𝜆 = 𝐶 𝑁 ( ℎ ) Δ ℎ = . 𝐸 − [ m / ] ) at the propagation distance of1km, 10km, and 20km is shown as blue solid, orange dashed, and green dottedline, respectively. In this subsection, the principle of MASS is summarized followingKornilov et al. (2003) and Tokovinin et al. (2003). Let us considerthat light with a wavelength of 𝜆 passes through a single layer ofatmospheric turbulence whose altitude is ℎ and thickness is Δ ℎ .The phase of the light is distorted by refractive index fluctuationin the turbulence layer according to the turbulence strength which isindicated by structure constant of refractive index fluctuation 𝐶 𝑁 ( ℎ ) .Then, assuming the Kolmogorov’s atmospheric turbulence model,the spatial power spectrum of the phase fluctuation Φ 𝜙 [ m ] is writtenas follows, Φ 𝜙 ( 𝑓 𝑥 , 𝑓 𝑦 ) = . 𝜆 − 𝑓 − / 𝐶 𝑁 ( ℎ ) Δ ℎ, (1)where 𝑓 𝑥 , 𝑓 𝑦 are spatial frequencies and 𝑓 = √︃ 𝑓 𝑥 + 𝑓 𝑦 . Besides,by assuming the Fresnel propagation, and the weak perturbation(Rytov’s) approximation, which is applicable for astronomical obser-vation at not very large zenith angle 𝑧 , the spatial power spectrum ofthe intensity fluctuation Φ 𝐼 [ m ] is written as follows, Φ 𝐼 ( 𝑓 𝑥 , 𝑓 𝑦 ) = (cid:16) 𝜋𝜆ℎ sec ( 𝑧 ) 𝑓 (cid:17) Φ 𝜙 ( 𝑓 𝑥 , 𝑓 𝑦 ) . (2)Here, ℎ sec ( 𝑧 ) represents apparent altitude of the turbulence layer,or light propagation distance from the layer. Generally, multiple tur-bulence layers at different altitude and thickness affect the intensityfluctuation. Total contribution from these multiple layers can be writ-ten as linear combination of that from each layer, Φ 𝐼 ( 𝑓 𝑥 , 𝑓 𝑦 ) = 𝑁 layer ∑︁ 𝑖 . 𝑓 − / (cid:26) sin ( 𝜋𝜆ℎ 𝑖 sec ( 𝑧 ) 𝑓 ) 𝜆 (cid:27) 𝐶 𝑁 ( ℎ 𝑖 ) Δ ℎ 𝑖 , (3)where 𝑁 layer is the number of layers and 𝑖 is the index of eachturbulence layer.Fig.1 shows the power spectrum of intensity fluctuation. In thisplot, Φ 𝐼 ( 𝑓 𝑥 , 𝑓 𝑦 ) in Eq.2 integrated at constant 𝑓 , i.e. 2 𝜋 𝑓 Φ 𝐼 ( 𝑓 𝑥 , 𝑓 𝑦 ) MNRAS , 1–11 (2020) tmospheric turbulence profiling with a SH-WFS is shown as a function of 𝑓 . Single turbulence layer with a constantturbulence strength and monochromatic ( 𝜆 =
500 nm) observationare assumed. The spatial frequency at which the power of scintillationhas its peak depends on the propagation distance of the turbulencelayer. The layer is higher, the lower spatial frequency accounts forlarge amount of the power. Hence, detecting scintillation at differentspatial frequencies makes it possible to discern the contributionsfrom turbulence layers at different altitudes. The amplitude of thepower of scintillation also depends on the altitude of the turbulencelayer, and contributions from lower layers are weaker. Scintillationhardly has the information of turbulence at ground layer or domeseeing because Φ 𝐼 in Eq.2 becomes 0 at ℎ = 𝐼 𝑋 , the intensity variance of the X-th annulus, referred asnormal scintillation index, is defined as follows, 𝑠 𝑋 = Var (cid:20) 𝐼 𝑋 (cid:104) 𝐼 𝑋 (cid:105) (cid:21) , (4)where (cid:104)(cid:105) represents time-average, and Var means variance. Likewise,the intensity covariance of the X-th and Y-th annuli, referred as adifferential scintillation index, is defined as 𝑠 𝑋𝑌 = Var (cid:20) 𝐼 𝑋 (cid:104) 𝐼 𝑋 (cid:105) − 𝐼 𝑌 (cid:104) 𝐼 𝑌 (cid:105) (cid:21) = 𝑠 𝑋 + 𝑠 𝑌 − (cid:20) 𝐼 𝑋 (cid:104) 𝐼 𝑋 (cid:105) , 𝐼 𝑌 (cid:104) 𝐼 𝑌 (cid:105) (cid:21) , (5)where Cov means covariance.These SIs can be expressed using the power spectrum of intensityfluctuation Φ 𝐼 ( 𝑓 𝑥 , 𝑓 𝑦 ) as follows, 𝑠 𝑋 = ∬ Φ 𝐼 ( 𝑓 𝑥 , 𝑓 𝑦 )|F [ 𝐴 𝑋 ( 𝑥, 𝑦 )]| 𝑑𝑓 𝑥 𝑑𝑓 𝑦 , (6) 𝑠 𝑋𝑌 = ∬ Φ 𝐼 ( 𝑓 𝑥 , 𝑓 𝑦 )|F [ 𝐴 𝑋 ( 𝑥, 𝑦 ) − 𝐴 𝑌 ( 𝑥, 𝑦 )]| 𝑑𝑓 𝑥 𝑑𝑓 𝑦 , (7)where F means Fourier transformation and 𝐴 ( 𝑥, 𝑦 ) is the normalizedaperture function; a function which returns value of 1 divided by thearea of the aperture for ( 𝑥, 𝑦 ) inside the aperture and value of 0for others. By substituting Eq.3 into Φ 𝐼 ( 𝑓 𝑥 , 𝑓 𝑦 ) in Eq.6 and Eq.7,following formula are obtained. 𝑠 𝑋 = 𝑁 layer ∑︁ 𝑖 𝑊 𝑋,𝑖 𝐽 𝑖 , (8) 𝑠 𝑋𝑌 = 𝑁 layer ∑︁ 𝑖 𝑊 𝑋𝑌 ,𝑖 𝐽 𝑖 , (9)where 𝑊 𝑋,𝑖 = ∬ . 𝑓 − / (cid:26) sin ( 𝜋𝜆ℎ 𝑖 sec ( 𝑧 ) 𝑓 ) 𝜆 (cid:27) |F [ 𝐴 𝑋 ( 𝑥, 𝑦 )]| 𝑑𝑓 𝑥 𝑑𝑓 𝑦 , (10) 𝑊 𝑋𝑌 ,𝑖 = ∬ . 𝑓 − / (cid:26) sin ( 𝜋𝜆ℎ 𝑖 sec ( 𝑧 ) 𝑓 ) 𝜆 (cid:27) |F [ 𝐴 𝑋 ( 𝑥, 𝑦 ) − 𝐴 𝑌 ( 𝑥, 𝑦 )]| 𝑑𝑓 𝑥 𝑑𝑓 𝑦 . (11) 𝐽 𝑖 = 𝐶 𝑁 ( ℎ 𝑖 ) Δ ℎ 𝑖 . (12) 𝑊 𝑋,𝑖 and 𝑊 𝑋𝑌 ,𝑖 are called normal weighting functions (WFs) anddifferential WFs, respectively, which can be calculated from the in-formation of aperture geometry and measurement wavelength. Bysolving Eq.8 and Eq.9, turbulence strengths 𝐽 𝑖 = 𝐶 𝑁 ( ℎ 𝑖 ) Δ ℎ 𝑖 ofmultiple layers are estimated. Typical MASS instrument has fourconcentric annular apertures whose diameters are 2.0, 3.7, 7.0, and13.0 cm, respectively. Then, ten SIs (four normal SIs plus six differ-ential SIs) are derived in order to reconstruct a turbulence profile of6 turbulence layers (0.5, 1, 2, 4, 8, and 16 km above the aperture).Each SI is calculated every minute from the photon counting with ∼ 𝑊 𝑋,𝑖 = ∬ . 𝑓 − / (cid:26)∫ sin ( 𝜋𝜆ℎ 𝑖 sec ( 𝑧 ) 𝑓 ) 𝜆 𝐹 ( 𝜆 ) 𝑑𝜆 (cid:27) |F [ 𝐴 𝑋 ( 𝑥, 𝑦 )]| 𝑑𝑓 𝑥 𝑑𝑓 𝑦 , (13) 𝑊 𝑋𝑌 ,𝑖 = ∬ . 𝑓 − / (cid:26)∫ sin ( 𝜋𝜆ℎ 𝑖 sec ( 𝑧 ) 𝑓 ) 𝜆 𝐹 ( 𝜆 ) 𝑑𝜆 (cid:27) |F [ 𝐴 𝑋 ( 𝑥, 𝑦 ) − 𝐴 𝑌 ( 𝑥, 𝑦 )]| 𝑑𝑓 𝑥 𝑑𝑓 𝑦 , (14)where 𝐹 ( 𝜆 ) is normalized spectral function which contains all thespectral characteristics mentioned above and is normalized to satisfy ∫ ∞−∞ 𝐹 ( 𝜆 ) 𝑑𝜆 = Because SH-WFS effectively divides the entrance pupil into gridpattern subapertures, it can be used to measure scintillation in manyspatial patterns by multiple combinations of subapertures. Therefore,SH-WFS is applicable for conducting MASS.Fig.2 shows comparison of the definition of spatial patterns intraditional MASS (top two panels) and SH-MASS (bottom two pan-els). In traditional MASS case, concentric annular spatial patternsare used in order to extract scintillation at a specific frequency whichcorresponds to the diameters of the annulus. For example, a normalSI is defined as a intensity variance measured by the red annulusin the left top panel of Fig.2 while a differential SI is defined as aintensity covariance between the red and blue annulus in the righttop panel of Fig.2.On the other hand, we define a subaperture pair, which consistsof two subapertures, as one spatial pattern of SH-MASS so thatwe can effectively extract a certain spatial frequency componentof scintillation which is characterized by the distance of the twosubapertures. Then, total intensity of a subaperture pair is used asa measured value to calculate a SI. For example, a normal SI isdefined as a fluctuation variance of total intensity measured by thetwo red subapertures in the left bottom panel of Fig.2. Based on thisdefinition of normal SI, Eq.4 can be rewritten as, 𝑠 𝑋 = Var (cid:20) 𝐼 𝑋 (cid:104) 𝐼 𝑋 (cid:105) (cid:21) = Var (cid:20) 𝐼 𝑖 + 𝐼 𝑗 (cid:104) 𝐼 𝑖 + 𝐼 𝑗 (cid:105) (cid:21) = Var [ 𝐼 𝑖 ] + Var [ 𝐼 𝑗 ] + [ 𝐼 𝑖 , 𝐼 𝑗 ]((cid:104) 𝐼 𝑖 (cid:105) + (cid:104) 𝐼 𝑗 (cid:105)) , (15)where 𝑖, 𝑗 are indices of subapertures which constitutes aperture 𝑋 , MNRAS000
500 nm) observationare assumed. The spatial frequency at which the power of scintillationhas its peak depends on the propagation distance of the turbulencelayer. The layer is higher, the lower spatial frequency accounts forlarge amount of the power. Hence, detecting scintillation at differentspatial frequencies makes it possible to discern the contributionsfrom turbulence layers at different altitudes. The amplitude of thepower of scintillation also depends on the altitude of the turbulencelayer, and contributions from lower layers are weaker. Scintillationhardly has the information of turbulence at ground layer or domeseeing because Φ 𝐼 in Eq.2 becomes 0 at ℎ = 𝐼 𝑋 , the intensity variance of the X-th annulus, referred asnormal scintillation index, is defined as follows, 𝑠 𝑋 = Var (cid:20) 𝐼 𝑋 (cid:104) 𝐼 𝑋 (cid:105) (cid:21) , (4)where (cid:104)(cid:105) represents time-average, and Var means variance. Likewise,the intensity covariance of the X-th and Y-th annuli, referred as adifferential scintillation index, is defined as 𝑠 𝑋𝑌 = Var (cid:20) 𝐼 𝑋 (cid:104) 𝐼 𝑋 (cid:105) − 𝐼 𝑌 (cid:104) 𝐼 𝑌 (cid:105) (cid:21) = 𝑠 𝑋 + 𝑠 𝑌 − (cid:20) 𝐼 𝑋 (cid:104) 𝐼 𝑋 (cid:105) , 𝐼 𝑌 (cid:104) 𝐼 𝑌 (cid:105) (cid:21) , (5)where Cov means covariance.These SIs can be expressed using the power spectrum of intensityfluctuation Φ 𝐼 ( 𝑓 𝑥 , 𝑓 𝑦 ) as follows, 𝑠 𝑋 = ∬ Φ 𝐼 ( 𝑓 𝑥 , 𝑓 𝑦 )|F [ 𝐴 𝑋 ( 𝑥, 𝑦 )]| 𝑑𝑓 𝑥 𝑑𝑓 𝑦 , (6) 𝑠 𝑋𝑌 = ∬ Φ 𝐼 ( 𝑓 𝑥 , 𝑓 𝑦 )|F [ 𝐴 𝑋 ( 𝑥, 𝑦 ) − 𝐴 𝑌 ( 𝑥, 𝑦 )]| 𝑑𝑓 𝑥 𝑑𝑓 𝑦 , (7)where F means Fourier transformation and 𝐴 ( 𝑥, 𝑦 ) is the normalizedaperture function; a function which returns value of 1 divided by thearea of the aperture for ( 𝑥, 𝑦 ) inside the aperture and value of 0for others. By substituting Eq.3 into Φ 𝐼 ( 𝑓 𝑥 , 𝑓 𝑦 ) in Eq.6 and Eq.7,following formula are obtained. 𝑠 𝑋 = 𝑁 layer ∑︁ 𝑖 𝑊 𝑋,𝑖 𝐽 𝑖 , (8) 𝑠 𝑋𝑌 = 𝑁 layer ∑︁ 𝑖 𝑊 𝑋𝑌 ,𝑖 𝐽 𝑖 , (9)where 𝑊 𝑋,𝑖 = ∬ . 𝑓 − / (cid:26) sin ( 𝜋𝜆ℎ 𝑖 sec ( 𝑧 ) 𝑓 ) 𝜆 (cid:27) |F [ 𝐴 𝑋 ( 𝑥, 𝑦 )]| 𝑑𝑓 𝑥 𝑑𝑓 𝑦 , (10) 𝑊 𝑋𝑌 ,𝑖 = ∬ . 𝑓 − / (cid:26) sin ( 𝜋𝜆ℎ 𝑖 sec ( 𝑧 ) 𝑓 ) 𝜆 (cid:27) |F [ 𝐴 𝑋 ( 𝑥, 𝑦 ) − 𝐴 𝑌 ( 𝑥, 𝑦 )]| 𝑑𝑓 𝑥 𝑑𝑓 𝑦 . (11) 𝐽 𝑖 = 𝐶 𝑁 ( ℎ 𝑖 ) Δ ℎ 𝑖 . (12) 𝑊 𝑋,𝑖 and 𝑊 𝑋𝑌 ,𝑖 are called normal weighting functions (WFs) anddifferential WFs, respectively, which can be calculated from the in-formation of aperture geometry and measurement wavelength. Bysolving Eq.8 and Eq.9, turbulence strengths 𝐽 𝑖 = 𝐶 𝑁 ( ℎ 𝑖 ) Δ ℎ 𝑖 ofmultiple layers are estimated. Typical MASS instrument has fourconcentric annular apertures whose diameters are 2.0, 3.7, 7.0, and13.0 cm, respectively. Then, ten SIs (four normal SIs plus six differ-ential SIs) are derived in order to reconstruct a turbulence profile of6 turbulence layers (0.5, 1, 2, 4, 8, and 16 km above the aperture).Each SI is calculated every minute from the photon counting with ∼ 𝑊 𝑋,𝑖 = ∬ . 𝑓 − / (cid:26)∫ sin ( 𝜋𝜆ℎ 𝑖 sec ( 𝑧 ) 𝑓 ) 𝜆 𝐹 ( 𝜆 ) 𝑑𝜆 (cid:27) |F [ 𝐴 𝑋 ( 𝑥, 𝑦 )]| 𝑑𝑓 𝑥 𝑑𝑓 𝑦 , (13) 𝑊 𝑋𝑌 ,𝑖 = ∬ . 𝑓 − / (cid:26)∫ sin ( 𝜋𝜆ℎ 𝑖 sec ( 𝑧 ) 𝑓 ) 𝜆 𝐹 ( 𝜆 ) 𝑑𝜆 (cid:27) |F [ 𝐴 𝑋 ( 𝑥, 𝑦 ) − 𝐴 𝑌 ( 𝑥, 𝑦 )]| 𝑑𝑓 𝑥 𝑑𝑓 𝑦 , (14)where 𝐹 ( 𝜆 ) is normalized spectral function which contains all thespectral characteristics mentioned above and is normalized to satisfy ∫ ∞−∞ 𝐹 ( 𝜆 ) 𝑑𝜆 = Because SH-WFS effectively divides the entrance pupil into gridpattern subapertures, it can be used to measure scintillation in manyspatial patterns by multiple combinations of subapertures. Therefore,SH-WFS is applicable for conducting MASS.Fig.2 shows comparison of the definition of spatial patterns intraditional MASS (top two panels) and SH-MASS (bottom two pan-els). In traditional MASS case, concentric annular spatial patternsare used in order to extract scintillation at a specific frequency whichcorresponds to the diameters of the annulus. For example, a normalSI is defined as a intensity variance measured by the red annulusin the left top panel of Fig.2 while a differential SI is defined as aintensity covariance between the red and blue annulus in the righttop panel of Fig.2.On the other hand, we define a subaperture pair, which consistsof two subapertures, as one spatial pattern of SH-MASS so thatwe can effectively extract a certain spatial frequency componentof scintillation which is characterized by the distance of the twosubapertures. Then, total intensity of a subaperture pair is used asa measured value to calculate a SI. For example, a normal SI isdefined as a fluctuation variance of total intensity measured by thetwo red subapertures in the left bottom panel of Fig.2. Based on thisdefinition of normal SI, Eq.4 can be rewritten as, 𝑠 𝑋 = Var (cid:20) 𝐼 𝑋 (cid:104) 𝐼 𝑋 (cid:105) (cid:21) = Var (cid:20) 𝐼 𝑖 + 𝐼 𝑗 (cid:104) 𝐼 𝑖 + 𝐼 𝑗 (cid:105) (cid:21) = Var [ 𝐼 𝑖 ] + Var [ 𝐼 𝑗 ] + [ 𝐼 𝑖 , 𝐼 𝑗 ]((cid:104) 𝐼 𝑖 (cid:105) + (cid:104) 𝐼 𝑗 (cid:105)) , (15)where 𝑖, 𝑗 are indices of subapertures which constitutes aperture 𝑋 , MNRAS000 , 1–11 (2020)
H.Ogane et al.
Normal aperture pattern of MASSNormal aperture pattern of SH-MASS Differential aperture pattern of MASSDifferential aperture pattern of SH-MASS c m on a p e r t u r e p l a n e Figure 2.
In the top panels, examples of MASS spatial pattern whose mea-surements are used to compute SIs are shown. For example, a normal SI iscalculated as the intensity variance measured by red annular aperture while adifferential SI is calculated as the intensity covariance between the measure-ment by red aperture and that by blue aperture. In bottom panels, examplesof spatial pattern for SH-MASS are shown. In our definition, a normal SI iscalculated as the intensity variance observed by red subapeture pair while adifferential SI is computed as the intensity covariance between the observa-tion by red subaperture pair and blue subaperture pair. In the calculation ofdifferential SIs, only two subaperture pairs which have common mid pointare used. and 𝐼 𝑖 represents spot intensity (or counts) observed in i-th subaper-ture. Thanks to a large number of SH-WFS subapertures, there aremany subaperture pairs which have common separation distance.Then, we calculated normal SIs for all subaperture pairs which havea common spatial distance and regarded average and standard devi-ation as a normal SI and its measurement error, respectively.Likewise, a differential SI is defined as a fluctuation covariancebetween total intensity measured by the two red subapertures andthat measured by the two blue subapertures in the right bottom panelof Fig.2. Then, Eq.5 can be rewritten as, 𝑠 𝑋𝑌 = 𝑠 𝑋 + 𝑠 𝑌 − (cid:20) 𝐼 𝑋 (cid:104) 𝐼 𝑋 (cid:105) , 𝐼 𝑌 (cid:104) 𝐼 𝑌 (cid:105) (cid:21) = 𝑠 𝑋 + 𝑠 𝑌 − (cid:20) 𝐼 𝑖 + 𝐼 𝑗 (cid:104) 𝐼 𝑖 + 𝐼 𝑗 (cid:105) , 𝐼 𝑘 + 𝐼 𝑙 (cid:104) 𝐼 𝑘 + 𝐼 𝑙 (cid:105) (cid:21) = 𝑠 𝑋 + 𝑠 𝑌 − [ 𝐼 𝑖 , 𝐼 𝑘 ] + Cov [ 𝐼 𝑗 , 𝐼 𝑘 ] + Cov [ 𝐼 𝑖 , 𝐼 𝑙 ] + Cov [ 𝐼 𝑗 , 𝐼 𝑙 ]((cid:104) 𝐼 𝑖 (cid:105) + (cid:104) 𝐼 𝑗 (cid:105))((cid:104) 𝐼 𝑘 (cid:105) + (cid:104) 𝐼 𝑙 (cid:105)) (16)where 𝑖, 𝑗 are indices of subapertures which constitutes aperture 𝑋 , while 𝑘, 𝑙 means indices of subapertures for aperture 𝑌 . Here,we calculated differential SIs for two subaperture pairs which havecommon mid point. That corresponds to taking concentric two annuliin the traditional MASS. By these definitions of spatial patterns inSH-MASS, 51 normal SIs and 234 differential SIs are obtained with10 ×
10 SH-WFS.Fig.3 shows the comparison of WFs of the traditional MASS andSH-MASS which are calculated assuming the aperture geometriesshown in Fig.2 and measurement wavelength of 500 nm. Each rowof the WF matrix represents WF of each spatial pattern, i.e. Eq.10or Eq.11. Here, values of WFs are normalized in each row so thatthe weight value transition in the direction of propagation distance N o r m a l D i ff e r e n ti a l N o r m a l D i ff e r e n ti a l Figure 3.
Top: WF matrix of the traditional MASS. Ten spatial patterns aremade by four concentric annular apertures whose diameters are 2.0, 3.7, 7.0,and 13.0 cm, respectively. Bottom: WF matrix of the SH-MASS. As the ge-ometry of SH-WFS, 10 ×
10 subapertures with each subaperture size scales1 . [ m − / ] normalized in each row so that characteristic prop-agation distance where weight value reaches 50%-ile can be seen as lightgreen. can be easily recognized. In traditional MASS case, the number ofspatial patterns is only ten (four normal plus six differential) and thetransition occurs at longer propagation distance for normal spatialpatterns while at shorter propagation distance for differential spatialpatterns. While the traditional MASS WF has a small number ofspatial patterns with a discontinuity at ∼ MNRAS , 1–11 (2020) tmospheric turbulence profiling with a SH-WFS the number of spatial patterns reaches ∼
300 and their transition dis-tances are continuous from the ground to 20 km high, which impliesSH-MASS’s aperture geometry gives sufficient number of constraintsto estimate a turbulence profile with higher altitude resolution.
Reconstruction of a turbulence profile is solving an inverse problemdescribed as Eq.8 and Eq.9 with analytically-derived WF matrix forobserved SIs. If we simply apply a linear reconstruction withoutconsideration of the parameter range, negative turbulence strengthoften appears. In order to avoid the situation, Tokovinin et al. (2003)applies the 𝜒 minimization with 𝑦 𝑖 of 𝐽 𝑖 = 𝑦 𝑖 as variable. 𝐽 𝑖 is theturbulence strength of each layer.In this study, the turbulence profile and associated uncertainty areevaluated based on Bayesian inference with Marcov Chain MonteCarlo (MCMC) method. As a prior function of MCMC, we appliedtop-hat filter to limit the parameter space as follows, 𝑃 ( (cid:174) 𝐽 ) = (cid:40) − < log 𝐽 𝑖 [ m / ] < −
11 is satisfied by all 𝐽 𝑖 ,0 otherwise, (17)where (cid:174) 𝐽 is the turbulence profile. The strength range of each tur-bulence layer i.e. − < log 𝐽 𝑖 [ m / ] < −
11 corresponds to2 . × − < 𝑟 [ m ] < . × in the Fried parameter assum-ing measurement wavelength of 500 nm and zenith direction. It isexpected that this parameter range covers the possible turbulencestrength of a single layer. As a likelihood function of MCMC, we usethe probability that the observed SIs are obtained from a Gaussiandistribution with a mean of the expected SIs and a standard deviationof observation errors as follows, 𝐿 ((cid:174) 𝑠 | (cid:174) 𝐽 ) = 𝑀 (cid:214) 𝑚 = (cid:169)(cid:173)(cid:173)(cid:171) √︃ 𝜋𝜎 𝑚 exp (cid:34) − ( 𝑠 𝑚 − ( 𝑊 (cid:174) 𝐽 ) 𝑚 ) 𝜎 𝑚 (cid:35)(cid:170)(cid:174)(cid:174)(cid:172) , (18)where 𝑀 is the number of spatial patterns, (cid:174) 𝑠 and (cid:174) 𝜎 are SIs and theirerrors, respectively, 𝑊 is WF matrix, and (cid:174) 𝐽 is the turbulence profile.By sampling the multi-dimensional space of the (cid:174) 𝐽 effectively usingthe Bayesian inference technique, the strength of turbulence and itserror is reconstructed for each layer. The reconstruction procedurewas conducted utilizing emcee , a MCMC tool for Python . In order to investigate the SH-MASS’s performance of atmosphericturbulence profiling quantitatively, we examine the response of SH-MASS to a single turbulence layer. In this calculation, at first, wecreate a turbulence profile which consists of single turbulence layerat a certain propagation distance. Then, we calculate theoretical SIsby multiplying WF to the turbulence profile assuming the measuredwavelength of 500 nm. The profiling is conducted with a predefinedset of layers. In this study, we select the combinations so that thepropagation distance of the lowest layer should be 0.5 km, propa-gation distance of the highest layer should be 20.0 km and the ratioof the propagation distance of the n-th lowest layer and that of then+1-th lowest layer should be constant. Finally, turbulence profile isreconstructed by using the MCMC method described in Sec.2.3.We repeat this procedure with changing the propagation distance of input single turbulence layer in order to obtain a response func-tion, which is defined as estimated turbulence strength of each re-constructed layer as a function of propagation distance of the inputsingle turbulence layer. Here, we assume that measurement errorsof SIs are 5% of SIs uniformly for all spatial patterns. Accordingto Tokovinin et al. (2003), typical errors of SIs are approximately2% except for 3-7% error for the smallest differential SI. The sizeof the error depends on photon flux. Thus, our assumption of 5%errors would be suitable in order to explore the worst case of typicalperformance of SH-MASS.In Fig.4, the response function of the traditional MASS setup (toppanels) and that of SH-MASS setup (bottom panels) are compared.The assumed spatial pattern for both of the MASS and SH-MASSare same as Fig.2 and Fig.3. Different column represents differentnumber of reconstructed layers. In each panel in Fig.4, x-axis repre-sents the propagation distance of input single turbulence layer whiley-axis represents the ratio of estimated turbulence strength to theinput turbulence strength, hence the sensitivity of MASS to a givenpropagation distance. Each coloured line represents the sensitivityof each reconstructed layer, and errors are defined as the standarddeviation of solutions from 10000 MCMC steps after convergence.Total sensitivity as a sum of the response of all the reconstructionlayers is represented by black dashed line.By comparing the top panels with bottom panels in Fig.4, wecan see that the size of estimation errors are smaller in SH-MASS.In addition, the triangular shape of SH-MASS response function iskept even with the number of reconstructed layers of ten while thetriangular shape is broken and the altitude resolution is poorer in thetraditional MASS case. The black dashed lines are also closer to unityin SH-MASS cases, which means that the estimation of integratedturbulence strength is improved. These results indicate that a largenumber of spatial patterns realized by SH-WFS subaperture geometryis effective in reconstructing turbulence profiles with high altituderesolution and sufficient accuracy.Fig.5 and Fig.6 show how the shape of response function changesif parameters of SH-WFS are changed. In Fig.5, the size of sub-aperture is varied, 1.3 cm, 2.0 cm, and 3.0 cm -diameter cases areshown in top, middle, and bottom panels, respectively. Here, formatof SH-WFS is fixed to 10 ×
10 in all cases. If the size of subaper-ture is increased, fine spatial structure of scintillation is no longerdetectable. Then, scintillation data does not have any informationof turbulence at short propagation distance which is associated withhigh spatial frequency scintillation pattern. For this reason, the prop-agation distance of lowest reconstructed layer is changed so that thepropagation distance ℎ l sec ( 𝑧 ) and the size of subaperture 𝑥 subap satisfy 𝑥 subap ∼ √︁ 𝜆ℎ l sec ( 𝑧 ) , while the propagation distance of thehighest reconstructed layer is fixed to ∼
20 km. These figures in-dicates that the size of subaperture is one of the most importantparameter of SH-MASS, which determines the dynamic range ofthe turbulence profiling. Though SH-MASS with small subaperturesmakes it possible to estimate atmospheric turbulence down to closeto the ground, small subaperture can be suffered from the problemof small number of photons. The diameter of subaperture should bedetermined carefully considering the altitude range of estimation,the magnitude of available star, or the availability of highly sensitivedetector with high pixel read out rate such as an electron multiplyingCCD (EM-CCD).In Fig.6, the format of SH-WFS is varied, 5 ×
5, 10 ×
10, and 15 × MNRAS000
10, and 15 × MNRAS000 , 1–11 (2020)
H.Ogane et al.
Propagation distance from input layer [km]0.00.20.40.60.81.01.2 T r a d i t i o n a l M A SS Propagation distance from input layer [km]
Propagation distance from input layer [km]
10 layers reconstructed Propagation distance from input layer [km]0.00.20.40.60.81.01.2 S h a c n - H a r t m a nn M A SS Propagation distance from input layer [km] 2 Propagation distance from input layer [km]
Figure 4.
Top panels: The response function for the traditional MASS setup in which turbulence strengths of 6, 8, and 10 layers are reconstructed usingscintillation data observed by 4 concentric annular apertures. Bottom panels: The response function for classical SH-MASS setup in which turbulence strengthsof 6, 8, and 10 layers are reconstructed using scintillation data observed by 10x10 SH-WFS whose subaperture diameter corresponds to 1.3 cm on the primarymirror. In each panels, each triangular line represents the response of each reconstruction layer to the input single turbulence layer. In the cases of 6, 8, and 10layers are reconstructed, the propagation distances from the reconstruction layers are [0.5, 1.0, 2.2, 4.6, 9.6, 20.0] km, [0.5, 0.8, 1.4, 2.4, 4.1, 7.0, 11.8, 20.0]km, and [0.5, 0.8, 1.1, 1.7, 2.6, 3.9, 5.8, 8.8, 13.3, 20.0] km, respectively. Black dashed line represents the total sensitivity as a sum of the response of all thereconstruction layers. Gray vertical lines stand for the propagation distance from reconstructed layers, while gray horizontal line stands for the sensitivity whenall input turbulence strength are sensed.
Propagation distance from input layer [km]0.00.20.40.60.81.01.2 . c m s u b a p e r t u r e Propagation distance from input layer [km]
Propagation distance from input layer [km]0.00.20.40.60.81.01.2 . c m s u b a p e r t u r e Propagation distance from input layer [km]2 Propagation distance from input layer [km]0.00.20.40.60.81.01.2 . c m s u b a p e r t u r e Propagation distance from input layer [km]
Figure 5.
Comparison of the SH-MASS response function for the subaperturediameter of 1.3 cm (top), 2.0 cm (middle), and 3.0 cm (bottom) in the caseof 6 (left) and 8 (right) reconstructed layers. Here, the format of SH-WFS is10 ×
10 for all cases. The propagation distance from the lowest reconstructionlayer is changed to 1 km and 2 km for the subaperture size of 2 cm and 3 cm,respectively.
Propagation distance from input layer [km]0.00.20.40.60.81.01.2 × S H - W F S Propagation distance from input layer [km]
Propagation distance from input layer [km]0.00.20.40.60.81.01.2 × S H - W F S Propagation distance from input layer [km]2 Propagation distance from input layer [km]0.00.20.40.60.81.01.2 × S H - W F S Propagation distance from input layer [km]
Figure 6.
Comparison of the SH-MASS response function for the SH-WFSformat of 5 × ×
10 (middle), and 15 ×
15 (bottom) in the caseof 6 (left) and 8 (right) reconstructed layers. Here, the diameter of SH-WFSsubaperture is 1.3cm for all cases.
This is because of large number of constraints realized by larger for-mat. The number of constraints is 43 (15 normals + 28 differentials),285 (51 normals + 234 differentials), and 1740 (106 normals + 1634
MNRAS , 1–11 (2020) tmospheric turbulence profiling with a SH-WFS differentials) for 5 ×
5, 10 ×
10, and 15 ×
15 SH-WFS, respectively.On the other hand, in 15 ×
15 case, MCMC solver does not convergewell in some calculations. This would be because too many con-straints result in complicated posterior probability distribution andcauses poor convergence in Monte Carlo sampling. Considering theresults and computational cost, the format of 10 ×
10 is suitable sizeof SH-WFS.
We investigate the required signal to noise ratio (SNR) in order toconduct the SH-MASS by simulating scintillation observations byMonte Carlo method. In our calculation, we assume one subapertureof SH-WFS with the diameter of 2.0 cm attached to a general tele-scope and wavefront sensor system whose total optical transmissionis assumed to be ∼ 𝑚 V , sky = . [ mag · arcsec − ] , 𝑚 R , sky = . [ mag · arcsec − ] . Here, the number of photons arechanged stochastically so that it should follow Poisson distribution.Additionally, in stellar photon case, the effect from scintillation is alsoconsidered. Assuming SI observed by the subaperture which takesa value from 0.2, 0.5 or 1.0, we randomly vary the number of stel-lar photon following log-normal distribution (Zhu & Kahn 2002).For read-out noise of the detector, we assume stochastic variablewhich follows a Gaussian distribution with the standard deviation of1 . [ ADU / pixel ] . The exposure time is fixed to 2 milliseconds. 𝑛 , thenumber of photon counting, was taken from 3000, 30000, or 300000times, which corresponds to 6 seconds, 1 minute, 10 minutes, re-spectively. After 𝑛 -time photon counting, we calculate observed SIand signal to noise ratio (SNR) using the series. The definitions ofthe SI and SNR in this calculation are as follows,SI = Var [ 𝑛 star ](cid:104) 𝑛 star (cid:105) , (19)SNR = (cid:104) 𝑛 star (cid:105) √︃ (cid:104) 𝑛 star (cid:105) + Var [ 𝑛 sky ] + Var [ 𝑛 ron ] , (20)where (cid:104) 𝑛 star (cid:105) and Var [ 𝑛 star ] are average and variance of count froma simulated star, respectively. Var [ 𝑛 sky ] is variance of count frombackground sky, and Var [ 𝑛 ron ] is variance of count from read-outnoise. We repeat this Monte Carlo simulation of SI and SNR 100times, and evaluated the dispersion of the SIs.Fig.7 shows the ratio of standard deviation of the 100 SIs to averageof the 100 SIs as a function of SNR. In this plot, the number of photoncounting is fixed to 30000. Hence one-minute observation for mea-suring SI is simulated here. Lines with different colours correspondto different input SI, and we can see how the SI measurement errordecreases when the magnitude of observed star increases dependingon the strength of atmospheric turbulence. The curves become flatat SNR > = .
0, the saturated value is lessthan 5%, which means that one-minute observation with SNR > 𝑛 value and fixedinput SI of 0.5. In this plot also, the ratio of the standard deviationto the average of 100 measured SIs become saturated at SNR higherthan ∼
3. The saturated value depends on how many frames are usedto estimate the SI. This is just because more samples make it possibleto more accurately estimate the properties of the SI. Though moreframes make it possible to measure SI with lower measurement error, Signal to noise ratio10 E rr o r o f s c i n t ill a t i o n i n d e x [ % ] SI=0.2, n=3e+4SI=0.5, n=3e+4SI=1.0, n=3e+4
Figure 7.
The SI measurement error (the ratio of the standard deviation to theaverage of 100 measured SIs) as a function of the measured SNR. Lines withdifferent colour means different input SI. As SNR increases, SI measurementerror decreases and saturates to a certain value which correlates with real SIvalue. The saturation occurs at around SNR>3, which would be the requiredSNR value for accurate SI measurement. Signal to noise ratio10 E rr o r o f s c i n t ill a t i o n i n d e x [ % ] SI=0.5, n=3e+3SI=0.5, n=3e+4SI=0.5, n=3e+5
Figure 8.
Same with Fig.7, but lines with different colour means differentnumber of photon counting. As SNR increases, the SI error decreases andsaturates to a certain value which anti-correlates with the number of photoncounting. The saturation occurs at SNR>3, which would be the required SNRvalue for accurate SI measurement. these frames should be obtained within the characteristic timescaleof atmospheric structure evolution. Considering the typical timescaleis ∼
10 minutes and the result of Fig.7, 𝑛 = In order to demonstrate the SH-MASS, we conducted a scintillationobservation using a wavefront sensor system attached to the 50 cmtelescope, IK-51FC in Tohoku University (see Fig.9). Our wavefrontsensor system is consist of a collimator, a Bessel’s R band filter,whose central wavelength is 630 nm, a 150 𝜇 m pitch lenslet array MNRAS000
10 minutes and the result of Fig.7, 𝑛 = In order to demonstrate the SH-MASS, we conducted a scintillationobservation using a wavefront sensor system attached to the 50 cmtelescope, IK-51FC in Tohoku University (see Fig.9). Our wavefrontsensor system is consist of a collimator, a Bessel’s R band filter,whose central wavelength is 630 nm, a 150 𝜇 m pitch lenslet array MNRAS000 , 1–11 (2020)
H.Ogane et al.
Figure 9.
SH-WFS system attached to the 50 cm telescope in Tohoku Uni-versity. A lenslet array is in the aluminium cylindrical tube at the center ofthe figure. The dewar attached to the tube is an EMCCD camera. The camerais cooled down using Peltier device and a liquid cooling system. (Thorlabs, MLA150-5C) which has a chromium mask for blockinglight that reaches outside of the circular aperture of each microlens,relay lenses and an EMCCD camera with E2V CCD60 128 ×
128 24 𝜇 m pixel detector and custom made readout electronics provided byNuvu Cameras. The primary mirror of 50 cm telescope is effectivelydivided into 20 ×
20 by the lenslet array, in other words, the effectivediameter of a subaperture correspond to 2.5 cm on the primarymirror. The field of view of each subaperture is ∼
35 arcsec. Thedetector was cooled down to − ◦ C so that dark current noise can benegligible. One pixel of the detector corresponds to 4.9 arcsec on thesky. Amplification signal of 42.6 V is applied to achieve factor 300multiplication gain of the EMCCD. High speed imaging of 500Hzwas repeated 30000 times targeting Deneb ( 𝑚 𝑅 = . [ mag ] ). Thisprocedure was repeated nine times in one hour on a clear night,October 16th, 2019 in Japan Standard Time. In the one hour, theelevation of the star changed from 46 ° to 34 ° . First step is measurement of count value of each SH-WFS spot ineach frame. The spot reference frame is created by averaging the30000 frames in each dataset, and spot size and locations of spots aremeasured. By fitting each spot of the reference frame with a Gaussianfunction, the diameter of Airy disk is measured to be 4.45 pixels(FWHM is 1.88 pixels). Then, for each frame, we define inside ofcircles with a center of each spot location and a diameter of 5.0 pixelsas spot region and others as background region. Background count isestimated as mean count of the background region. After backgroundsubtraction from all the pixels, each spot count is calculated as thetotal counts of each circle in the spot region. By this procedure,background-subtracted spot counts are calculated for all the spots inthe 30000 frames.Count fluctuation seen in the 30000 frames is ascribed to photonnoise and atmospheric scintillation. Variance of the fluctuation con-tributed from photon noise is the mean photon count, while that con-tributed from scintillation is proportional to the square of the meanphoton count. Therefore, considering the observed photon count of F r a c t i o n [ % ] A = = = Figure 10.
Histogram of count values of one spot of SH-WFS measuredwith 500 Hz in 1 minute. Horizontal axis represents photon counts aftermultiplication by the EM-CCD.The distribution is well fitted by a log-normal distribution, which implies thatthe detected intensity fluctuation is scintillation. ∼ 𝑓 ( 𝑥 ) = 𝐴 √ 𝜋𝜎𝑥 exp (cid:18) ( ln 𝑥 − 𝜇 ) 𝜎 (cid:19) , (21)where 𝜇 and 𝜎 are shape parameters of the distribution and 𝐴 is anormalization parameter. Fig.10 shows the histogram of count valuesof a spot in 1 minute. The histogram is well fitted by a log-normaldistribution function with parameters of 𝐴 = . × , 𝜇 = .
26 and 𝜎 = .
47. All other spots also follow similar shape of the histogram,which supports that the observed intensity fluctuation of SH-WFSspots are caused by atmospheric turbulence. However, it should benoted that the log-normal distribution of the count values is not asufficient condition for concluding that the fluctuation is due to thescintillation.Then, SIs of all the spatial patterns are calculated from the countfluctuations of spots. The mean, variance and covariance of eachspot’s count fluctuation are computed and SIs are calculated follow-ing Eq.15 and Eq.16.Finally, effects from finite exposure time are corrected. In thisobservation, the images are obtained with 2 milliseconds exposure,which means that fluctuation components which have the timescaleof less than 4 milliseconds is smoothed out. Therefore, real SIswhich can ideally observed by 0 millisecond exposure time shouldbe larger than that measured by the above-mentioned procedure. Inthis study, we follow the method described in Tokovinin et al. (2003)in which 0 millisecond SIs are estimated from linear extrapolation ofSIs measured by 𝜏 milliseconds exposure and that measured by 2 𝜏 milliseconds exposure, 𝑠 = 𝑠 𝜏 − 𝑠 𝜏 . (22)Here, datasets with 2 𝜏 milliseconds exposure is effectively obtainedby averaging two adjacent images in the data of 𝜏 milliseconds ex-posure.Top panel of Fig.11 shows the observed normal SIs as a functionof separation of two subapertures which constitutes a normal spatial MNRAS , 1–11 (2020) tmospheric turbulence profiling with a SH-WFS N o r m a l s c i n t ill a t i o n i n d e x N u m b e r o f p a tt e r n s Figure 11.
Top: The observed normal SIs are plotted as a function of a separa-tion distance of two subapertures which constitutes a spatial pattern. Differentcolour represents different observation time. The value of SIs decreases andbecome flattened as the separation increases, which indicates that the typicalcorrelation length of scintillation is shorter than 15 cm. And the trend thatSIs become larger as time goes by can be explained by variation of the ele-vation angle of the star Bottom: The number of subaperture pairs which havecommon separation distance. The normal SIs and their errors in top figurehave been calculated as the mean and standard deviation of these number ofstatistics. pattern. Each colour of plots represents the difference of observedtime. The error on each normal SI is small enough to discern thedifference between each scintillation state of each observation time.Bottom panel of Fig.11 shows the number of subaperture pairs whichhave common spatial pattern. At all observation time, the normal SIdecreases as a function of subaperture separation and get flattened atspatial length of 10-15 cm and longer. This feature reflects that thereare little atmospheric turbulence at higher than ∼
20 km. In fact, thespatial scale of 10-15 cm is consistent with the typical spatial scaleof scintillation created by turbulence layer i.e. √︁ 𝜆ℎ sec ( 𝑧 ) ∼ . 𝜆 ∼
600 nm, altitude ℎ ∼
20 km,and airmass sec ( 𝑧 ) ∼ .
5. In addition, this feature can be understoodusing Eq.15. For null separation case, Cov [ 𝐼 𝑖 , 𝐼 𝑗 ] becomes Var [ 𝐼 𝑖 ] ,and SI becomes Var [ 𝐼 𝑖 ]/(cid:104) 𝐼 𝑖 (cid:105) . Whereas for very long separationcase, Cov [ 𝐼 𝑖 , 𝐼 𝑗 ] becomes 0, and SI becomes Var [ 𝐼 𝑖 ]/ (cid:104) 𝐼 𝑖 (cid:105) , halfof the SI for the null separation case. In Fig.11 actually, normalSIs for longer separation than 15 cm is almost half of normal SIfor 0 cm separation. Besides, there is a trend that the value of SIbecome larger as time goes by. This can be explained by the changeof the elevation of the star. As the elevation becomes lower, apparentaltitude of turbulence layer becomes higher. Both effects accountfor the increase of SIs. These properties also support that detectedfluctuation of stellar intensity is due to the atmospheric turbulence. Fig.12 shows the atmospheric turbulence profile reconstructed fromthe observed SIs by the MCMC estimation method described in A l t i t u d e [ k m ] Measurement time: 22:39 = 3.098-layer fit (mcmc), = 3.01 Measurement time: 22:46 = 6.218-layer fit (mcmc), = 6.14 Measurement time: 22:51 = 3.388-layer fit (mcmc), = 3.35 A l t i t u d e [ k m ] Measurement time: 22:56 = 6.858-layer fit (mcmc), = 6.76 Measurement time: 23:11 = 4.778-layer fit (mcmc), = 4.75 Measurement time: 23:27 = 6.868-layer fit (mcmc), = 6.72 Cn dh[10 m ] A l t i t u d e [ k m ] Measurement time: 23:33 = 2.078-layer fit (mcmc), = 1.98 Cn dh[10 m ] Measurement time: 23:39 = 3.048-layer fit (mcmc), = 2.97 Cn dh[10 m ] Measurement time: 23:50 = 2.278-layer fit (mcmc), = 2.19 Figure 12.
The atmospheric turbulence profile reonstructed by SIs measuredat Tohoku university. These profiles are reconstructed by the MCMC estima-tion method mentioned in Sec.2.3. The effect of elevation angle of the staris corrected. Different panel corresponds to different observation time whichis described in each panel’s title. Blue lines are a profiles which are recon-structed assuming six layers while orange lines assume eight layers. Reduced 𝜒 values of each profile estimation are shown in the legend. Sec.2.3. Different panel corresponds to different observation timewhich is described in each panel’s title. The propagation distances ofreconstructed layers are varied so that the reconstruction altitudes ofthe layers should be same among observation times. Blue lines areprofiles which are estimated assuming six layers (altitude of 1.0, 1.8,3.3, 6.0, 11.0, 20.0 km above the telescope aperture) while orangelines assume eight layers (altitude of 1.0, 1.5, 2.4, 3.6, 5.5, 8.5, 13.0,20.0 km above the telescope aperture). The reduced 𝜒 values whichare less than 10 in all datasets and the small uncertainties whichrepresent one sigma values of turbulence strengths after MCMCconvergence implies that the observed normal and differential SIs aredescribed well with the scintillation model. In addition, the overallshape of the profiles shows that the strongest turbulence exists atthe lowest layers and second strongest peak distributes at roughly 10km. The profiles are consistent with that expected from the typicalcharacteristics of the Earth’s atmosphere such as the ground turbulentlayer and the tropopause, respectively. MNRAS000
The atmospheric turbulence profile reonstructed by SIs measuredat Tohoku university. These profiles are reconstructed by the MCMC estima-tion method mentioned in Sec.2.3. The effect of elevation angle of the staris corrected. Different panel corresponds to different observation time whichis described in each panel’s title. Blue lines are a profiles which are recon-structed assuming six layers while orange lines assume eight layers. Reduced 𝜒 values of each profile estimation are shown in the legend. Sec.2.3. Different panel corresponds to different observation timewhich is described in each panel’s title. The propagation distances ofreconstructed layers are varied so that the reconstruction altitudes ofthe layers should be same among observation times. Blue lines areprofiles which are estimated assuming six layers (altitude of 1.0, 1.8,3.3, 6.0, 11.0, 20.0 km above the telescope aperture) while orangelines assume eight layers (altitude of 1.0, 1.5, 2.4, 3.6, 5.5, 8.5, 13.0,20.0 km above the telescope aperture). The reduced 𝜒 values whichare less than 10 in all datasets and the small uncertainties whichrepresent one sigma values of turbulence strengths after MCMCconvergence implies that the observed normal and differential SIs aredescribed well with the scintillation model. In addition, the overallshape of the profiles shows that the strongest turbulence exists atthe lowest layers and second strongest peak distributes at roughly 10km. The profiles are consistent with that expected from the typicalcharacteristics of the Earth’s atmosphere such as the ground turbulentlayer and the tropopause, respectively. MNRAS000 , 1–11 (2020) H.Ogane et al.
Although MCMC-based profile estimation method enables us to eval-uate the estimation error, it requires large calculation cost, whichtypically takes a few tens of minutes for six-layer reconstruction witheight-core parallel processing using
Intel ® Core ™ i7-4790K CPUand highly depends on the number of reconstructed layers. On theother hand, atmospheric turbulence profile as a prior information fortomographic reconstruction matrix has to be updated in a timescale oftens of minutes. Then, we try faster profile calculation based on Broy-den–Fletcher–Goldfarb–Shanno (BFGS) algorithm, which is one ofiterative solvers for non-linear optimization problem. This algorithmcan be utilized with scipy.optimize.minimize module for
Python . Weimpose the same condition of − < log 𝐽 𝑖 [ m / ] < −
11 for allthe component of (cid:174) 𝐽 as MCMC-based method, and minimize the 𝜒 function directly.However, in iterative calculation method, the solution is not nec-essarily the global minimum. Actually, in the current case, one-timeiterative calculation does not give an identical solution. Hence, weconduct the BFGS algorithm 1000 times from 1000 different randominitial turbulence profile. Then, we pick out 100 final turbulence pro-files whose 𝜒 values are the smallest and calculated the mean andstandard deviation of the 100 profiles.In Fig.13, we compare the turbulence profile obtained by BFGSiterative method with that obtained by the MCMC method. In allobservation time, both estimation methods reproduce the same tur-bulence profile. The consistency suggests that 1000-time iterativeminimization from the random initial profiles is sufficient to findout the global minimum. Because the calculation time for the 1000-time iterative minimization is typically a few minutes, the iterativeBFGS method can be used for a faster profile reconstruction. Then,we conduct the profile estimation with 10 layers (1.0, 1.4, 1.9, 2.7,3.8, 5.3, 7.4, 10.3, 14.3, 20.0 km), which takes the timescale of dayswhen MCMC-based method is used. The result is shown in Fig.14.Here, higher altitude resolution with 𝑑ℎ / ℎ = . 𝑑ℎ / ℎ = . As mentioned in previous works (e.g. Avila et al. 1997), turbulence atless than several hundred meter high is undetectable by scintillation-based profiling methods. This is because the variance of observedintensity, or scintillation is proportional to the propagation distancewith the power of 5/6th. This characteristic can be seen in our resultof response function (Fig.4, Fig.5 and Fig.6) in which the estimationerror of turbulence strength becomes larger for shorter propagationdistance of input layer. In order to overcome this limitation, threesolutions can be considered as future improvement plans. The first isgeneralized mode. Likewise G-SCIDAR, turbulence at low altitudecan be measured by placing SH-WFS at some distance away fromthe pupil plane. The second is DIMM mode. In DIMM, two imagemotions of single bright star is measured by two apertures whosecenter separation is typically 20-30 cm. By using SH-WFS spots, starimage motion is observed by various separation length. Therefore, A l t i t u d e [ k m ] Measurement time: 22:39 = 3.098-layer fit (mcmc), = 3.016-layer fit (bfgs), = 3.088-layer fit (bfgs), = 3.00 Measurement time: 22:46 = 6.218-layer fit (mcmc), = 6.146-layer fit (bfgs), = 6.208-layer fit (bfgs), = 6.13 Measurement time: 22:51 = 3.388-layer fit (mcmc), = 3.356-layer fit (bfgs), = 3.388-layer fit (bfgs), = 3.35 A l t i t u d e [ k m ] Measurement time: 22:56 = 6.858-layer fit (mcmc), = 6.766-layer fit (bfgs), = 6.848-layer fit (bfgs), = 6.75 Measurement time: 23:11 = 4.778-layer fit (mcmc), = 4.756-layer fit (bfgs), = 4.778-layer fit (bfgs), = 4.75 Measurement time: 23:27 = 6.868-layer fit (mcmc), = 6.726-layer fit (bfgs), = 6.858-layer fit (bfgs), = 6.70 Cn dh[10 m ] A l t i t u d e [ k m ] Measurement time: 23:33 = 2.078-layer fit (mcmc), = 1.986-layer fit (bfgs), = 2.068-layer fit (bfgs), = 1.97 Cn dh[10 m ] Measurement time: 23:39 = 3.048-layer fit (mcmc), = 2.976-layer fit (bfgs), = 3.048-layer fit (bfgs), = 2.96 Cn dh[10 m ] Measurement time: 23:50 = 2.278-layer fit (mcmc), = 2.196-layer fit (bfgs), = 2.268-layer fit (bfgs), = 2.19 Figure 13.
Comparison of atmospheric turbulence profile reconstructed bythe MCMC based method (blue and orange lines) and iterative method (greenand red lines). The effect of elevation angle of the star is corrected. It canbe understood that iterative method reproduces same profile as the MCMCmethod. Different panel corresponds to different observation time which isdescribed in each panel’s title. Reduced 𝜒 values of solutions from eachmethod are shown in the legend. DIMM can be conducted by using SH-WFS system and turbulencestrength at ground layer can be estimated by the difference betweenDIMM measurement and MASS measurement. The third solutionis a combination with SLODAR. SLODAR uses two SH-WFSs tomeasure the correlation of wavefront distortion from a double star.Because triangulation-based profiler like SLODAR does not haveany sensitivity to turbulence at high altitude, some SLODAR systemis optimized for profiling of the ground layer (Butterley et al. 2020).Combining SH-MASS with SLODAR enables us to profile wholeatmospheric turbulence using a single optical system.
In this study, we investigate a new MASS-based atmospheric tur-bulence profiling method called SH-MASS, which reproduces theprofile from scintillation observed by a Shack-Hartmann wavefrontsensor. By evaluating the response function of the SH-MASS incomparison with those of the traditional MASS, it is shown that SH-MASS theoretically has higher altitude resolution than traditionalMASS under the assumption that the scintillation measurementshave 5% error. This high altitude resolution is enabled by a largenumber of spatial patterns realized by the grid pattern of SH-WFS.
MNRAS , 1–11 (2020) tmospheric turbulence profiling with a SH-WFS A l t i t u d e [ k m ] Measurement time: 22:39 = 3.088-layer fit (bfgs), = 3.0010-layer fit (bfgs), = 3.00 Measurement time: 22:46 = 6.208-layer fit (bfgs), = 6.1310-layer fit (bfgs), = 6.15 Measurement time: 22:51 = 3.388-layer fit (bfgs), = 3.3510-layer fit (bfgs), = 3.36 A l t i t u d e [ k m ] Measurement time: 22:56 = 6.848-layer fit (bfgs), = 6.7510-layer fit (bfgs), = 6.78 Measurement time: 23:11 = 4.778-layer fit (bfgs), = 4.7510-layer fit (bfgs), = 4.75 Measurement time: 23:27 = 6.858-layer fit (bfgs), = 6.7010-layer fit (bfgs), = 6.75 Cn dh[10 m ] A l t i t u d e [ k m ] Measurement time: 23:33 = 2.068-layer fit (bfgs), = 1.9710-layer fit (bfgs), = 1.98 Cn dh[10 m ] Measurement time: 23:39 = 3.048-layer fit (bfgs), = 2.9610-layer fit (bfgs), = 2.97 Cn dh[10 m ] Measurement time: 23:50 = 2.268-layer fit (bfgs), = 2.1910-layer fit (bfgs), = 2.21 Figure 14.
Atmospheric turbulence profile at Tohoku university. These pro-files are reconstructed from the observed SIs by the iterative estimationmethod mentioned in Sec.5.1. The effect of elevation angle of the star iscorrected. Different panel corresponds to different observation time whichis described in each panel’s title. Green, red and purple lines are a profileswhich are reconstructed assuming six, eight and ten reconstructed layers,respectively. Reduced chi squared values of each profile estimation are de-scribed in the legend.
By investigating the behaviour of response functions with changingthe parameters of SH-MASS, the larger size of subaperture makeslower sensitivity to the low altitude turbulence, therefore smaller sizeof subaperture is better as long as the signal to noise ratio of eachspot in a SH-WFS image is larger than ∼ ∼
10 minutes, we confirmthat faster iterative method can also reproduce the same profile as theMCMC-based method.
ACKNOWLEDGEMENTS
The authors thank Drs. Yosuke Minowa, Kazuma Mitsuda, and KokiTerao for helpful discussions. HO is supported by Graduate Programon Physics for the Universe (GP-PU), Tohoku University. MA issupported by JSPS KAKENHI (17H06129). Part of this work wasachieved using the grant of Joint Development Research supportedby the Research Coordination Committee, National Astronomical Observatory of Japan (NAOJ), National Institutes of Natural Sciences(NINS).
DATA AVAILABILITY
The data underlying this article will be shared on reasonable requestto the corresponding author.
REFERENCES
Arsenault R., et al., 2012, in Adaptive Optics Systems III. p. 84470JAvila R., Vernin J., Masciadri E., 1997, Applied Optics, 36, 7898Beckers J. M., 1988, in European Southern Observatory Conference andWorkshop Proceedings. p. 693Butterley T., Wilson R., Sarazin M., Dubbeldam C., Osborn J., Clark P., 2020,Monthly Notices of the Royal Astronomical Society, 492, 934Costille A., Fusco T., 2012, in Adaptive Optics Systems III. p. 844757Farley O., Osborn J., Morris T., Fusco T., Neichel B., Correia C., Wilson R.,2020, Monthly Notices of the Royal Astronomical Society, 494, 2773Fusco T., Costille A., 2010, in Adaptive Optics Systems II. p. 77360JGendron E., Morel C., Osborn J., Martin O., Gratadour D., Vidal F., Le LouarnM., Rousset G., 2014, in Adaptive Optics Systems IV. p. 91484NGilles L., Wang L., Ellerbroek B., 2008, in Adaptive Optics Systems. p.701520Guesalaga A., Perera S., Osborn J., Sarazin M., Neichel B., Wilson R., 2016,in Adaptive Optics Systems V. p. 99090HGuesalaga A., Ayanc ¥ ’an B., Sarazin M., Wilson R., Perera S., Le LouarnM., 2020, Monthly Notices of the Royal Astronomical SocietyHammer F., et al., 2004, in Second backaskog workshop on extremely largetelescopes. pp 727–736Kornilov V., Tokovinin A. A., Vozyakova O., Zaitsev A., Shatsky N., PotaninS. F., Sarazin M. S., 2003, in Adaptive Optical System Technologies II.pp 837–845Kornilov V., Tokovinin A., Shatsky N., Voziakova O., Potanin S., Safonov B.,2007, Monthly Notices of the Royal Astronomical Society, 382, 1268Lardi ¥ ‘ere O., et al., 2014, in Adaptive Optics Systems IV. p. 91481GMarchetti E., et al., 2007, The Messenger, 129Minowa Y., et al., 2017, Adaptive Optics for Extremely Large Telescopes V(AO4ELT5)],(Jun 2017)Rigaut F., 2002, in European Southern Observatory Conference and Work-shop Proceedings. p. 11Rigaut F., Neichel B., 2018, Annual Review of Astronomy and Astrophysics,56, 277Rigaut F., et al., 2014, Monthly Notices of the Royal Astronomical Society,437, 2361Rocca A., Roddier F., Vernin J., 1974, JOSA, 64, 1000Sarazin M., Roddier F., 1990, Astronomy and Astrophysics, 227, 294Saxenhuber D., Auzinger G., Le Louarn M., Helin T., 2017, Applied Optics,56, 2621Stone J., Hu P., Mills S., Ma S., 1994, JOSA A, 11, 347Tallon M., Foy R., 1990, Astronomy and Astrophysics, 235, 549Tokovinin A. A., 2003, JOSA A, 20, 686Tokovinin A., 2004, Publications of the Astronomical Society of the Pacific,116, 941Tokovinin A., Kornilov V., Shatsky N., Voziakova O., 2003, Monthly Noticesof the Royal Astronomical Society, 343, 891Vidal F., Gendron E., Rousset G., 2010, JOSA A, 27, A253Wilson R. W., 2002, Monthly Notices of the Royal Astronomical Society,337, 103Zhu X., Kahn J. M., 2002, IEEE Transactions on communications, 50, 1293This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000