Synthetic Modeling of Astronomical Closed Loop Adaptive Optics

We present an analytical model of a single natural guide star astronomical adaptive optics system, in closed loop mode. The model is used to simulate the long exposure system point spread function, using the spatial frequency (or Fourier) approach, and complement an initial open loop model. Applications range from system design, science case analysis and AO data reduction. All the classical phase errors have been included: deformable mirror fitting error, wavefront sensor spatial aliasing, wavefront sensor noise, and the correlated anisoplanatic and servo-lag error. The model includes the deformable mirror spatial transfer function, and the actuator array geometry can be different from the wavefront sensor lenslet array geometry. We also include the dispersion between the sensing and the correction wavelengths. Illustrative examples are given at the end of the paper.


INTRODUCTION
Astronomical adaptive optics (AO) are complex opto-electro-mechanical systems, designed to correct random aberrations generated by optical turbulence in earth's atmosphere, and improve the performance of astronomical telescopes -see Roddier et al. [1] and figure 1 for an illustration. Modeling the performance of such a system -for science cases analysis, AO data reduction and system design -requires sophisticated simulations tools. The most intuitive approach to construct an AO simulation tool is to break-down the system in its fundamental components, build a physical model for each of these components, and link all of them following a block-diagram architecture. The turbulent phase is then propagated through the model, mimicking a real system, up to the system's output, the focal plane, where the system's point spread function (PSF) is measured. The capacity of such end-to-end models to produce accurate performance predictions is in principle limited only by our capability to accurately code the behavior of each sub components and the system's input disturbance (optical turbulence, noise, other) -see for instance Carbillet et al. [2] and Le Louarn et al. [3].
End-to-end models have a limitation, though. Indeed, as the input of the system (optical turbulence) is of stochastic nature, the model needs to be run on a large number of instantaneous turbulent optical waves, typically more than a thousand to converge to a Figure 1: Top: basic elements of an AO closed loop system. The incoming turbulent wavefront is transmitted by the telescope optics to the entrance of the AO system. The deformable mirror surface shape (DM) is set by the control computer to compensate for the wavefront error. A fraction of the corrected light beam is transmitted by the beamsplitter (BS) to the wavefront sensor (WFS) where the residual wavefront is measured. The computer reconstruct the residual wavefront from the WFS measurement and computes/updates the DM surface to keep the wavefront residual as small as possible. The other part of the corrected beam is sent toward the science instrument. Bottom: equivalent block diagram of the AO closed loop system, indicating the main elements of the loop. R indicates the wavefront reconstruction operation, C is the control algorithm (an integrator in this paper).
pseudo-long exposure PSF from which performance metrics can be estimated. Such a procedure takes a lot of time: depending on order of the AO system, which defines the size of the command matrix, and the level of sophistication of the end-to-end model, getting a PSF without too many residual speckles can take several hours or even days. As a consequence, performance analysis using end-to-end tools is generally limited to a few well selected representative cases, and extensive studies of the system parameter space is rarely undertaken. Also, because all the error sources are naturally merged in an endto-end model, it is not easy to disentangle the impact of the individual sources of error on the overall performance, unless one has a deep understanding of how the end-to-end model was built. End-to-end models are therefore not a good choice for AO engineers for the broad analysis of a system performance, nor for astronomers interested in exploring the capability of a given AO system (existent or planned) for their science programme.
To suppress the limitation of end-to-end models, Rigaut et al. [4] proposed a totally different method, that we refer to here as the synthetic approach. Instead of letting the optical wave propagate through the system components and observe what comes at the output, we built a model for the system's output itself. The synthetic approach is based on an understanding of the system's behavior, and its accuracy is only limited by this understanding. A priori knowledge is critical for the synthetic approach, while this is not needed with the end-to-end approach. We can also said that the synthetic approach models the behavior of the system, while the end-to-end approach models the structure of the system.
The starting/central point in the synthetic approach is the construction of an analytical model for the long exposure (or average) AO-corrected phase spatial power spectrum density (s-PSD). AO correction is actually seen as a spatial filter applied on the turbulent phase s-PSD. This approach is therefore also called the Fourier method or spatial frequency method in the AO literature. From this residual phase s-PSD it is shown in this paper how the long exposure PSF can be computed, in a few steps. Getting the long exposure PSF is therefore very fast and exploring in detail the AO system parameter space becomes possible. Also, wavefront error budgets are easily build, because a s-PSD is attributed to each error source, from which we can get the wavefront error variance, by numerical integration in the spatial frequency domain. Finally, by nature, there are no residual speckles in a synthetic PSF: performance metrics (Strehl ratio, PSF width, integrated energy) are therefore not "noisy", and science case analysis are greatly facilitated.
This being said, the synthetic approach has its limits: non linear effects cannot be modeled, neither transient temporal behaviors. Also, the fundamental assumption is that the corrected phase is stationary within the pupil (which is not true near the edge of the pupil and for low order aberrations) and this can produce pessimistic performance predictions. For these reasons, end-to-end and synthetic models are to be seen as complementary rather than competitive methods: broad exploration of the system's performance is the domain of synthetic models, while detailed analysis of specific aspects of the system definitely requires end-to-end models. Actually, inclusion of our synthetic model within an end-to-end Monte-Carlo code has already been tried successfully by Carbillet et al. [5].
In an earlier work [6], we complemented Rigaut et al.'s work, explaining the foundations of the method in great details, and including PSF modeling in two dimensions. This initial model was open loop, where the wavefront sensor (WFS) measures the wavefront aberration before the deformable mirror (DM) correction, while the vast majority of current systems operate in a closed loop mode, i.e. the WFS measures the residual wavefront error after compensation by the DM -see figure 1.
Closed loop systems control the wavefront error and are therefore relatively insensitive to external disturbances (like noise) and internal variations of system parameters. Open loop systems on the contrary are very sensitive to modeling errors: it is critical that the system's components behave the way they are supposed too, as the quality of the correction is not controlled. On the other hand, feedback of the wavefront error in a closed loop system can generate diverging instabilities, if the error is overcompensated, or if time delay in the loop is too large. Open loop systems do not have this stability issue. See Ogata [7] for an introduction on control systems.
In bright guide star conditions, the AO system can be run at a high loop rate: the servo lag error (due to the time lag between the measurement and the actual correction) is low, and the noise level is negligible. In this case there are no differences between open and closed loop performance. For a dim guide star, the WFS exposure time is increased to gather more photons and keep the signal/noise ratio of the wavefront error measurement at an acceptable level: servo-lag error increases, and because of differences in open and closed loop transfer functions, the system performance significantly changes between the two modes. Because of this different behavior, open loop models cannot be used to predict close loop performance in high noise regime. As predicting the limiting magnitude of a system is of great importance, in particular for science cases studies, we have developed further our initial model and included closed loop modeling.
We also took this opportunity to introduce a DM spatial frequency transfer function, allowing the analysis of different influence function structures and actuators grid architecture, and the dispersion error, i.e the error induced by the air's refractive index dispersion, which makes that the wavefront measured at a given wavelength is slightly different from the wavefront corrected at any other wavelength, generating a non negligible error for tight error budget AO systems. Also, a few conceptual errors that appeared in our initial open loop paper are corrected.
Our closed loop model is developed for a single natural guide star (NGS) Shack-Hartmann WFS based system (SH-WFS), and a least square error (LSE) wavefront reconstruction algorithm. This case covers the vast majority of current NGS-based AO systems design, and in any case our model can be easily adapted to other schemes. Examples of usage of the synthetic method are given at the end of this paper, and we show how this method can be used to optimally dimension a system.
We finally note that the synthetic approach has been used and developed by other authors as well, increasing the diversity of views and understandings of the limits and strengths of this method. We must mention first the work of Ellerbroek [8] who basically developed the same method using a more general albeit potentially less detailed approach; the work of Rigaut [9], Tokovinin [10] and Jolissaint et al. [11] for a ground-layer AO mode (but not in closed loop and without wavefront sensor noise model); and more recently Neichel et al. [12] who introduced multiple guide star tomographic reconstruction, a very useful extension of the method. What our model brings to these recent developments is essentially the closed loop mode, and some useful sophistications like the DM transfer function. To finish, it is fair to mention that R. Conan and Ch. Verinaud (private communication) both independently developed a synthetic closed loop model, yet unpublished, using another but equivalent approach than the one presented here, namely the equivalence between the spatial and temporal frequency through the Taylor hypothesis (see text).

