Temporal statistics of residual wavefront variance of an adaptive optics system

We present the probability density function (PDF) for the residual wavefront variance of an adaptive optics system that includes control error, the fitting error of a deformable mirror, and Hartmann sensor detecting noise. The PDF is directly connected to adaptive optics system parameters and the spatiotemporal strength parameters of atmospheric turbulence, and it can be described as a generalized Chi square distribution. Our results provide a more precise theory for adaptive optics systems compared to the current theory based on the ensemble average. Thus, this study can contribute to the development of high-resolution and high-stability adaptive optics systems for astronomy and optical communications in the atmosphere.


Introduction
Adaptive optics (AO) systems are widely applied in astronomical imaging and atmospheric optical communications because of their ability to correct wavefront aberrations [1][2][3][4][5][6]. The current theory for AO systems is based on the ensemble average, and it provides the relationship between system parameters and average residual wavefront variance. Tyson and Hardy [7,8] presented most of the ensemble equations for AO systems in their books. The effectiveness of ensemble equations relies on a single fact: if we design an AO system with a mean wavefront variance of 1 rad 2 , 50% of images probably have a quality better than 1 rad. Moreover, there is a certain probability (∼10 −4 ) of obtaining excellent images with a Strehl ratio of 0.9. However, the exact probability for the distribution of wavefront variance or image quality can never be determined using only current ensemble average results.
A Strehl ratio of over 0.9 is required for the AO systems used in large telescopes for extrasolar planet observations [9,10]. However, building an AO system with an average Strehl ratio of 0.9 that uses visible wavelength, as expected in next-generation AO (NGAO) systems, is a technical challenge [11]. Lucky imaging may be required for AO; this would require research on how 'lucky' an AO system is. The current ensemble equations of AO cannot provide the temporal variation information of wavefront spatial variance, which makes the stability and accuracy required by NGAO difficult to analyze [11].
The stability of an AO system becomes extremely important when it is applied in optical communications through the atmosphere. Optical communications systems require high availability or stability. A given data package cannot be discarded as in astronomical imaging. If the AO characteristics for imaging and communications are compared, the 1000th-qualified (for a given threshold) AO wavefront compensation for imaging may be acceptable for a quasi-static astronomical object. In contrast, such a performance would mean an unacceptable bit error rate of 10 -3 for an optical communications link. Therefore, an AO system for imaging can have low availability, but an AO system for optical communications must have extremely high stability.
Generally, the large difference in information bandwidths for imaging and communications requires different analysis methods for AO systems in different scenarios. The frequency of astronomical imaging has a magnitude of 1 Hz. The closedloop bandwidth of an AO system is hundreds of hertz, and the bit rate for optical communications can be in gigahertz. Imaging quality is approximately the average of hundreds of frames with AO compensation. Thus, the ensemble method is applicable to AO design. On the contrary, millions of signal bits are contained in one sample interval of AO compensation. Hence, communications performance can only be precisely calculated with the probability distribution of the AO compensation metric, which has a nonlinear correlation with the bit error rate [12].
In 2008, Gladyz et al investigated the temporal variability and statistics for the instantaneous Strehl ratio of AO compensation and presented a three-parameter gamma distribution model for AO wavefront variance [13], which they used to discuss the distribution of the Strehl ratio. Canales and Cagigal investigated the speckle statistics of partially compensated AO system [14,15]. Their works gave the patterns of the AO correction performance indexes such as phase variance, Strehl ratio and received speckle on image plane. However, the relationship between the distribution model and AO system parameters in the form of equations has not been reported yet. In 2016, we presented the probability density function (PDF) of an AO system as a function of the closedloop bandwidth and atmospheric Greenwood frequency [16]. In this study, we integrate control error, fitting error, and Hartmann sensor detecting error to produce a comprehensive PDF equation for an AO system with determined parameters.
The principle of this work is based on a linear system description of an AO system, where the Zernike coefficients of a reconstructed wavefront are a linear combination of subaperture centroid vectors from different error sources. We can analyze the probability distribution of wavefront error from the randomness of the centroid errors corresponding to residual wavefront error.

