## Abstract

Adaptive optics systems are used to compensate for distortions of the wavefront of light induced by turbulence in the atmosphere. Shack-Hartmann wavefront sensors are used to measure this wavefront distortion before correction. However, in turbulence conditions where strong scintillation (intensity fluctuation) is present, these sensors show considerably worse performance. This is partly because the lenslet arrays of the sensor are designed without regard to scintillation and are not adaptable to changes in turbulence strength. Therefore, we have developed an adaptable Shack-Hartmann wavefront sensor that can flexibly exchange its lenslet array by relying on diffractive lenses displayed on a spatial light modulator instead of utilizing a physical microlens array. This paper presents the principle of the sensor, the design of a deterministic turbulence simulation test-bed, and an analysis how different lenslet arrays perform in scintillation conditions. Our experiments with different turbulence conditions showed that it is advantageous to increase the lenslet size when scintillation is present. The residual phase variance for an array with 24 lenslets was up to 71% lower than for a 112 lenslet array. This shows that the measurement error of focal spots has a strong influence on the performance of a Shack-Hartmann wavefront sensor and that in many cases it makes sense to increase the lenslet size. With our adaptable wavefront sensor such changes in lenslet configurations can be done very quickly and flexibly.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## Corrections

18 November 2020: A typographical correction was made to the author affiliations.

## 1. Introduction

Shack-Hartmann wavefront sensors (SHWFS) are used extensively in astronomical adaptive optics (AO) [1]. However, in more challenging propagation scenarios such as long-range, terrestrial free-space optical communications, where strong scintillation and branch points in the wavefronts are common, these sensors produce increasingly erratic measurements with increasing turbulence strength [2,3]. Scintillation causes light intensity to fluctuate across the pupil. Due to the finite dynamic range of the detector, the spots generated by microlenses receiving high intensities might be saturated, while the spots under dark microlenses might disappear below the noise level. The centroid algorithm therefore has problems finding all spot positions with good accuracy which, in turn, reduces the reliability of the wavefront measurements of the SHWFS. These faulty measurements can then cause stability issues for the feedback loops involved in AO systems [3].

Scintillation occurs on very short scales, both spatially and temporally. In weak-to-moderate scintillation conditions, the intensity correlation width is controlled by the Fresnel scale $\sqrt {L/k}$, where $L$ is the path length and $k$ is the wavenumber. The characteristic correlation time is of the order of $\sqrt {L/k}/V$, where $V$ is the transverse wind speed, for very small receivers, and $D/2V$, where $D$ is the receiver aperture, for receivers significantly larger than the Fresnel scale [4]. For the application considered in this paper, i.e. terrestrial laser-based communications, the wavelength is usually 1550 nm, path lengths are between 1 and 10 km and expected wind speeds are between 1 and 10 m/s. Using these values we can arrive at the typical spatial scales of scintillation of 1 to 5 cm and typical correlation times between 1 and 50 ms.

To improve the performance of a SHWFS, and to gather more light, its exposure time could be increased [5]. However, in closed-loop AO systems this might not be beneficial due to temporal evolution of phase aberrations, which the SHWFS should measure. Another possibility to improve performance is to change the lenslet array of the SHWFS to one with larger individual lenslets that spatially average over the intensity fluctuations, an effect called aperture averaging. Using larger lenslets, however, constitutes a trade off between spatial sampling of the wavefront (higher fitting error) and the ability to locate focal spot positions reliably (lower measurement error thanks to aperture averaging). With larger lenslet apertures, the local $d/r_0$ increases which leads to a breakup of spots (Fig. 1(d), right column displays an enlargement of a broken up spot). The measurement error for photon-noise-limited SHWFSs in seeing-limited operation is proportional to $n_{\mathrm {ph}}^{-1/2} (d/r_0)$, where $n_{\mathrm {ph}}$ is the number of photoelectrons per subaperture of diameter d and per exposure time, and $r_0$ is Fried’s coherence length [6]. In electronic-noise-dominated SHWFSs, the measurement error is proportional to $n_{\mathrm {ph}}^{-1} (d/r_0^2)$. It can be easily shown that by increasing the lenslet size and transitioning from diffraction-limited operation ($d \approx r_0$) to seeing-limited operation ($d > r_0$) the measurement error due to photon noise stays constant and errors due to electronic noise components can actually decrease by an order of magnitude. Both these effects can be explained by the fact that $n_{\mathrm {ph}} \propto d^2$. When scintillation can be neglected, then it is possible to derive an optimal $d$ analytically, by optimizing the fitting, the measurement, and the time-delay error [7].

In a typical SHWFS, lenslet arrays are chosen for the expected turbulence conditions beforehand, and since they are physical elements, they are not easily exchangeable as turbulence changes. Another possibility is to to implement diffractive lenslet arrays with a spatial light modulator (SLM) instead of using physical lenslet arrays. Rha et al. [8–10] and Seifert et al. [11] were the among the first to use liquid crystal devices to display diffractive lenslets as a replacement for physical lenslet arrays in a SHWFS. Rha et al. [8–10] built a small $4\times 4$ lenslet array for a SHWFS which could be changed to a $2\times 2$ configuration when a neutral density (ND) filter reduced the incoming light intensity. They mentioned that lenslets could be locally combined if there was too little light, but, to our knowledge, did not follow up on this idea. Seifert et al. [11] adapted individual lenslets such that focal spots would not be distorted during measurements of optical elements.

In an effort to increase the flexibility and performance of a SHWFS, we built on their ideas and developed a SHWFS with adaptable lenslet arrays displayed on an SLM. This allows our sensor to freely adapt its lenslet array to prevailing turbulence conditions. To this end, we built a turbulence and scintillation test-bed in the laboratory where we can generate deterministic conditions and understand how different lenslet configurations influence the performance of our sensor. We examined two different ways to perform an adaption of the lenslets to scintillation: local and global adaption. Figure 1 shows the lenslet configurations that were used in the left column, the same configurations with an added phase screen ($D/r_0 = 30$) in the middle column and in the right column the resulting inverted camera images when scintillation was applied. The reference lenslet array with 112 lenslets inside the pupil (Fig. 1(a)) shows how scintillation causes spot intensities to be quite low in certain areas. In the local adaption scheme (Fig. 1(b)), four small lenslets are replaced with one large lenslet in these low-light areas and the spot positions can more likely be measured correctly. With global adaption (Figs. 1(c), (d)) all lenslets are enlarged yielding a smaller spot count of 68 or 24 with higher individual spot intensities. Local adaption preserves the spatial sampling of the wavefront in areas with high enough intensities, while at the same time reducing the measurement errors in low-light areas.