FOUNDATIONS OF THE SPATIAL FREQUENCY METHOD: A SUMMARY
A detailed description of the foundations of the method is given in Jolissaint et al. [6]. A summary is given here for convenience.
The method is based on the relationship between the phase spatial frequency power spectrum and the phase structure function (SF) in one side, and the SF and the average (long exposure) AO system optical transfer function (OTF) on the other 2 . Let us review the different steps from the phase s-PSD to the long exposure PSF.
Analytical expressions for the s-PSD will be developed later. The starting point is the OTF of the whole system made of (1) the column of air above the telescope, seen as an optically transparent medium with a turbulent field of refractive index, (2) the telescope optics, possibly with static aberrations, (3) the AO system optics, (4) the science instrument optics, that can be merged with the telescope optics. Splitting the phase aberration into a static, constant part ϕ and a turbulent, zero mean part δϕ, ϕ(r, t) = ϕ(r) + δϕ(r, t) where r is the position vector in the pupil plane and t is the time, and remembering that the OTF is also given by the autocorrelation of the phasor exp (iϕ) in the pupil plane [13], we get, for the long exposure OTF (averaged over an infinite number of realization of the 2 remember that OTF and PSF are Fourier transforms of each other random AO corrected turbulent phase) where ν is the angular frequency vector in the focal plane, associated to the spatial shift ρ = λν in the pupil plane (λ is the optical wavelength), S p is the pupil area and P(r) is the pupil transmission (1 inside the pupil, 0 outside), and · indicates a time or ensemble average. Assuming a Gaussian statistics for the phase aberration, which is a very good assumption for the uncorrected as well as for the AO corrected phase, it is shown in Roddier [14] that the average can be moved into the exponential argument, and we get defines the phase structure function, a measure of the variance of the phase difference between two points separated by a vector ρ in the pupil. We see that the structure function depends not only on the separation vector ρ, but also on the location r where the phase difference is measured. Therefore, if we want to compute the long exposure OTF, we need a model in (r,ρ) of the structure function, which is not necessarily difficult to obtain, analytically, but what is more annoying is that we need then to perform a numerical integration of Eq. (3) over r for each angular frequency ν = ρ/λ. This can be a very time consuming effort, and goes against the very objective of the synthetic approach.
Now, it is demonstrated in [14] that the optical turbulence phasebefore AO correction -is stationary over the pupil, i.e. its statistical properties do not depend on the location r inside the pupil. The corrected phase, on the other hand, is not stationary, and its residual variance increases from the center to the edge of the pupil. This being said, this non stationarity affects mostly the first orders -piston, tip-tilt, defocus ... -and the corrected phase can be considered to be stationary for a moderate to high order AO system (i.e. moderate to high Strehl). If the phase is stationary, which we will assume from now on, the phase structure function can be written D ϕ (r, ρ) = D ϕ (ρ) and the structure function exponential can be extracted from the integral in Eq. (3), so we get The system's OTF can therefore be seen as the telescope OTF filtered by an AO OTF filter, exp[−D ϕ (ρ)/2]. Separating the OTF this way is equivalent to splitting the overall optical system into two independent systems: the optical turbulence plus AO system on one side, and the telescope plus instrument optics on the other. The first system is therefore not related to the pupil optics in any way, and its description does not include the pupil boundaries anymore. For this very reason, synthetic modeling is sometime referred to as infinite aperture modeling.
We discuss now the relationship between the stationary SF and the phase s-PSD. Thanks to the stationary assumption, the AO system can be considered as an optical system applying a correction on a turbulent phase, regardless of any beam boundaries, all over an hypothetical plane perpendicular to the telescope optical axis. Everything looks as if the phase was pre-corrected by the AO system before being intercepted by the telescope beam. Now, it is shown in [15] that the stationary structure function D ϕ is related to the spatial correlation of the phase, B ϕ , which is itself equal to the Fourier transform of the phase s-PSD (written Ξ ϕ here), and we get The integral of the phase s-PSD gives the phase variance, and the cosine term is there instead of the usual complex exponential we have in a Fourier transform, because as the phase s-PSD is even, the sine component of the FT would be naturally null. With this last equation, we have completed the link between the phase s-PSD and the long exposure PSF.
To summarize, the procedure to get the long exposure PSF from the phase s-PSD is the following:

