1. Introduction

The turbulence in the Earth’s atmosphere affects the propagation of light by causing images of point sources to appear to rapidly change in brightness (scintillation) and break up into speckles. This optical turbulence is caused by the mixing of air of different temperatures, and hence density, developing spatial and temporal variations in refractive index. Atmospheric optical turbulence is a major limitation to several mature and emerging applications, such as; Free Space Optical Communications (FSOC) [1,2], either horizontally or to and from space and including long-range Quantum Key Distribution (QKD) [3,4]; imaging of space objects for Space Surveillance and Tracking (SST) [5]; space debris deorbiting via acceleration due to laser ablation or photon pressure [6]; and both daytime and night-time ground-based astronomy (for example, [7]).

For each of these applications knowing the strength and variability of the optical turbulence enables us to select the best possible sites, model how each system can perform, optimise the optical design, and develop novel turbulence mitigation schemes, for example using Adaptive Optics (AO) technology.

Continuous 24-hour turbulence monitoring is required to support FSOC and SST, where facilities will need to operate continuously, night and day. FSOC will also need to operate in urban and suburban locations, environments that are likely to have stronger turbulence than the pristine mountain top sites of most astronomical observatories where turbulence monitoring instrumentation tends to be located.

Even for astronomy the ability to measure the turbulence continuously enables data assimilation into turbulence forecasting tools. This continuous data assimilation has been shown to significantly improve the performance of turbulence forecasting tools enabling a much more efficient smart scheduling of observations [8], ensuring that the most sensitive observations are carried out at the optimum time, maximising the probability of success.

Here, we demonstrate for the first time that a system that exploits a Shack-Hartmann wavefront sensor (WFS) operating in the short-wave infrared (SWIR) and utilises techniques for wavefront-sensing in high noise can provide continuous 24-hour vertical turbulence monitoring even in strong turbulence conditions. The system is built around a commercial off-the-shelf 11" telescope and the low weight and size make it ideal for portable campaigns anywhere in the world. The instrument, which we call the 24-hour Shack-Hartmann Image Motion Monitor (24hSHIMM), monitors a single bright star and can extract the atmospheric turbulence coherence length, $r_0$, coherence angle, $\theta _0$, coherence time, $\tau _0$, the Rytov variance, $\sigma _R^2$, and a four-layer vertical profile. This system is a development of the SHIMM [9]. However, until now, the SHIMM has been limited to night-time operations only. The work presented here describes a significant upgrade to existing systems in both the hardware and analysis and represents a necessary advance in the field of turbulence monitoring. Current techniques are limited to either day or night-time, leading to differential instrument bias between the two windows and a significant gap in measurement during twilight hours. They are also generally restricted to weak turbulence, limiting their use for monitoring conditions in urban environments, for example near data and population centres, a requirement for free-space optical communications networks.

2. 24hSHIMM

In selecting equipment for the 24hSHIMM, an effort was made to maintain portability while maximising the contrast of daytime WFS images with the vast majority of equipment obtained off-the-shelf. Figure 1 shows the assembled instrument. For the telescope a Celestron C11 279.4mm SCT was chosen to maximise mirror collecting area while retaining portability and ease of setup. This was paired with a low-cost but programmable Celestron CGX mount. The daytime brightness of the sky background decreases steeply with increasing wavelength as a result of Rayleigh scattering, therefore optical components were optimised for SWIR observations. For autoguiding, a SkyWatcher Evoguide 50ED guidescope was paired with a ZWO ASI462MC CMOS camera (60% QE at 850 nm) and an 850 nm longpass filter to block visible wavelengths. The WFS detector was a First Light Imaging CRED-3 Indium Gallium Arsenide (InGaAs) camera [10] with a bandpass of 900-1700 nm. The optical design of the 24hSHIMM is shown in Fig. 2. We used a NIR-II coated achromatic doublet for the collimating lens, 500 µm-pitch square lenslet array and a 0.8 mm minimum diameter (corresponding to a minimum field-of-view of 59^″) iris to function as a field stop minimising background light. The collimating lens and microlens array focal lengths were chosen to produce six 4.7 cm-wide sub-apertures across the pupil. This microlens array sub-aperture size coupled with our choice of camera ensured that WFS focal spots were at least Nyquist-sampled and motions due to wavefront distortion could be measured accurately. This was verified in simulation as detailed in section 4.1.

 

Fig. 1. Left, a photo of the 24hSHIMM operating during daytime on La Palma, Spain. Right, a typical WFS frame from the 24hSHIMM taken at approximately 2pm local time in La Palma. The red squares indicate the sub-aperture focal spots used for data analysis.