In global adaption, the fitting error is increased over the whole pupil, while the measurement errors are reduced for every spot. In theory, local adaption should achieve the best results as it compromises the resolution only in selected areas. As will be made evident in the experiment section, this turned out not to be the case.

This paper starts with the description of the experimental test-bed. It continues with an explanation of the implementation of diffractive lenses, the wavefront reconstruction algorithm and the generation of simulated phase and scintillation screens. Then, the experimental results and the comparison of different lenslet configurations are presented and discussed.

## 2. Methods

In this section, we describe the implementation of our wavefront sensor in detail. We present the laboratory setup, the diffractive lenses, and the wavefront reconstruction method. We also explain our approach to generate turbulent phase and scintillation screens.

#### 2.1 Experimental setup

We first introduce the two devices used for phase and amplitude modulation of the laser, that is, for simulation of phase delay and scintillation, before we give a detailed overview of the experimental setup.

A phase-only SLM is used to modulate the phase of light. The specific model is a PLUTO 1 SLM from HOLOEYE with a resolution of 1920$\times$1080 pixels and a pixel pitch of 8 µm. In order to ensure reliable, linear phase modulation up to 2$\pi$, the SLM was previously characterized and calibrated using an interferometric approach as described by Bergeron et al. [12] and in HOLOEYE’s manuals. By overlaying the phase functions of turbulent phase screens with those of diffractive lenslets, the SLM doubles as both the wavefront aberrator and the lenslet array of the SHWFS that we describe in this paper.

A digital micromirror device (DMD) is responsible for the amplitude modulation of the laser beam and thus for the simulation of scintillation. We use a DLP7000 from Texas Instruments with an array size of 1024x768 micromirrors and a mirror pitch of 13.6 µm that is controlled by a DLP Discovery 4100 Development Platform module. While DMDs are binary amplitude modulators that turn individual pixels on and off by flipping mirrors in or out of the beam path, they are also very fast with frame rates of several tens of kilohertz (32 552 Hz for our device). Due to this speed, they can also display 8-Bit gray level amplitude patterns using a pulse width modulation scheme. During image acquisition, we make sure that the camera exposure time is an integer multiple of the individual frame time of the DMD. Every acquired image then averages over at least 10 DMD frames to ensure that gray values are reproduced reliably. Initially, we intended to use a second SLM and an analyzer as an amplitude modulator. We then realized, that it is not possible to operate a SLM as amplitude-only modulator and that the phase of light is always changed as well, which introduces additional wavefront aberrations.

The laboratory setup depicted in Fig. 2 consists of a DMD and a SLM for scintillation and phase simulation and a SHWFS that is formed by the SLM and a camera that records the focal spots.

We use a 17 mW HeNe laser from Research Electro Optics at 633 nm as the light source of our setup. An ND filter is used to reduce the light intensity of the laser during normal operation. It can also be flipped out of the beam path so that lens back reflections can be used for alignment purposes. The 1 mm laser beam is collimated and expanded by a factor of 10 and a 50 µm pin hole is inserted in the focal plane of the first lens to suppress higher order modes of the laser.

The beam then hits the DMD which imprints the intensity profile of the scintillation screens in the beam. The individual square micromirrors in most DMDs flip around their diagonal and not around their vertical axis. Therefore, the DMD has to be rotated by 45° if the beam is to stay in the horizontal plane. This setup constraint is of no concern for our experiment, however, as only the central circular and thus rotationally symmetric region of the DMD is used for applying scintillation to the beam.

The DMD plane is then conjugated onto the SLM plane with a 4-f system. An iris is inserted in the focal plane to filter out higher orders stemming from the grid of micromirrors on the DMD. The iris is opened to half the distance between the zeroth and the first order focus as we want to filter as little information as possible from the zeroth order. An alignment camera is used on one side of the 50:50 beam splitter to fine-tune the conjugation of the DMD to the SLM plane. When a checkerboard pattern is displayed on the DMD, the alignment camera is moved until the pattern is in focus. Then, the distance of the SLM to the center of the beam splitter is adjusted to be the same as the distance of the camera to that point.

The SLM displays a diffractive lenslet array that focuses light onto a camera that is located behind the beam splitter. The arrays all have a focal length of 122.7 mm in our current setup. All lenslet parameters can be freely changed in the control software as the lens patterns are generated on the fly. A linear polarizer that is aligned to the horizontal axis of the SLM is added before the beam splitter to ensure phase-only modulation of the SLM. Due to its large sensor, we use a Bonito CL-400B camera as our detector. The camera also has a very high framerate, but since we integrate over several DMD frames this feature is of no concern in our setup. According to the turbulent phase screens that are displayed, the focal spots will move in the detector plane allowing a reconstruction of the wavefront aberration.

In a typical measurement series, we keep one scintillation screen constantly displayed on the DMD and acquire one image for each of the 100 phase screens. For each of these images, an average of three images is taken to reduce measurement noise. This measurement is repeated for each of 50 scintillation screens at the chosen turbulence strength. When $50 \times 100$ images are taken and one measurement series is completed, it is repeated for different lenslet configurations. Before each series of 100 phase screens, a dark frame and a reference image without added turbulence phase or scintillation screen is acquired.

The local adaption to scintillation (Fig. 1(b)) which combines lenslets only in low-light areas was done manually. We displayed each scintillation screen and edited the lenslet configurations with the GUI of our LabVIEW control program. We set a noise threshold similar to the one used in the centroid algorithm (mean plus three times the standard deviation of several noise samples) and adapted the lenslets so that their focal spots were above the noise level. The resulting lenslet configuration was then saved and read during measurements. In the future this adaption will be automated.