The fundamental equation of adaptive optics
The starting point for the development of the phase s-PSD is the so-called fundamental equation of adaptive optics, which states that at any instant t, the residual wavefront error w e is the difference between the incoming atmospheric turbulent wavefront w a and the mirror command 3 w c -see figure 1, where r is the location in the pupil plane, θ is the angular separation between the science object and the guide star (assumed on-axis without loss of generality), λ s is the science observation wavelength, and λ m is the wavefront sensing wavelength. Note that in our initial paper, we used the phase instead of the wavefront in the fundamental equation, but we believe now that using the wavefront formulation is more appropriate because it is actually the wavefront which is corrected in an AO system. We will therefore develop equations for the residual wavefront s-PSD, which will be transformed into the phase s-PSD by multiplication with the usual factor (2π/λ s ) 2 . Polychromatic PSF will be modeled by computing and averaging the phase s-PSD over the chosen optical bandwidth.
Including air refractivity The air's refractive index depends (slowly) on the wavelength, and measuring the wavefront at a different wavelength than the science observation channel introduce a small but noticeable error for systems with a tight wavefront error budget. Formally, the air's refractivity (N=n-1) is given by the sum of the refractivity of the air's constituents (nitrogen, oxygen, water, carbon dioxide etc.) Practically, though, it is shown in [16] that N can be written as the sum of a continuum and anomalous terms. The anomalous terms are associated with the excitation/absorption lines of water vapor and carbon dioxide (others constituents have a negligible impact), and because the atmosphere is naturally opaque at theses wavelengths, the anomalous terms are of no interest to us. The origin of the continuum term is actually not different from the anomalous terms: it is a sum of the wings of the strong nitrogen/oxygen/ozone excitation lines in the ultraviolet which extend to the visible and infrared wavelengths. There are very good theoretical/empirical models for this continuum [17,18] in the visible and infrared up to 10 µm, that are function of the air temperature, pressure, humidity and carbon dioxide content. We will not dig here into these models. What is of interest for us is that these models can be written as the product of a chromatic term which depends only on the wavelength, and a non-chromatic term which depends on the other variables (temperature, pressure etc.) Therefore, within the transmission windows of the atmosphere (i.e. inside the photometric bands), the wavefront error is simply proportional to the air's refractivity, and we can rewrite the fundamental equation of AO as where we define ν(λ m , λ s ) ≡ N(λ m )/N(λ s ) as the dispersion factor. The later formulation allows us to develop the model of the mirror command w c for a single wavelength -here the science wavelength -and correct for the fact that the wavefront is actually measured at another wavelength. The effect of dispersion in discussed with more details in Jolissaint and Kendrew [19].