Hartmann sensor detecting error
The typical structure of an AO system is shown in figure 1. A wavefront aberration is sensed by a Hartmann sensor. A wavefront controller controls a tilt mirror (TM) and a deformable mirror (DM) to correct atmospheric tilt and highorder aberration separately in real time. The four main components constitute the independent TM and DM closed-loops, and the residual wavefront error of the AO system is mainly determined by the error factors in these four components.
In general, only two main types of noise are analyzed in the Hartmann sensor detecting error: photon shot noise and readout noise [17]. We can neglect high-order aberrations and assume that only tilt exists on the subapertures of the Hartmann sensor when the AO is in a closed loop. Thus, all spots on the subapertures have an ideal Airy profile with different center shifts. As shown in figure 2, the coordinate on a subaperture is OXY after the Hartmann sensor is calibrated. The other coordinate, oxy, is located at the center of the exact Airy spot. Hence, r 0 denotes the centroid shift on the subaperture. Here, we only calculate the x value of centroid error. The centroid error on the y axis is statically identical to that on the  x axis. Therefore, the centroid error due to detecting error can be calculated from the sums on all X coordinates as where K denotes the received photon count at coordinate X, K 0 denotes the corresponding photon count without noise, and K − K 0 is the photon fluctuation caused by noise. If we define the spot centroid without noise as X 0 , then where the zero point of the x coordinate is the centroid of the spot, then equation (1) can be rewritten as Clearly, the centroid error due to the Hartmann sensor is determined by the signal-to-noise ratio and spatial distribution of noise.
Here, we want to focus the effect of noise, therefore, only tilts on Hartmann subaperture are considered by assuming the Hartmann sensor has adequate spatial sampling frequency. Otherwise, the sampling error of Hartmann subapertures should be included.

Photon shot noise
For photon shot noise, the received photon count, K, is randomly distributed as a Poisson distribution with expectation The expectation of the photon count on the subaperture is centrosymmetric because it is proportional to the intensity profile of the spot. Therefore, we can calculate the centroid error caused by shot noise in polar coordinates, as shown in figure 3. Then, we can express centroid error as where I is intensity, hv is photon energy, F is the frame frequency of the Hartmann camera, I r ( ) is the intensity function of the Airy spot, and where I 0 is the peak intensity, λ is the wavelength, f and d l are the focal length and diameter, respectively, of the microlens of the Hartmann sensor, and J 1 (·) is the first-order Bessel function of the first kind. For each subaperture, the denominator in equation (4) is the total photon count in one frame image. It follows the Poisson distribution, and its average is To calculate the value of the numerator in equation (4), we consider the random variable series in the inner part of the sum q -K K cos , ) are independent of each other, but they have different PDF functions and a corresponding expectation,  In equation (9), q | is always satisfied when ring width > dr 0. Thus, the limit of equation (9) is zero when positive integer p approaches infinity. Therefore, according to the Lyapunov condition of the central limit theorem, the sum of the inner term of the numerator in equation (4) where R is the radius of the area on the Hartmann sensor corresponding to one subaperture. Considering that equation (11)   For two independent random variables A and B, the variance [18]. In equation (4), the numerator obeys Gaussian distribution with zero mean value; the denominator obeys Poisson distribution, its variance and expectation å å . Therefore, if the inequality of the total photon count among all subapertures due to scintillation can be neglected, every element in the centroid error vector, z ph , caused by photon shot noise has the same variance. Its value can be calculated by connecting equations (4), (5), (7), and ) for the total photon count of the subaperture spot on Hartmann sensor, the variation of å K is about several å K , 0 which means å K could be approximated as a constant (i.e. å K 0 ) for a practical AO system in which the subaperture spots have over hundreds of photons. Therefore, the centroid error Dx in equation (4) approximately obeys Gaussian distribution.
Note that the equivalent Gaussian width of the Airy spot is s l = f d 0.45 .
A l Hence, equation (13) provides mostly the same results as Cao and Yu [19]. The centroid error variance of photon shot noise is inversely proportional to the total photon count, and it is affected by the parameters for the microlens of the Hartmann sensor according to equation (13).