We keep the exposure time constant throughout the experiment. To determine the exposure time, we calculate the average gray value of all scintillation screens, which corresponds to the average light intensity of the screens, and display that gray level on the DMD. Then we set the exposure time so that the focal spots produced by the non-adapted 112 lenslet array (Fig. 1(a)) do not saturate the sensor (70 ms for $D/r_0 = 10$, 90 ms for $D/r_0 = 20$, 150 ms for $D/r_0 = 30$; see Section 2.4 for an explanation of the simulated turbulence conditions). When the actual scintillation screens with differing light intensities across the pupil are displayed and the fixed exposure time is used, some focal spots will be overexposed due to the finite dynamic range of the camera sensor. However, allowing saturation actually increases the dynamic range of a SHWFS essentially improving the wavefront reconstruction performance as was shown in simulations by Plett et al. [5].

#### 2.2 SLM as diffractive lenslet array and wavefront aberrator

In order to display lenses on the phase screen of the SLM, we replicate the phase function of a spherical glass lens on the SLM. When a plane wave hits a plano-convex lens, it is converted into a spherical wave that converges to a focal point. Depending on the varying thickness of glass, a lens adds a local phase to the wavefront of light. For a spherical lens, the phase profile relative to the on-axis phase is a function of the thickness profile of the lens and is described by [13,14]

with $\lambda$ as the wavelength of light, $f$ the focal length of the lens, $r$ the radial coordinate and $x$ and $y$ the corresponding Cartesian coordinates for an easier description in pixel coordinates.There is also a minimum focal length for displaying individual diffractive lenslets that is defined by [11]

This depends on the number of pixels on one side of a square lens $n$, the pixel pitch of the SLM $p$ and the wavelength $\lambda$. Lenslets with shorter focal lengths than this minimum are considered under-sampled and should be avoided due to additional higher diffraction orders that appear in the same focal plane. For our arrays with lenslet side lengths of 170, 105 and 85 pixels, the minimum focal lengths are 17.2 mm, 10.6 mm and 8.6 mm respectively. The minimum focal length in our setup is set by the distance between the reflective SLM and the camera which are located on two sides of a beam splitter. With a resulting focal length of 122.7 mm we are far above the theoretical minimum focal length. Shorter focal lengths could also be used by relaying the focal plane with another lens.Our HOLOEYE SLM was calibrated for the specific experimental wavelength which means that there are 256 phase levels between 0 and 2$\pi$. However, fewer phase levels are usually available due to, for example, phase flicker, voltage discretization, temperature variations or the desire to use multiple wavelengths with one SLM without recalibration. Some SLMs are also not user configurable and have a factory calibration for a specific wavelength range, like Hamamatsu’s SLMs. However, also SLMs with a lower actual phase resolution can be used for displaying lenses. Theoretical diffraction efficiencies of up to 99% can be achieved for diffractive lenses with only 16 phase levels as demonstrated by O’Shea et al. [13]. This is in line with our measurements which only showed a small increase of around 5% in diffraction efficiency when displaying lenses with 16 or 256 phase levels. Thus, the available phase resolution is usually not a limiting factor when displaying diffractive lenses.

In order to display a lens phase profile as described in Eq. (1) on a SLM, it is wrapped at a phase of 2$\pi$ and converted to gray values according to the phase modulation function of the SLM. Such a resulting lens is often called a phase Fresnel lens. The combination of multiple such lenses forms the lenslet arrays used in our experiments. Our arrays use 12, 10 and 6 lenses on a side which results in 112, 68 and 24 lenslets inside the circular pupil. When the SLM is also used as a wavefront aberrator as in our experiments, simulated turbulent phase screens are simply added to the phase profile of the lenslet arrays and wrapped and converted accordingly. The resulting lenslet array overlaid with phase screens can be seen for the turbulence case of $D/r_0=30$ in the middle column of Fig. 1.

#### 2.3 Wavefront reconstruction algorithm

The modal Zernike wavefront reconstruction algorithm that we used in this work is based on the freely available MATLAB-based Shack-Hartman wavefront sensor toolbox [15] developed by Antonello [16]. It was heavily modified to suit our needs of locally-adapted lenslet arrays where differently sized subapertures are present.

First, the displacements of focal spots due to a phase distortion are computed by a thresholding centroid algorithm [6,17,18]. A previously recorded dark frame is subtracted from each image and a threshold of the mean plus three times the standard deviation of several noise samples in each image is applied to reduce the influence of noise on the centroid search. The displacements in $x-$ and $y-$direction $\Delta x$ and $\Delta y$ are then calculated in comparison to the reference positions. From this, the average slopes $s_x$ and $s_y$ of the wavefront over a lenslet $l$ with focal length $f$ are determined from the small angle approximation by [6]

The vector $\mathbf {s}$ contains the measured wavefront slopes $s_{x_l}$ and $s_{y_l}$ for each lenslet, the matrix $\mathbf {E}$ consists of the average derivatives of Zernike polynomials $\left . \frac { \partial Z_j} { \partial x } \right \rvert _{_l}$ and $\left . \frac { \partial Z_j} { \partial y } \right \rvert _{_l}$ evaluated over the subaperture area of lenslet $l$ and $\mathbf {a}$ is the unknown vector of Zernike coefficients $a_j$. Each row in $\mathbf {s}$ and $\mathbf {E}$ corresponds to one specific lenslet that is numbered from 1 to $L$.

In order to solve the system of linear equations in (5), Singular Value Decomposition is used to calculate the generalized inverse $\mathbf {E^+}$ of $\mathbf {E}$. The Zernike coefficient vector $\mathbf {a}$ is then found by $\mathbf {a = E^+s}$ [19].

When four lenslets are combined to one larger lenslet in local adaption, the total number of lenslets reduces by three. Thus, the number of rows in the matrix $\mathbf {E}$ and the slope vector $\mathbf {s}$ reduces by $3\times 2 = 6$ ($\times 2$ as there is one entry each for $x-$ and $y-$direction). For example, if there are 10 adapted lenslets in a 112 lenslet array, the number of rows of $\mathbf {E}$ would be $(112-10\times 3)\times 2 = 164$ instead of $112\times 2=224$. The average derivatives of the Zernike polynomials for an adapted lenslet in $\mathbf {E}$ are then taken over the corresponding increased subaperture area. Therefore, each lenslet configuration has its own matrix with a smaller row dimension in local adaption and the lenslets are numbered differently for each configuration. It is possible to precompute the generalized inverse $\mathbf {E^+}$ for each possible lenslet configuration and store these inverses for efficient computation of the coefficients of the Zernike polynomials.