The residual wavefront error in closed loop mode
Continuous process assumption AO control is a discretized process: the DM shape is updated periodically at a loop period ∆t, with a time delay t lag following the end of the wavefront sensor exposure. In-between these updates, the DM shape is kept constant (at least classical systems are working this way). AO control is therefore an integral and hold process, in a sense that the updated DM shape is equal to the previous one plus a weighted estimate of the wavefront residual. Now, as our approach is stationary in nature, no specific instant can be privileged: in other words, the equations we are writing need to be applicable at any instant, therefore the discrete integral control equation where c k represents the DM command at instant t k and g × e k the error signal weighted by the loop gain g, needs to be replaced by the continuous integral . This approximation is equivalent to assume that the closed loop control is applied continuously, as if at any instant t, the DM is updated with a command computed from the WFS measurement an instant t − t lag earlier.
The deformable mirror spatial transfer function The command applied to the DM is made of the projection of the WFS measurement onto the DM space, whose basis is the N-dimensional set of the DM influence functions I i=1...N . The equivalence of this projection in our stationary approach is the convolution of the reconstructed wavefront with the DM spatial response, or, in the spatial frequency domain, the multiplication of the wavefront Fourier transform with the DM spatial transfer function. The DM spatial response is defined by the projection of the Dirac impulse onto the influence function basis The coefficients p i are computed from the minimization of the quadratic distance between the Dirac impulse and its projection, and we find where S is the covariance matrix of the DM influence functions.
Note that as the DM response, by definition, is not supposed to vary across the DM surface, while actually the projection of the Dirac impulse does depend on its location u within the actuators grid, we will simply consider the average DM response over all possible locations u. We have found that this ad-hoc procedure generates DM transfer functions that better represent transfer functions measured on real systems. Practically, as the DM actuator array is periodic, we find that b i as the average of the influence function I over the square space centered on the optical axis and of width equal to the inter-actuator pitch. For another actuator grid geometry, the space over which the DM response is averaged would be different.
The DM spatial transfer function is given by the Fourier transform of the DM spatial response, and as the i-th influence function at a position r i can be written as the central influence function I 0 shifted by −r i , we find where I 0 (f) is the Fourier transform of the central influence function. There are numerous models for the DM influence function, all depending on the type of DM technology -see for instance [20]. We have developed our own empirical models for a Xinetics, Inc. 177 actuators DM model and a Boston MEMS 144 actuators model, with a modeling error of less than about 0.1 % in amplitude, and used these to compute the DM transfer function. A cut through these two DMs response and transfer functions is shown in Figures 2 and is compared with Gaussian and pyramid influence function models.
It is important to note that any other DM basis of independent functions can be used to build the DM transfer function: the choice we made here of using the influence function basis was dictated by the fact that we had at our disposal accurate empirical models of several influence functions types. Now, as pointed out by a reviewer of this paper, influence functions are not always a good choice to model the DM surface, and we certainly agree that it is particularly true for low order modes (for instance tilt is badly represented by an addition of influence functions). If correction of the low order aberrations is of particular concern for the modeling of a particular system, it might therefore be worth to use a basis particularly adapted to the representation of these modes (for instance, as proposed by the reviewer, a DM surface generated by a cubic spline interpolation).
Now, let us develop the DM command w c , including the continuous assumption and the DM response. The DM command w c at any instant t is the update of the previous DM command (an instant ∆t earlier) with the residual wavefront error, weighted by the loop gain g loop where * is the convolution product, R is the operator associated with the wavefront reconstruction from the WFS measurement m, and n is the WFS measurement noise. Note that this equation is not specific to any type of WFS nor reconstruction operator.
The WFS output m is defined by the measure (direct, gradient or laplacian) of the wavefront residual error w e averaged over the WFS integration time ∆t, delayed in time by the lag t l due to the WFS readout time and the command computation. Using the notation q to indicate a time average with a time lag, we write the WFS measurement as the application of a wavefront measurement operator M on the instantaneous residual error w e in the direction of the NGS, which is assumed on-axis (i.e. θ=0), Separating the deformable mirror space and the wavefront analysis space For further developments, we need now to split the atmospheric wavefront into the component which is corrected by the DM, and the component which is simply reflected off the DM surface, Beside, the WFS samples the wavefront with a spatial interval Λ WFS equal to the lenslet separation distance. This naturally split the spatial frequency domain into a low spatial frequency domain, i.e. the frequencies that can be seen by the WFS, and therefore corrected, below the WFS spatial cutoff frequency and a high spatial frequency domain, above f WFS . For a Shack-Hartmann WFS, the lenslet array has a square geometry, therefore the low spatial frequency domain is defined by the inequalities |f x | ≤ f WFS and |f y | ≤ f WFS , a square, and the high spatial frequency domain by the complement of this square. One finally gets, with Eq. (16), The later formulation allows the separation of the DM actuators and WFS lenslet array grid architectures, which can be now studied independently.
From this point, our model development is done in the spatial frequency domain only. Our final objective is indeed to write an equivalent spatial frequency power spectrum filter for the AO correction. In the spatial frequency domain, the fundamental equation of AO becomes, whereq indicates the Fourier transform of q, the DM command becomes and the WFS measurement Inserting Eq. (21) into Eq. (20), with Eq. (18), and as the reconstruction and measurement operators become spatial filters in the frequency domain, it comes (we drop momentarily the common variables to shorten the notation) In our approach, the product Γ DM (1 − Γ DM ) would describe the projection onto the orthogonal of the DM space, followed by the projection onto the DM space. This product naturally has to be replaced by the null operator. Beside, the wavefront reconstruction action is to revert the WFS measurement, therefore, for a wavefront q LF strictly limited to the WFS low frequency space, the cumulated operation of wavefront measurement and wavefront reconstruction is the identity operator, i.e. R{M{q LF }} = q LF . With the later remarks, and noting that the DM command, by nature, belongs to the low frequency WFS space, Eq. (22) simplifies to A note on the loop gain Nothing prevent us from setting the loop gain g loop here as a free variable in f, too. This allows optimization of the loop gain frequency-by-frequency, which is equivalent -in our stationary approach -to modal gain optimization. We will therefore keep the notation g loop (f) even if loop gain optimization is not discussed further in this paper.
Independence of the turbulent layers and Taylor hypothesis In the atmosphere, optical turbulence is distributed in thin independent layers, each being characterized by (1) the so-called refractive index structure constant C 2 N,i , a measure of the variance of the refractive index spatial fluctuation within the layer, (2) the apparent velocity of the turbulent layer -see for instance [21]. As seen from the pupil of the telescope, the time scale over which the wavefront associated to each layer evolves significantly is generally longer than the time it takes for the wind to push the layer across the telescope beam. Therefore, in first approximation, everything looks as if the optical turbulence profile was made of a certain number of frozen wavefront screens translating across the telescope beam with the layers wind speed and directions. This assumption of frozen optical turbulence layers is referred to as the Taylor hypothesis in the literature 4 , and has the nice consequence that it is possible to transpose, within the layer, a shift in time ∆t into a shift in space ∆r = v∆t, with v the layer's wind velocity.
Assuming independence of the turbulent layers, the correction of the total wavefront summed over the N l turbulent layers 5 is equivalent to the correction of each wavefront from each layer taken individually, as in any case, cross terms between the layers will vanish on average in the computation of the long exposure residual phase s-PSD. Let us 4 departure from the Taylor hypothesis is discussed in [22] 5 about 10 layers are generally needed to model a turbulent profile therefore compute the wavefront error spectrum w e,l associated with the layer l.
As shown in [6], in the Fourier domain, the time average of a wavefront q over a time interval ∆t, followed by a time lag t lag becomes, for the layer l with a wind velocity v l , where ∆t/2 + t lag , equaling the time interval between the middle of the WFS exposure time and the application of the new DM command, represents the overall time lag. With the later, the DM command spectrum Eq. (23) becomes, for the layer l, where we have replaced the time shift ∆t within w c,l (t−∆t) by its equivalent phase change exp (2πi ∆t f · v l ) in the Fourier domain, using the Taylor hypothesis.
The noise term, being added to the wavefront slope measurement, is not linked in any way to the turbulent layers, therefore it does not make any sense to write a noise term for each layer. The noise term is independent from the other error terms and needs to be treated separately, so we do not include any noise term in the last equation. The overall WFE will simply given by the sum of the servo-lag contribution and the noise contribution.
Regrouping the terms in w c,l , we end up with an expression for the DM command Inserting Eq. (26) into Eq. (19), we get, again for the layer l, As an angular shift θ is seen, at the layer altitude h l , as a spatial shift ∆r = h l θ, using the shift theorem of the Fourier transform we get w a,LF,l (f, θ, λ s , t) = exp (2πi h l f · θ) w a,LF,l (f, 0, λ s , t) We can now rewrite Eq. (27), w e,l (f, θ, λ s , t) = w a,HF,l (f, θ, λ s , t) high order WFS "fitting" error + F AS,l (f) w a,LF,l (f, 0, λ s , t) aniso-servo error + F AL,l (f) w a,HF,l (f, 0, λ s , t) WFS aliasing error (29) and identify the four fundamental terms of the residual wavefront error: 1. the high order WFS error, usually named the "DM fitting error" in the AO literature -which is actually for us the part of the atmospheric wavefront which is not seen by the WFS, therefore cannot be corrected by the DM, so we think that calling this error the high order WFS error is more appropriate, 2. the angular anisoplanatic AND loop servo-lag error, identified here as the "anisoservo" error, as anisoplanatism and servo-lag error are correlated (with the Taylor hypothesis, a wavefront shift in time can be compensated by a negative wavefront shift is space), 3. the WFS aliasing error: the high order wavefront error is seen by the WFS as a low spatial frequency error and reconstructed as such, therefore the AO system is compensating an error which actually is non-existant, 4. and finally the WFS noise term, discussed later.
F AS,l is defined as the aniso-servo spatial filter for the layer l, and F AL,l is the WFS aliasing spatial filter for the layer l It is interesting to examine the limits of the aniso-servo and WFS aliasing spatial filters when there is no angular separation between the science object and the NGS, i.e. θ=0, and when the WFS integration time and loop lag are set to zero, ∆t=t lag =0. We find which indicates that in the absence of aniso-servo error, the residual low frequency wavefront error is generated by (1) the refractive index dispersion -and we see as we would expect that this error is proportional to the dispersion factor, and (2) the aberrations seen by the WFS that the DM cannot correct, which are never null because even if the DM actuator pitch is equal to the WFS lenslet pitch, a perfect correction of the low spatial frequencies would require a DM with sinus cardinal influence function (Fourier transform of a sinus cardinal is a door function), which is only approximated by actual influence functions -see Figure 2. The same is true for the WFS aliasing error. This behavior corresponds well to what we would have expected.
A note on the correlation between the error terms Computing the s-PSD associated with the four fundamental wavefront error terms above essentially consists in computing the modulus square of Eq. (29), averaged over the time t. Therefore, cross products appear between the four error terms. Now, the noise term is naturally not correlated with the other errors and can be treated separately. The high order WFS error is not seen by the system and transmitted to the output of the system, unaffected. We will assume in this paper, without discussing it further, that the cross products between the low and high order spatial frequencies are negligible relative to the main error terms.
Therefore, in what follows, the four terms above will be discussed independently from each others.
We will now make use of the expressions developed above for the wavefront error terms to develop the analytical expressions of the residual phase s-PSD of the four fundamental errors we have identified. As the residual phase s-PSD is given by spatial filtering of the optical turbulence phase s-PSD, we will start by recalling the expression of the later, as given in the literature.