Readout noise
The readout noise on Hartmann sensor pixels follows a Gaussian distribution. Thus, the centroid error from readout noise can be calculated according to equation (3): where P ij is the ideal photon count on pixel, n ij is the noise equivalent photon count on pixel ij, and x ij is the coordinate of pixel ij. The numerator in equation (14) is a linear combination of Gaussian random variables. The sum also follows a Gaussian distribution. Moreover, the centroid error caused by readout noise can still be considered as a Gaussian random variable, similar to the analysis of photon shot noise. Note that the zero point of the coordinate is located at the center of the spot. Thus, the mean value of the numerator in equation (14) is zero, and its variance is [19] where s rd 2 denotes the readout noise variance at each pixel and L and H are the numbers of pixels in the X and Y directions, respectively. The centroid variance caused by readout noise is inversely proportional to the square of spot power and proportional to readout noise variance.

Fitting error
An AO system is typically analyzed as a linear system. The temporal characteristics of the fitting error of a DM can be modeled in the same manner. This implies that (1) the fitting error of a DM is the linear sum of each Zernike mode when mode coupling is neglected, (2) the spatial shape is similar when single modes with different magnitudes are fitted, and (3) the magnitude of fitting error is proportional to the magnitude of the fitted wavefront.
To examine these three assumptions, we modeled a 61cell DM to fit the Zernike mode, Z x y , . m ( ) The fitting error for where H x y , k ( ) indicates actuator influence functions and b k denotes the related coefficients of actuators. The actuator arrangement used in the simulation is shown in figure 4(a). The white points denote the centers of the actuators. We used this DM to fit the fifth-order and 20th-order Zernike modes, Here, N is the number of subapertures of the Hartmann sensor and f is the focal length of the microlens. á ñ n · represents the spatial average on subaperture n. Here, we assume that the spatial frequency of the Hartmann sensor is unlimited; hence, aliasing error is neglected.
According to Noll's analysis, the Zernike coefficients of Kolmogorov turbulence a m have a Gaussian distribution with an average of zero [21]. Therefore, the total centroid error vector is the linear combination of all Zernike modes: where G(·) is the gamma function, D is the telescope diameter, and r 0 is the atmosphere coherent length (or Fried constant). We assume that the fluctuations of the Zernike modes are independent. Equations (18) and (19) indicate that the centroid error vector caused by the DM follows a Gaussian distribution. The variance vector corresponding to , , .