#### 2.4 Generation of phase and scintillation screens

We investigated three turbulence levels in simulations: $D/r_0 = 10$, 20 and 30 where $D$ is the receiver aperture diameter and $r_0$ is Fried’s coherence length. We assumed $D = 10\;\textrm{cm}$ and $\lambda = 633\;\textrm{nm}$, to generate 100 pupil-plane phase screens using a covariance matrix of Zernike modes in Kolmogorov turbulence [20]. We diagonalized this matrix using singular value decomposition, leading effectively to a Karhunen-Loéve projection of the Zernike modes. This is the most accurate method of simulating atmospheric phase screens, albeit also the most time-consuming [21]. We chose the number of Zernike orders in the simulations to be 60. This translates to 1890 actual Zernike polynomials in the covariance matrix and a very accurate representation of the turbulent phase. The simulations generated phase screens on a grid of $1080\times 1080$ pixels to match the resolution of the SLM. We validated the average, long-exposure point-spread function against the theory of Fried [22] and found a very good match.

We generated the scintillation patterns using wave-optics simulations [23]. In order to accurately represent scintillation, several phase screens along the propagation path and Fresnel propagation between them have to be implemented. Here, the aforementioned, Zernike-based method is not viable as it is quite slow. Thus, independent phase screens were generated using the fast Fourier transform (FFT) based method of McGlamery [24], whereby an array of random numbers is filtered according to the Kolmogorov spectrum and then inverse Fourier-transformed. This method suffers from inadequate representation of low-spatial-frequency aberration, and therefore, our simulation implemented subharmonics correction [25] with eight levels of subharmonics. This yielded good accuracy when simulated and theoretical phase structure functions were compared. A collimated Gaussian beam of wavelength 633 nm and beam radius $w_0= 4\;\textrm{cm}$ was numerically propagated from a transmitter of diameter 20 cm through a 1.2 km path. For $D/r_0 = 10$, the average refractive index structure parameter $C_n^2$ along the path was set to $6 \times 10^{-14}\;\textrm {m}^{-\frac {2}{3}}$, for $D/r_0 = 20$ to $1.2 \times 10^{-13}\;\textrm {m}^{-\frac {2}{3}}$ and for $D/r_0 = 30$ to $1.67 \times 10^{-13}\;\textrm {m}^{-\frac {2}{3}}$. The value of $r_0$ for each screen was optimized to match the global, plane-wave $r_0$ as well as the theoretical on-axis scintillation index [23,26]. As recommended by Schmidt [23], the strength of each phase screen was optimized in order to match analytical quantities, which assume propagation through a continuous medium as opposed to the approximate, discrete representation that the split-step method uses. This must be done as separating the medium into discrete slabs can result in discrete versions of some quantities like Fried parameter, $r_0$, or cause the isoplanatic angle to deviate from the assumed values. In this work, we decided to match continuous and discrete values of two parameters: global, plane wave $r_0$ and on-axis scintillation index of a collimated Gaussian beam. This does not mean that the simulation is limited to plane-wave propagation. On the contrary, by optimizing the plane wave $r_0$ and the Gaussian beam scintillation index expressions, we have a sound propagation setup for collimated beams (the beams in these simulations do not expand significantly to be approximated by point sources).

Naturally, in a real free-space optical communications system a different wavelength, for example 1550 nm, would have been employed, but here we chose the wavelength to match the experimental wavelength of 633 nm. The propagation grid was $2048\times 2048$ pixels from which, on the receiver side, we cut out a central $768\times 768$ portion to match the resolution of the DMD. The physical grid size in meters was set to four times the theoretical value of the turbulence-induced beam spread [26] in order to contain the beam throughout the entire propagation path. The number of screens along the path was optimized so that the Rytov variance does not increase by more than 10% between two successive screens [27]. The resulting numbers were 14, 15 and 19 screens for $D/r_0 = 10$, 20 and 30, respectively. We implemented beam propagation between the phase screens using the Fresnel integral and the FFT [23]. The influence of the Gibbs phenomenon ("ringing" [23]) was alleviated using apodization with the Butterworth function of order 8. The on-axis scintillation indices measured in the plane of the receiver for $D/r_0=10$, $20$ and $30$ were 3.45, 2.36 and 2.16, respectively (the latter two cases are in the saturation regime), but speckles had significantly finer structure in the stronger turbulence simulations.

## 3. Experiments

This section is devoted to the results of our experiments on the comparison of different lenslet arrays in scintillation conditions with the developed adaptable SHWFS.

To illustrate how our system reconstructs a wavefront, an example of reconstruction with and without applied scintillation is given in Fig. 3 for a turbulence strength of $D/r_0 = 10$. The left column represents the input phase screen, the middle column the reconstruction and the right column the residual wavefront after correction, i.e. subtraction of the reconstructed wavefront from the ground truth phase screen. A comparison of the two lenslet configurations with 12 and 6 lenslets across shows that the configuration with more lenslets performs better when no scintillation is present, however, the one with fewer lenslets achieves better results when scintillation is added. Of course, this is not the case for every phase or scintillation screen, which is why we have to quantify the average reconstruction performance for many phase and scintillation screens.

The residual phase variance after the perfect removal of a number of Zernike modes is used as a metric to compare the different sensor configurations and thus the reconstruction performance of the sensor. As the induced wavefront aberrations are known, we can compare the reconstructed Zernike coefficients $a_{j\, \mathrm {rec}}$ with the ground truth coefficients $a_{j\, \mathrm {truth}}$ of the phase screens. The ensemble average over all measurements in one series (5000) yields the residual phase variance $\sigma _\phi ^2$ as

Figure 4 shows the experimental results of wavefront reconstructions with different lenslet arrays in three different turbulence conditions ($D/r_0 = 10$, 20 and 30). The normalized residual phase variance $\sigma _{\phi ,\mathrm {norm}}^2$ is plotted against the number of removed Zernike modes $N$. The plots start with the piston-removed phase variance (Zernike mode 1 in Noll’s Zernike numbering system [20]) and the ticks on the $x$-axis correspond to the radial Zernike orders ("rows" in the Zernike pyramid). As the removal of the first two Zernike modes (tilts) already significantly reduces the phase variance, we limit the $y$ axis in all cases to $\sigma _{\phi ,\mathrm {norm}}^2 = 0.35 \, \mathrm {rad}^2$ as it would otherwise be hard to discern the features of the graphs. The theoretical limit for the residual phase variance is the perfect sensing and correction of Zernike modes of a random wavefront [20] that is plotted in black in the sub figures. The ticks on the $y$-axis denote the minimum residual phase variance $\sigma _{\phi ,\mathrm {min}}^2$ for each reconstruction.