The turbulent phase spatial power spectrum
The s-PSD of the turbulent phase is discussed in Roddier [14]. In the atmosphere, the extension of optical turbulence is necessarily limited by the individual layers thickness, and an optical turbulence outer scale L 0 was included in the Kolmogorov s-PSD to account for this spatial limitation. This modified Kolmogorov s-PSD is called the von Karmàn s-PSD in the literature (see for instance Winker [23], and Maire et al. [24] for a few other s-PSD models) and is given, at the wavelength λ, by Ξ ATM (f) = 0.0229 r 0 (λ) −5/3 (|f| 2 + 1/L 2 0 ) −11/6 (32) where r 0 is the Fried parameter, a measure of the strength of turbulence, defined as the telescope diameter whose focal plane angular frequency cutoff would be the same than the optical turbulence cutoff frequency (see Fried [25]). r 0 is generally given at 500 nm in the literature, and we will follow this convention, unless indicated differently. Typical values for r 0 at 500 nm extend from 5 cm (bad observation site, day-light conditions) to 25 cm (excellent site). The optical outer scale is generally in the range 20 to 40 m, surprisingly with very few variations between the different sites where this quantity has been measured.
The Fried parameter is associated to the vertical profile of the optical turbulence structure constant C 2 N (h), following which can be written, in the case of N l independent optical turbulence layers, as a sum where ∆h l is the layer's thickness, and r 0,l defines the layer's Fried parameter. Consequently we can also define a phase s-PSD for each layer, which naturally sums up to the overall phase s-PSD, Using this notation, the phase s-PSD associated with the low and high order turbulent wavefront errors -w a,LF,l and w a,HF,l -will be written now where µ LF and µ HF are low and high spatial frequency masks, defined, for a Shack-Hartmann WFS, by the square domain and 3.4 The high order WFS spatial power spectrum -or "fitting error" The high order WFS phase error s-PSD is simply given by the atmospheric turbulence phase s-PSD limited to the high spatial frequency domain and is written where Ξ ATM is given by Eq. (32).