Download Full Size | PDF

 

Fig. 2. Simple 24hSHIMM optical design. $f_\mathrm {tel},f_\mathrm {col},f_\mathrm {mla}$ are the focal lengths of the telescope, collimating lens and microlens array respectively and the $d_\mathrm {col},d_\mathrm {mla}$ and $D$ are the diameters of the collimated beam, the microlens array lenslets and the telescope pupil respectively.

Download Full Size | PDF

The instrument was aligned with the aid of a 3d Autodesk Inventor model of the optomechanics integrating a model of the optics produced in Zemax. After achieving focus, alignment errors were accounted for by performing an image scale calibration on a double star with a known angular separation. This gave an average image scale of $\rho = {2.27\pm 0.01}^{\prime \prime }$ pixel$^{-1}$ compared to the predicted ${2.17}^{\prime \prime }$ pixel$^{-1}$. The effective wavelength of the system, a key variable for the calculation of atmospheric parameters, was found to be 1280 nm. The dependence of the effective wavelength on the stellar spectrum was ignored. The value of the effective wavelength was estimated by constructing a theoretical transmission curve based on manufacturer data for the optics and camera and a study into the Celestron corrector plate coating [11] which we use as an estimate of the transmission of the three surfaces in the telescope. The last assumption is likely to result in a systematic error in the effective wavelength calculation, however this was the best possible estimate in the absence of telescope transmission data at SWIR wavelengths.

The experimental setup consisted of a Linux workstation running the cameras, mount and data acquisition software and could be remotely configured and controlled. The process by which data were collected was as follows. Initially stars with a zenith angle greater than 40° were filtered out and the remaining stars were ordered by J-band magnitude. The best target was then judged primarily on brightness, but also on separation from the sun and total observation time. The 24hSHIMM was then made to track the best available star and the WFS camera set to run at a short exposure time of 2ms to ensure that wavefront measurements were not smeared over the course of the exposure [12]. For a typical wind speed of 20 m s⁻¹ this will produce a blur of 4cm, which is comparable to the sub-aperture size, and bias caused by the finite exposure time is ignored. 3600 images were recorded at approximately 120 Hz constituting a single data point or measurement. After each measurement, if necessary, the star could be guided back into the centre of the field of view using the autoguiding equipment. This process was repeated until either the target exceeded the zenith angle threshold or the mount would be forced perform a meridian flip and a new target had to be found. With a robust algorithm for determining optimal targets this process may be fully automated to gather large quantities of continuous data.

To investigate target availability for the 24hSHIMM, the minimum required signal-to-noise ratio was estimated by running simulations for a variety of shot signals using detector noise parameters quoted in section 4.1. This investigation returned a shot signal of approximately 5000$e^-$ corresponding to a maximum J band magnitude of 1.5 for the target star. Using data from a star catalogue [13] and this magnitude limit, the number of valid target stars at each hour was investigated for the date and location of the 36-hour experiment in La Palma. It was found that at any given time there would be a minimum of 20 targets available. However, due to limitations in the mount pointing accuracy, the guidescope, which was significantly more limited by sky background noise, was relied upon to acquire targets. The magnitude limit used in the experiment was zero during the day - still providing at least one usable target at all times.

3. Methods

3.1 Shack-Hartmann daytime image analysis

The centroids of pixel values of the focal spots in the Shack-Hartmann lenslet array are calculated via a version of the moving window centre-of-gravity algorithm [14] using a top-hat window in conjunction with a brightest pixel centre-of-gravity calculation [15] to find the centroid position within the moving window. The centroids not only give the wavefront slope across each Shack-Hartmann sub-aperture, but also enable accurate measurement of the intensity in the spots, which is here achieved by summing the pixel values in a small circle surrounding the centroid.

Without a careful background subtraction, the centroid and intensity of the Shack-Hartmann spot cannot be calculated accurately and the rapidly changing sky background at sunset and sunrise presents a significant challenge as any on-sky background images can quickly become obsolete. Background subtraction for the 24hSHIMM is therefore implemented in the real-time data analysis stage. The mean background level is measured in a thin annulus around each sub-aperture spot, which is then subtracted from each image. When calculating the spot centroids, a 2x2 median filter is applied to the images to reduce centroid bias from the residual noise and and hot pixels present. This filter is not applied to images when measuring the spot intensities.