We only consider up to 21 removed Zernike modes from the wavefront. When too many modes are reconstructed and removed, the residual phase variance actually increases due to errors in the reconstruction as was shown by Takato et al. [28]. In the case without scintillation, it is actually possible to reconstruct up to 45 modes without a subsequent increase in reconstruction error. However, as soon as scintillation is involved and the measurement error increases, the maximum number of reconstructable modes reduces, which also coincides with Takato’s calculations. It would also have been possible to reconstruct more modes than 21 for some cases, but the gain would be minimal as the reduction in residual phase variance is quite flat after 21 modes. For the 24 lenslet array, however, we stopped the reconstruction already after removing 15 modes, as afterwards the reconstruction errors become larger and the residual phase variance starts to grow. This is due to the higher fitting error for a small number of lenslets. Therefore, we only show the minimum residual phase variances for the respective configurations up to 21 modes.

Without scintillation, the wavefront reconstruction with 112 lenslets was very good and came close to the theoretical limit (light green curve in Fig. 4(a)-(c)). In previous work, we found that the more subapertures were used in our sensor, the better the performance of our sensor was [29], which is to be expected as the fitting error gets smaller. As we mentioned in the previous paragraph, we can reconstruct up to 45 modes reliably and remove them from the wavefront when no scintillation is present (modes higher than 21 are not shown in the plots in Fig. 4). The minimum residual phase variance in this case would be 32.7%, 37.8% and 26.7% better than with only 21 modes for the three turbulence cases $D/r_0 = 10$, 20 and 30.

In comparison to the case without scintillation, when scintillation is involved, the minimum residual phase variance $\sigma _{\phi ,\mathrm {min}}^2$ for the unadapted 112 lenlet array increases by a factor of 5 for $D/r_0 = 10$ (0.028 rad$^2$ to 0.140 rad$^2$), a factor of 8.3 for $D/r_0 = 20$ (0.025 rad$^2$ to 0.208 rad$^2$) and a factor of 9.1 for $D/r_0 = 30$ (0.034 rad$^2$ to 0.310 rad$^2$). The previous advantage of having a small fitting error vanishes. When 45 modes were reconstructed, the decrease in residual phase variance compared to 21 reconstructed modes is only 2.2%, 2.6% and 2.9%.

When we adapted the lenslet arrays and, through aperture averaging, reduced the influence of scintillation, the residual phase variance after correction became smaller and the performance of the sensor improved. Beside the curves in Fig. 4, an overview of the percentage decrease in minimum residual phase variance in comparison to the unadapted 112 lenslet array in scintillation conditions, is presented in Table 1 for each lenslet configuration. The locally-adapted lenslet array and the 68 lenslet array behaved similarly and reduced the residual phase variance. It is notable that the curves for $D/r_0 = 10$ and 20 lie almost on top of each other. We think this is mostly coincidental as the non-averaged data does not show such a similarity for the individual measurements. For local adaption, we had on average 10.3, 14.2 and 13.2 adapted, large lenslets in $D/r_0 = 10$, 20 and 30, respectively. The lowest residual phase variances in scintillation were achieved with the 24 lenslet array in all turbulence strengths for 15 removed modes (red curves in Fig. 4). We obtained a relative decrease in minimum residual phase variance of $-56.6\%$, $-64.0\%$ and $-71.3\%$ for $D/r_0 = 10$, 20 and 30, respectively.

## 4. Discussion

We attribute the worsening of the reconstruction performance with and without scintillation mainly to an increase in the centroid measurement error as underexposed focal spots due to scintillation are hard to locate for the centroid algorithm. Overexposed focal spots due to scintillation reduce the accuracy of finding the spot positions, but an approximate position can still be found reliably. As Plett et al. [5] have found, overexposure therefore can actually increase the dynamic range of a SHWFS in scintillation. The higher the turbulence strength, the worse the minimum residual phase variance becomes when scintillation is present (0.14 rad$^2$, 0.208 rad$^2$ and 0.31 rad$^2$ for $D/r_0 = 10$, 20 and 30). Without scintillation, this increase is not as significant. Wavefront reconstruction for $D/r_0 = 20$ is in fact slightly better ($\smash {\sigma _{\phi ,\mathrm {min}}^2 = 0.025\,\mathrm {rad}^2}$) than the one at $D/r_0 = 10$ ($\smash {\sigma _{\phi ,\mathrm {min}}^2 = 0.028\, \mathrm {rad}^2}$) and for $D/r_0 = 30$ it is only slightly worse ($\sigma _{\phi ,\mathrm {min}}^2 = 0.034\, \mathrm {rad}^2$). We suspect that the worsening performance in strong turbulence with scintillation is due to the fact that phase aberrations cause a distortion and broadening of the focal spots (breakup of spots, see Fig. 1(d), right column, enlargement). These turbulence-limited spots have a non-centrosymmetric shape and low peak intensities due to the break down of the focal spot. Together with scintillation, underexposed focal spots could be even harder to locate and therefore the centroid measurement error increases even more when scintillation is present. Furthermore, we found that due to the limited diffraction efficiency of the lenslets, there was always some undiffracted light around the produced focal spots, which adds noise to the centroid search.

Another factor affecting the strongly decreased reconstruction performance with scintillation might be that we used a Gaussian beam as a light source. The Gaussian intensity profile causes focal spots on the borders of the pupil to have lower intensities than the ones in the middle. When adding scintillation, these spots are more prone to fall below the noise level. We noticed this during the manual selection of locally-adapted lenslets where the outer lenslets would be adapted more often than the inner ones. Also, the simulated scintillation patterns come from a Gaussian beam, meaning that there would more likely be higher intensities in the middle of the pupil. In the future, we would flatten the intensity distribution of the light source with an intensity mask that could be applied on the DMD together with the scintillation screens.