The aniso-servo spatial power spectrum
The aniso-servo phase error s-PSD is given by the time average of the modulus square of the aniso-servo wavefront error, translated into a phase error, summed over the N l turbulent layers. From Eq. (29) and Eq. (37), we get where Ξ ATM,l is given in Eq. (35) and While this equation seems impressive, coding it into a computer program does not represent a particular challenge.

The WFS aliasing spatial power spectrum
The aliasing error is given, for the layer l, by the term -see Eq. (27) the product R M w a,HF,l was already developed in our initial paper, and is recalled here. A Shack-Hartmann WFS produces a measurement of the wavefront slope (gradient) in both x and y directions, with a spatial sampling given by the WFS lenslet array pitch Λ WFSwhich is also the lenslet width. In the spatial frequency domain, the slope measurement is given by the multiplication of the wavefront Fourier transform with the two components operator (one for each direction) where the product with the spatial frequency vector f stands for the derivative in the Fourier domain, the sinc function is for the wavefront average over the lenslet area, and III(Λ WFS f), the Dirac comb, represents the recurrence of the measured slope spectrum with a spacing 2f WFS in both x and y directions, and is responsible for the aliasing of the part of the spectrum above the WFS cutoff frequency f WFS inside the low spatial frequency domain -which always occurs, because the turbulent wavefront spectrum is not band limited at high spatial frequency.
The analytical expression for the reconstruction operator Fourier transform R is computed from the minimization of the quadratic distance between the slope measurement and the actual slope (least square error algorithm, LSE). The weakness of the LSE algorithm is that the WFS noise is reconstructed as a real signal, without penalty. Other algorithms have therefore been proposed that make use of a priori knowledge of the Kolmogorovstatistics based signal and noise statistics to minimize the contribution of the noise on the reconstructed signal -see for instance [12]. A discussion of the pros and cons of these different algorithms is beyond the scope of this paper, though, so we will stick with the LSE-based algorithm, as it is the most simple and most straightforward to implement.
As we have seen, the slope measurement operator in the Fourier domain is, basically, a multiplication with the spatial frequency vector f. The reconstruction operator is therefore the inverse operator, i.e. the inverse of the vector 6 f, but ignoring the Dirac comb, because the reconstruction does not extend beyond the WFS cutoff frequency f WFS , and we find (the factor Λ WFS disappears because it is actually part of the Dirac comb convolution product) We can now develop the term R M w a,HF,l from Eq. (44). Using the two equalities and we find where it is important to note that |m| + |n| > 0 because we do not want to include the low frequency part of the wavefront spectrum in the sum (as it is of course not aliased).
We can give now the expression for the WFS aliasing phase error s-PSD. Summed over the N l independent layers, it is given by where, for each layer, where we made the assumption that the correlation of the phase for frequencies separated by a 2f WFS interval is negligible. It is important to realize that the term T =