Although the use of an InGaAs camera minimises sky background noise during the daytime, these sensors have high levels of fixed pattern noise and dark noise, especially when uncooled. Indeed during experiments with the 24hSHIMM it was observed that the combined quadrature sum of detector readout noise variance and dark noise variance was $66.5 e^- \mathrm {pix}^{-1}$. This was greater than the daytime sky background noise. The CRED-3 camera was chosen due to its high linearity and suitability for wavefront sensing applications [10]. The camera is uncooled - reducing weight and expense, however by running at the very short exposure times required for wavefront sensing, the increased contribution of dark current to the detector noise as a result of higher sensor temperature is minimised. Hot pixels are calibrated through a factory-set mask provided with the camera. Furthermore to address the high levels of fixed pattern noise and dark noise, master-dark images are taken at regular intervals throughout experiments as the sensor temperature changes. Additionally the 24hSHIMM utilises a small 224$\times$224 pixel region of the sensor minimising the small effects of sensor non-linearity on centroid measurements.

3.2 Four-layer vertical profile

The vertical profile of the optical turbulence strength, $C^{2}_{n}(h)\textrm{d}h$, is estimated from the auto-covariance of the intensity in the Shack-Hartmann sub-apertures. This is similar to the SCO-SLIDAR (Single Coupled SLODAR SCIDAR) method [16,17] but only utilising the intensity covariance information, removing issues associated with telescope vibrations in the phase covariance. This is also the approach taken by the Shack-Hartmann MASS instrument [18]. Adopting here similar notation to [19], the vector of measured intensity covariances is denoted $\mathbf {C}$, the unbiased covariances $\mathbf {C'}$, the matrix of weighting functions $\mathbf {M}$ and the unknown $C^{2}_{n}(h)\textrm{d}h$ vector $\mathbf {S}$. The weighting functions $\mathbf {M}$ are computed using standard weak-scintillation theory for monochromatic light at the effective wavelength of the instrument, not accounting for the effects of finite spectral bandwidth and partial saturation. Additionally a vector, $\mathbf {E}$, containing the bias to the intensity covariances as a result of photometric noise is introduced and the inversion problem is given by

It is assumed that the photometric noise is independent in each sub-aperture, and therefore that intensity covariances between different sub-apertures are unbiased and the corresponding components of $\mathbf {E}$ are equal to zero. For intensity covariances between a sub-aperture and itself, denoted components $k$ in $\mathbf {C}$, the bias as a result of photometric noise, and therefore the corresponding components $k$ in $\mathbf {E}$, are

For a single exposure, $S$ is the sub-aperture shot signal from the target star, $n_{\mathrm {pix}}$ is the number of pixels used to measure the spot intensity and its typical value in this work is 29, $B$ and $D$ are the sky background and dark current counts per pixel, $\sigma _{\mathrm {Rd}}$ is the RMS readout noise per pixel and the angular brackets denote an average over all frames in the measurement sequence. Equation (1) is then solved for $\mathbf {S}$ using a non-negative least squares algorithm. The uncertainties in the values of the reconstructed vertical turbulence profile, $C^{2}_{n}(h)\textrm{d}h$, are estimated through the bootstrap method [20]. This involves randomly selecting (with replacement) frames from a measurement, performing the profile reconstruction from this sample and repeating this a large number of times. The standard deviation of the resulting sampled profiles is therefore the bootstrap estimation of the standard error in the reconstructed profile.

However, the key limitation of this approach is that the intensity is insensitive to the ground layer turbulence. Therefore, the ground layer strength is found by subtracting the sum of this profile from the integrated turbulence strength obtained from WFS slopes. This is generally a similar approach to [21] and the Multi-Aperture Scintillation Sensor - Differential Image Motion Monitor [22], except in this analysis, as explained in section 3.3, the integrated turbulence strength is found by comparing the auto-covariance of WFS angle of arrival slope measurements with a theoretical response. The ground layer is placed at 0 m and the vertical heights of the layers reconstructed by solving Eq. (1) are chosen to be 4, 12, and 20 km.

3.3 Atmospheric turbulence parameters

The 24hSHIMM takes advantage of the SLOpe Detection And Ranging (SLODAR) angle of arrival auto-covariance analysis [23] to measure the turbulence coherence length (or Fried parameter), $r_0$, which is related to the integrated turbulence strength. This approach involves calculating the theoretical centroid auto-covariance response of the WFS for imaging through Kolmogorov turbulence of a given strength and fitting it to measured auto-covariances from the 24hSHIMM [24,25]. The estimate of the coherence length is then further refined by scintillation correction using the four-layer vertical profile as described in section 3.4.

The remaining atmospheric turbulence parameters are calculated from the four-layer vertical $C^{2}_{n}(h)\textrm{d}h$ profile. Equations used to calculate the coherence angle, Rytov variance and coherence time are found in [7]. The coherence angle, $\theta _0$, represents the largest angle on-sky over which AO corrections may be considered valid. The Rytov variance, $\sigma _R^2$, defines the magnitude of the received intensity variance after passing through the atmosphere. This is equivalent to the scintillation index for a point in weak turbulence ($\sigma _R^2<0.3$). These two parameters may be calculated directly from the measured $C^{2}_{n}(h)\textrm{d}h$ profile. To estimate the optical turbulence coherence time, the vertical profile of the wind speed is also required as faster turbulence will lead to shorter coherence time. The vertical profile of the wind velocity, $V(h)$ is extracted from the reanalysis of the ERA5 meteorological forecast [26] with the ground layer being replaced by a measurement from a local anemometer. It has been shown that meteorological forecast wind velocity is consistent with optical turbulence velocity [27]. The ERA5 reanalysis has a spatial resolution of 0.25° in longitude and latitude and a time resolution of one hour. The error in the wind speed profile can be obtained from the ensemble standard deviation of the forecasts and it was found that the median fractional error for the 36-hour data set was 13%. To the best of our knowledge, this is the first time such a hybrid approach has been taken for turbulence monitoring.

3.4 Scintillation correction

During the process of computing the theoretical centroid auto-covariance response of the 24hSHIMM, following the analysis in [25], the theoretical centroid covariances between pairs of sub-apertures on the WFS are computed numerically. This involves evaluating the spatial covariance of the phase aberration as a result of the turbulence,

4. Results

4.1 Validation in simulation

Figure 3 shows the turbulence parameters estimated in an end-to-end Monte-Carlo simulation of the methods presented above carried out using the AOtools simulation package [29]. The input turbulence profiles are taken from the Stereo-SCIDAR database at Paranal, Chile [30]. The simulation is a full physical optics simulation capable of modelling strong turbulence propagation. The simulation assumes monochromatic light at the effective wavelength of the 24hSHIMM, 1280 nm, as this was found to be a very good approximation to running the simulation with multiple wavelengths weighted by the instrument transmission spectrum. The effects of finite exposure time have also been neglected. To include the stronger low-altitude turbulence expected during the daytime as a result of solar heating, an additional high-strength layer of turbulence has been placed at an altitude of 0 m. The simulations have been carried out for a mean target signal level per spot of 35000$\mathrm {e}^{-}$ and a mean sky background level of 2850 $\mathrm {e}^- \mathrm {pix}^{-1}$. Also included are the shot noise of the source flux and the sky background light based on the mean signal and background flux levels in addition to an RMS detector read noise of 66.5 $\mathrm {e}^- \mathrm {pix}^{-1}$. These values were obtained from WFS images taken during the 36-hour experiment in La Palma and give a signal-to-noise ratio of 70.6. The shot signal is equivalent to that of the faintest target observed during the day and the simulation sky background and the read noise are equivalent to the median values observed during the day. The gain of the camera is $G = 2.01 \mathrm {e}^- \mathrm {ADU}^{-1}$. It should be noted that the RMS read noise here is the quadrature sum of the detector dark noise variance and readout noise variance because a separate measurement of the dark count was impossible as only the dark-subtracted images were saved. For the analysis of simulations, a 2x2 median filter was applied when calculating the centroids to reduce the bias in measurements due to strong noise and the intensities were measured from unfiltered images.

 

Fig. 3. Parameter measurement methods are tested via end-to-end Monte Carlo simulation of the instrument using measured vertical turbulence profiles from the Stereo-SCIDAR instrument as input. The panels show the coherence length, coherence angle on the top row and the coherence time and Rytov variance on the bottom row. The red dashed line shows the linear best-fit for the data calculated through linear regression, whereas the black solid line indicates the perfect instrument response. r indicates the value of the Pearson correlation coefficient for the data. Simulations were carried out for monochromatic light with a wavelength of 1280 nm, however parameters reported in this figure have been corrected to their values at 500 nm.

Download Full Size | PDF

In Fig. 3, the "input profile" atmospheric coherence time is calculated using interpolated wind speed profiles from the stereo-SCIDAR data and a perfect knowledge of this $V(h)$ profile by the 24hSHIMM is assumed. For the reconstructed profile, the wind speed of the 0 m layer is set as $V(0)$ and for the three atmospheric layers the mean value of $V(h)^{5/3}$ in an 8 km bin around each layer is used. This analysis relies on the assumption that these bins are large enough to contain a uniform distribution of turbulence strength. However the uncertainty in the value of the $V(h)$ depends on the true distribution of $C_n^2(h)$ and $V(h)$ within each bin and so is not included as the 24hSHIMM cannot quantify this, although its effect does contribute to the spread of data points in the figure.

We see that in general, the turbulence strength is slightly underestimated, particularly in weaker conditions and the strongest conditions. The scintillation correction process described in section 3.4 reduced the bias on measurements of $r_0$ by approximately 3% and $\tau _0$ by 2%. The residual bias is primarily due to the pixel sampling in the detector being too large to accurately measure the small centroid motions induced by weak turbulence and speckling of the spots in very strong turbulence reducing the signal-to-noise ratio for centroiding. The coherence angle and Rytov variance however depend solely on the vertical turbulence profile calculated from the intensity covariances. The coherence angle displays little evidence of bias, whereas the Rytov variance appears to be slightly underestimated. From analysing the response of the profile reconstruction to a single turbulent layer, there is evidence that the Rytov variance will be underestimated for profiles with strong low-altitude turbulence and the effect of this is greater on this parameter than the coherence angle due to the relative scaling with height.

4.2 36-hour continuous monitoring

Figure 4 shows an example 36-hour sequence of turbulence parameter monitoring at the Roque de Los Muchachos observatory, La Palma, Spain, from the 14th to 16th May 2022. There is a gap from around 5am to 8am UTC on the 15th May due to extremely high levels of wind that temporarily curtailed measurements due to the spot motion as a result of wind shake being larger than the WFS field of view. Although the instrument was still running, measurements taken during this time were not reliable and are therefore excluded from the graph. All stars observed during the day were magnitude zero or brighter in the J band. Due to the large number of individual measurements, for presentation purposes only every second data point has been plotted. We see a strong contrast between night and day conditions in $r_0$ and $\tau _0$ but little difference in $\theta _0$ and $\sigma _R^2$. This suggests that solar heating during the daytime increases the surface layer turbulence strength but does not affect higher altitude turbulence.

 

Fig. 4. 36-hour sequence of turbulence parameters measured at the Observatorio de los Muchachos, La Palma, Spain from the 14th to 16th May 2022. From top to bottom panels show the coherence length, coherence angle, coherence time and Rytov variance. The background colours, dark grey, light grey and white, indicate night, twilight and day respectively. There is a gap around 5am UTC on the first night due to winds in excess of 50 km/hour. Parameters are corrected to zenith and a wavelength of 500 nm.

Download Full Size | PDF

5. Conclusion

We have demonstrated the first continuous vertical turbulence monitor, capable of working in the day and the night without any modifications or downtime. The instrument uses simultaneous measurements of slopes and intensities in a Shack-Hartmann wavefront sensor (WFS) to estimate a four-layer vertical optical turbulence profile, from which the coherence length, angle and Rytov variance are derived, with additional wind speed data from meteorological forecasts used to derive the coherence time. The instrument also implements a method of correcting the ground layer of the measured vertical profile for the effects of scintillation.

The 24-hour Shack-Hartmann Image Motion Monitor (24hSHIMM) is simple and robust, utilising a well-known and trusted Shack-Hartmann wavefront sensor design. Most of the analysis techniques used in this work are well documented in technical literature and have been validated here using detailed Monte-Carlo simulations. As such, the on-sky results have not been compared to alternative turbulence monitoring techniques as it was not possible to co-locate instruments in order to match the local turbulence conditions and there were no other instruments available to compare daytime and twilight measurements. This could be investigated in further experiments.

The 24hSHIMM is a significant development on previous capabilities. This is an important step forward in the field of ground based optical instrumentation as it will enable site selection, performance optimisation and validation through continuous monitoring for astronomical observatories and free-space optical communications. This development delivers the first instrument that can be used at any time of the day and is small, cheap and portable enough to be used in any location, including urban and sub-urban environments. This instrument will be able to build databases of 24-hour turbulence conditions for current and future optical ground station sites in urban areas, while also providing a wealth of data for assimilation into turbulence forecasting models which are critical for the operation of free space optical communications networks, solar and astronomical observatories. This instrument provides the field with the first truly continuous, day and night turbulence monitor for all but the most demanding of applications.