Before performing the experiments, we assumed that locally-adapted lenslet arrays would be the best approach to counter scintillation as they would keep the spatial sampling high in areas with high intensity while sacrificing resolution only where necessary. However, in our experiments this did not turn out to be the case: a simple global change to a 68 lenslet array yielded the same or better reconstruction performance as the locally-adapted lenslet array. This might be due to the fact that we did not adapt enough lenslets overall or that our adaption was limited to combine only 4, not 9 or more lenslets. Another factor that we did not consider in the experiment was whether the non-uniform lenslet geometries in local adaption actually affect the performance of the reconstruction algorithm. The non-uniform distribution of measurement points could lead to some modes being sensed worse than for a uniform geometry. For example vertical coma could be mistaken for secondary vertical astigmatism or defocus if lenslets were grouped densely only in the top part of the pupil. In a future experiment, we want to investigate how much of an impact non-uniform arrays have on reconstruction. Furthermore, we want to test other adaption schemes such as non-square or L-shaped lenslets as sometimes we could not adapt all lenslets due to the geometry of the array. We will also write an automated algorithm to adapt lenslets locally.

There are further drawbacks that make local adaption more of a proof-of-principle approach that is not practical in real systems. The most severe disadvantage is the speed of SLMs. Since scintillation changes on timescales of around 1000 Hz, the local lenslet adaption would have to be faster than that. In the first measurement, the system would have to find out which lenses are underexposed before it could actually change the lenslet configuration. However, as liquid-crystal-based SLMs are limited in speed by the viscosity of the crystals, their update rates are in the range of tens of Hz. Ferroelectric SLMs are faster with around 1000 Hz switching rate and could theoretically be used, but since they are binary phase modulators the diffraction efficiency of lenses would be much lower and even less light would be available for sensing the wavefront. Overall, we think that with current technology, local adaption of the lenslets is not a feasible approach.

Global adaption of lenslet arrays to changing turbulence conditions is a more viable method for practical AO. As we could see in our measurements, the performance of the adaptable SHWFS in scintillation improved when the lenslet size was increased due to aperture averaging, even though less subapertures were available for wavefront reconstruction. The advantage of being able to correct many modes with many small lenslets was undermined by the strong increase in measurement error that scintillation introduced. The measurement error thus has a higher influence on the performance characteristics of our SHWFS than the fitting error.

There is one more source of error, which "although present in our measurements" is not analyzed in detail here. This is the aliasing error, which arises because the spatial-frequency content of the Kolmogorov phase aberration is not band limited, but a SHWFS will sample it with a finite spatial period. An elegant way to reduce aliasing is to position a field stop at a focal plane before the SHWFS. This stop acts as a low-pass filter on the phase and will reduce the high-spatial-frequency content seen by the sensor [30]. The size of this spatial filter, which, when paired with the SHWFS proposed in the current paper should be square, depends on the geometry of the SHWFS: Big lenslets require a small stop, while small lenslets allow for a larger stop. In an adaptable SHWFS, the size of the field stop should be variable to account for the changing size of the subapertures. In the case of local adaption of the lenslet array, the size of the field stop cannot be set properly for both lenslet sizes, which are present at the same time. With the proposed global adaption of the lenslet array, a variable field stop could be realized with an electronically controlled iris. On the other hand, the benefits of spatially filtered SHWFS are seen only for very high Strehl ratios ("extreme adaptive optics"). For the scenarios discussed in this paper, which assume that the phase aberration is large and produce a large focal-plane spot, and where scintillation triggers a switch to larger lenslets in turn requiring a smaller field stop, a spatially filtered SHWFS might actually be detrimental in the same way scintillation is: the loss of light due to a small field stop can result in significant fluctuations in intensity across the pupil [30], which is precisely the effect we want to compensate for.

The underlying mechanism which improves SHWFS performance when a global (and local) adaption approach is used, is aperture averaging: larger lenslets average over more scintillation patterns than smaller lenslets. It would be useful to have in advance some guidance as to choosing the optimal lenslet size. It is known that aperture averaging becomes efficient only when d is increased beyond the dominant scintillation scale and this efficiency only persists until around ten times this scale [31]. Beyond this value, the aperture averaging curve flattens out. The dominant scintillation scale can be the Fresnel scale $\sqrt {L/k}$, the inner scale $l_0$ or $r_0$, depending on the turbulence regime. Without relying on any model, we have directly computed the intensity correlation width $\rho _c$ from the $1/e^2$ point of the normalized covariance function of the simulated scintillation screens. The values of $\rho _c$ were 10, 4.4 and 5.4 mm for $D/r_0$ = 10, 20 and 30, respectively (for reference note that $\sqrt {L/k}$ was constant for all simulations and equal to 17 mm). This means that for the original configuration with 112 lenslets (12 lenslets across) $d/\rho _c = 0.7$, 2.1 and 2 for the three turbulence strengths. After adapting the sensor to a pattern with 24 lenslets (6 lenslets across) this ratio doubled to 1.4, 4.2 and 4, respectively. This means that the lenslet diameter entered the region of most efficient aperture averaging and explains the significant performance gains listed in Table 1. We also observed that in our simulations $\rho _c$ was closer to $r_0$ than to $\sqrt {L/k}$. We recommend a global adaption scheme with $d$ in the range 1 to 10 times $\sqrt {L/k}$ in the weak fluctuation regime and in the range 1 to 10 times $r_0$ for longer propagation paths and/or high $C_n^2$.

Global adaption approaches are already used in some real AO systems where several microlens arrays are mounted on a motorized stage [10,32]. However, we think that SLMs offer more flexibility than changing physical microlens arrays. Depending on the application, it would be possible to change the lenslet configuration very quickly to respond to, for example, cloud cover and then back. There would also be no alignment issues as the SLM position stays fixed while the lenslets in a motorized system might move, cause vibrations and be more difficult to align to a detector [10]. The lenslet arrays could also be shaped to any desired form, for example, to match the actuator layout of a deformable mirror. Furthermore, the SLM could also be used to correct system aberrations at the same time as acting as the lenslet array.

However, in light-starved conditions, the light loss when using SLMs as lenslet arrays is a big drawback. There are multiple sources of light loss in our system: the required polarizer for the SLM, the limited diffraction efficiency of the lenslets and the beam splitter that is needed with a reflective SLM. The beam splitter is not necessary when utilizing a transmissive SLM, however, which also reduces the optical complexity to be similar to that of a system with traditional lenslet arrays. Also, diffractive lenslets only work for a specific wavelength which impairs the broadband applicability of our system. This might not be a problem in free-space optical communications systems which are optimized for a design wavelength, but in astronomical AO a narrow-band filter would have to be used which again causes light to be lost.