The WFS noise spatial power spectrum
From Eq. (23), it comes What is the noise n(r, t) made of, exactly ? it is a discrete quantity, in space and time, made of two components in x and y, n(r, t) = [n x (r, t), n y (r, t)] sampled on a grid of spacing Λ WFS . Its spatial spectrum is therefore necessarily limited to the domain |f x |, |f y | ≤ 1/(2Λ WFS ). As the noise over the lenslets is uncorrelated, all values are possible at any instant and location, therefore the noise s-PSD is necessarily white, i.e. it is constant in the domain |f x |, |f y | ≤ 1/(2Λ WFS ), and is the same for both x and y components. So, we can define such that |fx|,|fy|≤ 1/(2Λ WFS ) so, where σ 2 NEA, CL is the closed loop noise equivalent angle (NEA) variance, discussed later. The noise s-PSD is given by the average modulus square of w NS . With the reconstructor Fourier transform -given in Eq. (46), it comes therefore we get, from Eq. (53), Note that the loop gain does not appear anymore directly in the noise s-PSD. Indeed, the noise s-PSD has to be seen as a spatial filter, actually not different from its formulation in open loop, but where the NEA is now a closed loop NEA. In other words, it is the NEA which is affected by the closed loop noise transfer function, not the spatial properties of the noise s-PSD. Note also that there is no analogy between the spatial frequency white noise and this open loop white noise. Indeed, any type of temporal spectrum is possible for the NEA signal on the lenslets, and as the lenslet noise is decorrelated from a lenslet to another, all noise distribution have the same probability over the lenslet array, therefore the spatial noise distribution is white whatever the lenslet noise statistics. In other words, it is the independence of the noise from a lenslet to another which enable the separation of the noise s-PSD from the noise temporal PSD.
We discuss now the closed loop NEA variance. Let us consider a classical model of the loop architecture: a WFS with an integration time ∆t, followed by a delay t lag due to the WFS readout time and the command computation time, then a numerical integral controller with gain g loop , a digital to analog converter, and the DM. The noise rejection temporal transfer function, defined as the ratio between the NEA signal and the residual error, is given, in the steady state, by (see Demerle et al. [26] for a detailed discussion), where ν is the temporal frequency. The noise power transfer function is given by the modulus square of H n , and we find The closed loop NEA variance σ 2 NEA, CL is given by the integral of the filtered open loop temporal power spectrum, and as the later is a white noise limited to the domain |ν| < 1/(2∆t), we find The open loop NEA variance σ 2 NEA, OL depends on the number of NGS photons received per lenslets during the WFS integration time, the WFS geometry (lenslet width), the WFS integration time, the WFS detector read noise, and the NGS image size, which is tilt-compensated for a closed loop system. Several models have been developed in the literature for this term, and will not be reproduced here -see for instance Rousset [27] and Thomas et al. [28].
Our closed loop phase s-PSD model is now complete, and is given by the sum of the s-PSD of the four fundamental AO errors: the high frequency WFS error, the aniso-servo error, the WFS aliasing error, and the WFS noise error, Let us now illustrate the usefulness and usage of the synthetic model with a few examples.

ILLUSTRATIVE EXAMPLES
The synthetic method has been coded into our AO modeling code PAOLA 7 , a general purpose IDL-based toolbox for modeling the AO correction of segmented telescope static and optical turbulence aberrations. It includes open and closed loop single NGS mode, and a complete multiple NGS ground layer AO mode. We present in this section several studies undertaken with PAOLA, to illustrate the usefulness and usage of such a synthetic tool. An optical turbulence profile and standard telescope and AO system parameters are defined in the Tables 1 and 2 Tables 1 and 2). The wavefront RMS error for the four classical components for this example are given in Table 3. The s-PSD profile and the corresponding PSF profile are shown in Figure 3. It is well known [29,30], as can be seen with this example, that the PSF wings structure mimics the PSD shape. Figure 4 shows the four fundamental errors s-PSD in the spatial frequency plane. The noise and aniso-servo errors affects the central parts of the low frequency domain. Anisotropy of the aniso-servo error is a combined consequence of the wind direction and the off-axis location of the NGS. WFS aliasing, as it is expected, affects the highest spatial frequencies of the low frequency domain. One can therefore expect, in general, the noise errors to contribute essentially to a widening of the PSF core, the aniso-servo error to affect the PSF wings in the region between the core and the transition to the residual seeing halo (due to the high frequency error), while aliasing would affect mostly the transition region. In this example, noise clearly dominates the PSF structure, though.
A note on the spatial frequency pixel size We have seen that both aniso-servo and aliasing s-PSD equations include cosine functions of products in f · v and f · θ. These cosine terms need to be well sampled in the spatial frequency domain when building the numerical matrices f x and f y : an under-sampling would lead to an underestimate of the wavefront error variance, as the later is estimated from numerical integration of the s-PSD, and an incorrect representation of the PSF wings structures. The consequence of such an under-sampling is an over-optimist estimate of the Strehl ratio for large off-axis NGS angle (the Strehl would saturate above a certain value while it should absolutely converge to zero for larger and larger off-axis angles). The same is true for the servo-lag error, where the performance would be over-estimated for large WFS integration time and/or high wind speed. Practically, our experience with PAOLA shows that the cosine terms should be sampled with at least 10 samples over one period. Nyquist sampling is by far not sufficient here.