PDF for the residual wavefront variance
The Zernike coefficient vector, A t , due to insufficient control bandwidth has a Gaussian distribution [16], and the variance of each element is where f TG is the tilt Greenwood frequency of atmospheric turbulence, -f t 3dB is the tilt-control system bandwidth (3 dB cutoff frequency), f G is the Greenwood frequency of highorder Zernike aberrations, -f d 3dB is the bandwidth of the DM control loop, and M is the Zernike mode number in the wavefront reconstruction.
We can obtain the Zernike coefficient vector by adding control error, detecting error, and fitting error: where R is the constant reconstruction matrix of the AO system. All elements in the vectors on the right side of equation (22) follow a Gaussian distribution. Consequently, as the sum of a linear combination of Gaussian random variables, each element in A follows a Gaussian distribution with zero mean. Their variances, s ,   In equation (23), the atmospheric fluctuations of Zernike mode amplitudes are assumed to be approximately independent. Our simulation shows that the effect of the covariance of Zernike coefficients on the PDF of the residual wavefront is negligible. However, a strict expression can be obtained using Karhunen-Loève modes, in which equations (17) and (19) are replaced by the results of the Karhunen-Loève expansion [22,23]. If we use the standard Zernike polynomials with unit variance, the total residual wavefront (spatial) variance of an AO system is equal to the sum of the squares of the elements of A: According to probability theory, the temporal variation in the residual wavefront (spatial) variance of an AO system follows the generalized c 2 distribution. Its PDF is extremely complicated, and it could be approximated by a single-gamma distribution with the same mean and variance [24].
We can also rewrite equation (24) in the following form:  ( ) After the PDF of the wavefront variance of an AO system is obtained, the index of availability (or stability) could be calculated by the integral of the PDF of wavefront variance from zero to a given upper boundary, as discussed in our previous work [16].
We examined the PDF equations using the residual wavefront data of an AO system for a 1.8 m diameter telescope operated by the Adaptive Optics Laboratory of the Chinese Academy of Science [25,26]. This AO system consists of two closed loops to independently control its TMs and DMs with 127 actuators using a proportional integral controller. The Hartmann sensor runs at 0.55 μm with an image acquisition frequency of 2000 Hz. The data were measured for two stars with magnitudes of 2.9 and 1.6, and the corresponding atmospheric coherent lengths were 7.9 cm and 9.3 cm at 0.55 μm, respectively. Each measurement consisted of 10 000 frames of wavefront data acquired in 5 s.
The wavefront data were reconstructed from the first 35 Zernike polynomials. We found that the histogram of every Zernike coefficient strongly agreed with the Gaussian function with different variances. The PDFs of residual wavefront variance were calculated using the wavefront data and theoretical equations. The results are shown in figure 5.
As shown in figure 5, the theoretical analysis of the PDF of residual wavefront variance agrees well with the closedloop AO system data. Comparing the PDF curves between Mag=2.9 and Mag=1.6, their residual high-order aberrations are almost the same. The difference mainly arises from different tilt correction performances, which show a large tail on the right of the PDF curve. The large tail would considerably degrade long exposure image quality and receiving efficiency in optical communication. Therefore, the primary object of AO optimization is to eliminate the tail of the PDF curve of wavefront variance. In other words, we should improve the tilt correction ability of an AO system first by increasing the control bandwidth of the TM, as discussed in [14]. Then, we would expect the PDF curve to shift to the centrosymmetric shape along with decrease in residual tilt. This issue could be analyzed from the point of view of the skewness of the PDF curve.
The skewness of the PDF curve of AO residual wavefront variance can also be theoretically estimated from its PDF equations. The skewness of f A m If residual tilt can be neglected for a fine tilt correction loop, residual higher-order Zernike terms dominate the residual wavefront, and s m 2 is almost the same for higher-order terms, then the skewness of the PDF of residual wavefront variance s wf 2 is approximated by The PDF curve of residual wavefront variance approaches centrosymmetry for large M. However, because of the negative exponential calculation from wavefront variance to the Strehl ratio, the PDF of the Strehl ratio may have negative skewness for a small average of s .

Discussion and conclusions
In this work, we obtained the Gaussian characteristics of Zernike coefficients of AO residual wavefront from the randomness of error sources. Furthermore, we obtained the connection between the distribution model and AO structure parameters, which made it possible to predict the distribution model at the beginning of system design.
The AO theory is mainly constructed according to the ensemble average, which is used to estimate average AO system performance and for system design in the past. Through the linear system analysis for an AO system, we developed a group of equations that connected the PDF of AO residual wavefront variance with AO system parameters and turbulence parameters. These equations could directly guide AO design and optimization with higher precision. Moreover, the PDF was a complete mathematical description for a random process. Therefore, the results in this paper provide a helpful framework for investigating the principle of AO systems from the viewpoint of probability for all AO applications.
In future, more error sources, such as aliasing error, discrete sampling error, and isoplanatic error, must be considered for furthering investigating the probability characteristics of the random process of AO compensation. Certain other topics are also important, such as the errors due to nonintegral controllers and non-Hartmann sensing in AO systems. Moreover, the boundary conditions for each type of error should be clarified for different AO application scenarios and specific configuration of AO system. This will help to balance the error budget and optimize system performance, and finally to construct the complete probability theory of the AO system.