There is an important effect that the experiments described in this paper do not capture and that becomes very important in propagation over extended volumes of turbulence: Branch points occur in an optical wavefront when it is heavily distorted or the interference is such that zeros of intensity are seen in the plane of the receiver [33]. At these nulls of intensity, the phase of the optical wavefront is undefined and it is at these positions that branch points are formed. Branch points occur in pairs (which are of opposite rotation or sign) and are connected by wave dislocations called branch cuts. Along these branch cuts the phase of the wavefront undergoes a $2\pi$ jump. The branch point is at the origin of the $2\pi$ discontinuity that causes these wave dislocations. Due to branch points, the phase of the optical wavefront at the receiver is not continuous and the least-squares reconstructor used in this paper is unable to sense this hidden phase component. In order to combat the occurrence of branch points, they must first be detected and the associated branch cuts moved into positions inside the pupil plane where they pose least threat to wavefront reconstruction. Proposed methods for branch point detection include the vortex potential and the contour sum algorithms [33,34]. Regarding the subsequent wavefront reconstruction step, some success has been demonstrated with, for example, the complex exponential reconstructor [35]. In order to obtain a realistic estimate of the performance of an adaptable SHWFS in "deep" turbulence we will in the future use advanced, branch-point-tolerant reconstructors.

To avoid the reconstruction problem associated with gradient wavefront sensors in strong scintillation, wavefront sensors (WFSs) based on self-referencing interferometers (SRIs) have been proposed as a promising alternative. SRIs measure the local phase differences from differential intensity measurements making them potentially much less sensitive to the presence of branch points (by avoiding the gradient reconstruction problem), and scintillation across the aperture. A laboratory-based AO system with an SRI sensor has been shown to work well in the strong turbulence regime up to a Rytov variance of 3.3 [36]. Also, the shearing interferometer has been shown to perform better than SHWFSs under strong scintillation conditions [37]. A more recent development has shown that digital holography is able to resolve branch points associated with strong scintillation, because it provides an estimate of the complex-optical field and therefore gives access to the least-squares, and the hidden phase components [38]. The two other WFS common in astronomy, the curvature and the pyramid WFS [1], have not yet to our knowledge been tested for their robustness against scintillation. It can be expected though, that in the pyramid WFS, digitally binning the pixels in each of the four quadrants of the detector would lead to the same improvements as cited here, due to the same underlying process of aperture averaging.

## 5. Conclusion

The wavefront reconstruction of Shack-Hartmann wavefront sensors is increasingly unreliable in turbulence conditions where strong scintillation is present. This is due to the fixed microlens geometry that is typically used in SHWFSs which cannot be updated according to the prevailing turbulence conditions. This paper reviewed our approach to build a flexible, adaptable SHWFS based on diffractive lenslet arrays and the turbulence test-bed used to test the wavefront reconstruction performance of different lenslet configurations in scintillation conditions.

In experiments, we verified that in presence of scintillation, the wavefront reconstruction performance of our SHWFS worsened and that we could counter this effect by increasing the lenslet size. We initially assumed that a locally-adapted and thus mixed-size lenslet array would perform best, but found that global adaption is the more feasible and better performing approach. The best reconstruction results were obtained for the array with 24 lenslets with a reduction in residual phase variance of up to 71% compared to a 112 lenlet array in scintillation conditions.

In future measurements, we want to test our sensor under real turbulence conditions with the help of the 8 km free-space laser link that is currently being built between our institute in Ettlingen and the Karlsruhe Institute of Technology. There, we could verify the results from our deterministic turbulence simulator with real-world turbulence. Also, we want to explore new adaption and reconstruction approaches, such as noise-weighted reconstruction, to allow our sensor to be involved in a real-time AO system.

## Funding

Office of Naval Research Global (N62909-17-1-2037); Wehrtechnische Dienstelle 91 (LIMIS-UM II).

## Acknowledgments

We want to acknowledge help from the whole Adaptive Optics Group at Fraunhofer IOSB. We are especially indebted to Prof. Andrew Lambert for the helpful and lively discussions about turbulence theory, setup considerations and many other optics-related topics.

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## References

**1. **R. K. Tyson, * Principles of Adaptive Optics* (CRC, 2011), 3rd ed.

**2. **C. A. Primmerman, T. R. Price, R. A. Humphreys, B. G. Zollars, H. T. Barclay, and J. Herrmann, “Atmospheric-compensation experiments in strong-scintillation conditions,” Appl. Opt. **34**(12), 2081–2088 (1995). [CrossRef]

**3. **G. Marchi and C. Scheifling, “Turbulence compensation for laser communication, imaging on horizontal path and non-natural turbulence using direct and iterative procedures,” Proc. SPIE **8520**, 85200C (2012). [CrossRef]

**4. **H. T. Yura, “Optical downlink propagation from space-to-earth: aperture-averaged power fluctuations, temporal covariance and power spectrum,” Opt. Express **26**(21), 26787–26809 (2018). [CrossRef]

**5. **M. L. Plett, P. R. Barbier, D. W. Rush, P. Polak-Dingels, and B. M. Levine, “Measurement error of a Shack-Hartmann wavefront sensor in strong scintillation conditions,” Proc. SPIE **3433**, 211–220 (1998). [CrossRef]

**6. **F. Roddier, * Adaptive Optics in Astronomy* (Cambridge University, 1999).

**7. **C. Baranec and R. Dekany, “Study of a MEMS-based Shack-Hartmann wavefront sensor with adjustable pupil sampling for astronomical adaptive optics,” Appl. Opt. **47**(28), 5155–5162 (2008). [CrossRef]

**8. **J. Rha and M. K. Giles, “Implementation of an adaptive Shack-Hartmann sensor using a phase-modulated liquid crystal spatial light modulator,” Proc. SPIE **4493**, 80–87 (2002). [CrossRef]

**9. **J. Rha, “Wave Front Sensing and Reconstruction Using a Twisted Nematic Liquid Crystal Device,” Ph.D. thesis, New Mexico State University (2003).