Strehl ratio of the four fundamental wavefront errors
In this section we simply illustrate how the Strehl associated with the four fundamental errors varies with the main system parameter associated to each error: the WFS lenslet pitch for the WFS high frequency and WFS aliasing error ( Figure 5), the NGS off-axis angle for the anisoplanatic error ( Figure 6, left), and the NGS magnitude for the WFS noise error (Figure 6, right). It is worth noting that these curves were built in only a couple of seconds of CPU time (iMac computer, 3.06 GHz Intel Core 2 Duo processor). We tested also the Maréchal approximation, stating that for low to moderate phase variance σ 2 ϕ the Strehl ratio is given by S = exp(−σ 2 ϕ ). See the dashed curves in Figures 5 and 6. We find that this approximation is actually excellent for the high frequency and WFS aliasing error, relatively good for the WFS noise error, and more questionable for the anisoplanatic error -below a Strehl ratio of about 40% in the example given here.

Open loop versus closed loop performance
We claimed in the introduction that open and closed loop systems behave differently in dim NGS conditions, and limiting magnitude might be quite different for both modes. This claim came from the realization that the noise transfer function are quite different in the two cases, as well as the servo-lag error transfer functions. In order to illustrate this, we computed the Strehl at 1.25 microns for our standard conditions, and a NGS  Initially, we set the WFS integration time fixed at 1 ms, for both modes, and a closed loop gain of 0.5. See Figure 7. We find that the servo-lag error is higher for the closed loop mode, because the rejection transfer function has a lower bandwidth in closed loop than in open loop. The noise error on the other hand is higher in open loop, and this is because the noise is basically unfiltered in open loop, while the noise transfer function is a low pass in closed loop, filtering the high temporal frequency of the white noise spectrum. For given wind conditions, an open loop system can be run faster than a closed loop system because the servo-lag error is intrinsically lower. Therefore, we might think that a dimer NGS could be used in open loop. This is actually true only for bright NGS, where the increased open loop noise error is still low with respect to the servo-lag error. One can see for instance in Figure 7, right, that for a Strehl specification of 0.9 (as it would be for an  As a final experiment, we optimized the WFS integration time and the closed loop gain, for each value of the NGS magnitude. See Figure 8. Optimization has several consequences. First, the limiting magnitude gain is very significant, more than two magnitudes in this example. Second, it makes the open and closed loop servo and noise error converge: this is explained by the fact that the structure of the servo-lag and noise s-PSD are the same in both modes (see Figure 4), therefore optimization converge towards the same solution.
Note that the open loop advantage for bright stars disappears with optimization: the closed loop mode performs the same as the open loop mode at any magnitude. Therefore, contrary to what the intuition would tell, from a control efficiency point-of-view, we assert that there is no advantage of using an open loop rather than a closed loop scheme in NGS AO mode.

CONCLUSION
This paper presents a synthetic modeling method for closed loop astronomical adaptive optics, and complements earlier work on open loop modeling using the same approach. The concept of the synthetic method and its complementarity with the more classical endto-end modeling approach is discussed extensively. The main advantages of the method is that it allows a rapid and direct modeling of the long exposure PSF without going through a long and cumbersome Monte-Carlo process Then, we give the detailed analytical calculation of the spatial power spectrum of the residual AO corrected phase, as well as the steps to go from the power spectrum to the long exposure PSF, allowing the reader to write his/her own modeling code. Dispersion of the air refractive index is included in the model, as well as the deformable mirror spatial transfer function. This method has been implemented into our AO modeling toolbox PAOLA, and used to study a few illustrative examples of the usage of the synthetic method to explore the performance of closed loop AO systems. It is found for instance that when optimizing the WFS integration time, open loop and closed loop system have basically the same performance (same limiting magnitude). Finally, it is important to recall that the foundations of the method do not depend on the type of WFS neither on the type of wavefront reconstruction method, or control algorithm. Also, the method is in principle not limited to single NGS case but can be extended, as it has been done by others, to multi-NGS and multi-DM modes. In this paper, though, we have simply considered the case of an AO system with a single NGS, for a Shack-Hartmann type WFS, a classical LSE wavefront reconstruction and a simple integrator control.