**10. **J. Rha, D. G. Voelz, and M. K. Giles, “Reconfigurable Shack-Hartmann wavefront sensor,” Opt. Eng. **43**(1), 251–256 (2004). [CrossRef]

**11. **L. Seifert, J. Liesener, and H. J. Tiziani, “The adaptive Shack-Hartmann sensor,” Opt. Commun. **216**(4-6), 313–319 (2003). [CrossRef]

**12. **A. Bergeron, J. Gauvin, F. Gagnon, D. Gingras, H. H. Arsenault, and M. Doucet, “Phase calibration and applications of a liquid-crystal spatial light modulator,” Appl. Opt. **34**(23), 5133–5139 (1995). [CrossRef]

**13. **D. C. O’Shea, T. J. Suleski, A. D. Kathman, and D. W. Prather, * Diffractive Optics: Design, Fabrication, and Test* (SPIE, 2003).

**14. **D. Preece, R. Bowman, G. Gibson, and M. Padgett, “A comprehensive software suite for optical trapping and manipulation,” Proc. SPIE **7400**, 74001N (2009). [CrossRef]

**15. **J. Antonello, “Modal-Shack-Hartmann-Wavefront-Sensor toolbox,” https://github.com/jacopoantonello/mshwfs (2017).

**16. **J. Antonello, “Optimisation-based wavefront sensorless adaptive optics for microscopy,” Ph.D. thesis, Delft University of Technology (2014).

**17. **F. Kong, M. C. Polo, and A. Lambert, “Centroid estimation for a Shack–Hartmann wavefront sensor based on stream processing,” Appl. Opt. **56**(23), 6466–6475 (2017). [CrossRef]

**18. **S. Thomas, “Optimized centroid computing in a Shack-Hartmann sensor,” Proc. SPIE **5490**, 1238 (2004). [CrossRef]

**19. **G.-M. Dai, “Modal wave-front reconstruction with Zernike polynomials and Karhunen–Loève functions,” J. Opt. Soc. Am. A **13**(6), 1218–1225 (1996). [CrossRef]

**20. **R. J. Noll, “Zernike polynomials and atmospheric turbulence,” J. Opt. Soc. Am. **66**(3), 207–211 (1976). [CrossRef]

**21. **N. A. Roddier, “Atmospheric wavefront simulation using Zernike polynomials,” Opt. Eng. **29**(10), 1174–1180 (1990). [CrossRef]

**22. **D. L. Fried, “Optical Resolution Through a Randomly Inhomogeneous Medium for Very Long and Very Short Exposures,” J. Opt. Soc. Am. **56**(10), 1372–1379 (1966). [CrossRef]

**23. **J. D. Schmidt, * Numerical Simulation of Optical Wave Propagation* (SPIE, 2010).

**24. **B. L. McGlamery, “Computer Simulation Studies Of Compensation Of Turbulence Degraded Images,” Proc. SPIE **0074**, 225–233 (1976). [CrossRef]

**25. **R. G. Lane, A. Glindemann, and J. C. Dainty, “Simulation of a Kolmogorov phase screen,” Waves in Random Media **2**(3), 209–224 (1992). [CrossRef]

**26. **L. C. Andrews, * Field Guide to Atmospheric Optics* (SPIE, 2019), 2nd ed.

**27. **J. M. Martin and S. M. Flatté, “Intensity images and statistics from numerical simulation of wave propagation in 3-D random media,” Appl. Opt. **27**(11), 2111–2126 (1988). [CrossRef]

**28. **N. Takato, M. Iye, and I. Yamaguchi, “Wavefront reconstruction error of Shack-Hartmann wavefront sensors,” Publ. Astron. Soc. Pac. **106**, 182–188 (1994). [CrossRef]

**29. **D. Lechner, A. Zepp, and S. Gładysz, “Development of a Shack-Hartmann Sensor Based on Adaptable Diffractive Lens Arrays for Reduction of Scintillation Effects,” in Imaging and Applied Optics Conference, (2019).

**30. **L. A. Poyneer and B. Macintosh, “Spatially filtered wave-front sensor for high-order adaptive optics,” J. Opt. Soc. Am. A **21**(5), 810–819 (2004). [CrossRef]

**31. **J. H. Churnside, “Aperture averaging of optical scintillations in the turbulent atmosphere,” Appl. Opt. **30**(15), 1982–1994 (1991). [CrossRef]

**32. **G. Rousset, F. Lacombe, P. Puget, N. N. Hubin, E. Gendron, T. Fusco, R. Arsenault, J. Charton, P. Feautrier, P. Gigan, P. Y. Kern, A.-M. Lagrange, P.-Y. Madec, D. Mouillet, D. Rabaud, P. Rabou, E. Stadler, and G. Zins, “NAOS, the first AO system of the VLT: on-sky performance,” Proc. SPIE **4839**, 140 (2003). [CrossRef]

**33. **K. Murphy and C. Dainty, “Comparison of optical vortex detection methods for use with a Shack-Hartmann wavefront sensor,” Opt. Express **20**(5), 4988–5002 (2012). [CrossRef]

**34. **K. Murphy, D. Burke, N. Devaney, and C. Dainty, “Experimental detection of optical vortices with a Shack-Hartmann wavefront sensor,” Opt. Express **18**(15), 15448–15460 (2010). [CrossRef]

**35. **J. D. Barchers, D. L. Fried, and D. J. Link, “Evaluation of the performance of Hartmann sensors in strong scintillation,” Appl. Opt. **41**(6), 1012–1021 (2002). [CrossRef]

**36. **J. Notaras and C. Paterson, “Demonstration of closed-loop adaptive optics with a point-diffraction interferometer in strong scintillation with optical vortices,” Opt. Express **15**(21), 13745–13756 (2007). [CrossRef]

**37. **J. D. Barchers, D. L. Fried, and D. J. Link, “Evaluation of the performance of a shearing interferometer in strong scintillation in the absence of additive measurement noise,” Appl. Opt. **41**(18), 3674–3684 (2002). [CrossRef]

**38. **D. E. Thornton, M. F. Spencer, and G. P. Perram, “Deep-turbulence wavefront sensing using digital holography in the on-axis phase shifting recording geometry with comparisons to the self-referencing interferometer,” Appl. Opt. **58**(5), A179–A189 (2019). [CrossRef]