FAST: Fourier domain adaptive optics simulation tool for bidirectional ground-space optical links through atmospheric turbulence

: Free space optical links between the ground and space may be severely degraded by atmospheric turbulence. Adaptive Optics, a technique allowing partial correction of this degradation, is beginning to see use in the field with the potential to achieve more robust and higher bandwidth links. Here we present a simulation tool, FAST, which utilises an analytical Fourier domain Adaptive Optics model developed for astronomy. Using the reciprocity principle, the simulation may be applied either to downlink post-compensated or uplink pre-compensated beams. We show that FAST gives similar results to full end-to-end simulations with wave-optical propagation whilst being between 10 and 200 times faster, enabling the characterisation of optical links with complex Adaptive Optics systems in timely fashion.


Introduction
The use of light at optical wavelengths as a communications channel has several advantages over traditional radio frequency communications. Increases in bandwidth, security and power efficiency make free space optical communication (FSOC) [1] particularly attractive for Earth-space links between orbiting satellites and ground stations [2][3][4].
Working with visible and near-infrared wavelengths brings with it the limitation of significant signal degradation from atmospheric turbulence. Small random refractive index changes in the air through which the light is propagating cause distortions of the beam. This manifests as both a decrease in the average received power and variability of the power over millisecond timescales. For the communications channel this can result in bit errors, where a transmitted 0 is received as a 1 or vice versa, and signal fades where the link is lost entirely for a short period [5].
Adaptive optics (AO) is a mature technology designed to mitigate the effects of atmospheric turbulence. First developed for astronomy [6,7], it has found applications ranging from microscopy to vision science [8]. The fundamental principle is to measure the turbulence induced wavefront aberration using a wavefront sensor (WFS) and to apply a correction to this aberration in real time, usually with a deformable mirror (DM). For FSOC, AO can be applied both at the transmitter (pre-compensation or pre-distortion) as well at the receiver, and offers the possibility of a more robust, higher bandwidth link [9][10][11].
Modelling of such AO-corrected optical links is of critical importance in order to predict the performance of these systems and assess the potential gain of employing AO. Theory exists for gaussian beams propagating through uncorrected turbulence based on the model of Rytov perturbations, with modifications to account for beam wander [5]. A model of a tracked beam (low order tip/tilt modes removed) may also be derived from the same theory [12]. However, analytical expressions derived from Rytov theory are not sufficient in themselves to model an optical link, since they give the mean and variance (scintillation index) of the received flux but not the shape of the flux probability distribution. Several probability distribution models exist for the received flux that are applicable in different scenarios however they must be assumed a priori and do not necessarily employ parameters with a physical interpretation [5,13,14].
Depending on the turbulence strength and beam sizes, low-order turbulence (LOT) theory may provide a better model of beam wander [15], and similarly a tip/tilt corrected model may be derived [16]. Since LOT is by definition limited to low order turbulent modes, higher order AO cannot be modelled.
The most accurate simulation of an AO corrected optical links are obtained with Monte Carlo wave optical (MC-WO) simulation tools, where the turbulent atmosphere is modelled as a series of discrete random phase screens through which the beam propagates [11,17]. Between the screens, wave optical propagation (usually near-field i.e. Fresnel [18]) is used. Many random or temporally correlated turbulence iterations are used to characterise the link. High fidelity models of the various AO components (WFS, reconstructor, DM, etc.) may be included. Both WO propagation and the AO system model can add significant computational complexity leading to long simulation times. Simplified simulations employing a more basic AO model and avoiding WO propagation [19,20] are therefore useful in the cases where the high fidelity of full MC-WO simulation is not required.
In astronomy, many simulations exist that model the performance of AO in the vertical downlink propagation case [21][22][23][24][25][26]. A particularly interesting subset is those that operate entirely in the Fourier domain [27][28][29][30]. By making certain assumptions about the AO-corrected residual phase, the power spectral density (PSD) of this phase in the telescope pupil may be computed analytically. From the PSD, the long exposure point spread function (PSF) can be calculated with only a single Fast Fourier Transform (FFT). This is particularly attractive since it reduces computation times by several orders of magnitude when compared with full Monte Carlo simulation [31].
Here, we apply the Fourier domain AO simulation methods of [28,32] to the FSOC case. FAST (Fourier domain Adaptive optics Simulation Tool) takes a semi-analytical approach similar to [19,33], the difference being the computation of residual Monte Carlo phase screens directly from the analytically computed residual phase PSD. The Fourier AO model includes five fundamental single conjugate AO errors (DM fitting, anisoplanatism, servo-lag, WFS aliasing, WFS noise) and can operate in closed loop. Although here we consider a Shack-Hartmann WFS system, in principle different WFSs such as pyramid sensors can also be simulated [34]. Turbulence outer and inner scale effects are included in the PSD, and temporal sequences can be generated by translating phase screens as in standard MC AO simulation.
We show that FAST gives similar intensity distributions and fade statistics to MC-WO simulation in both downlink and uplink (via the reciprocity principle [33]) ground-space links to GEO, whilst being 10 times faster for the downlink and 200 times faster per iteration in the uplink cases. We believe this simulation will be a useful intermediate point between limited analytical formulae and computationally expensive MC-WO simulation.

Simulation fundamentals
Simulations for astronomical AO are designed for the vertical or slant path downward propagation regime, where the light from the distant astronomical source may be modelled as an infinite plane wave. Atmospheric turbulence, which extends to around altitude h = 20 km above the Earth's surface, is therefore considered to be close to the receiver.
In satellite-to-ground (downlink) FSOC, a laser beam with gaussian profile is launched from the satellite and received by a telescope at an optical ground station (OGS). The launched beam is usually no more than several centimetres in width due to mass and space requirements on the satellite. The altitude of low earth orbit (LEO) and geostationary orbit (GEO) satellites (hundreds to tens of thousands of kilometres) and the small beam size mean that the downlink beam at the top of the atmosphere can be considered an infinite plane wave to good approximation [5], and hence the standard astronomical AO downlink simulations are applicable.
For the ground-to-satellite (uplink) case, the turbulence is concentrated near the transmitter, which can result in significant beam wander and speckle at the satellite depending on the size of the launched beam [35]. The modelling of such a link is therefore fundamentally different to the downlink case. An additional challenge is that the transmitter must aim ahead of the satellite in order to illuminate it, due to the light travel time between the transmitter and receiver and the angular velocity of the satellite as seen from the ground. For GEO, this point-ahead angle (PAA) is of the order of a 4 arcseconds, whereas for LEO it can be up to 10 arcseconds. This poses a problem for AO if the downlink beam is used for wavefront sensing since the two paths may encounter different turbulence.
However, the principle of far-field reciprocity [33,36] implies that, regardless of the turbulence distribution along the path, if the roles of transmitter and receiver are interchanged, the fractional power transfer between the two is identical. In this case, the two reciprocal planes are the focal plane of the ground based receiver and the far-field pupil plane of the satellite (see Fig. 1). Therefore we can use the downlink simulation for the uplink, with a change of aperture function from a fully illuminated aperture to a truncated gaussian beam or any arbitrary launch beam profile. AO pre-compensation is then equivalent to AO post-compensation of the downlink beam, i.e. the standard astronomical case. The following will therefore be focused only on the downlink case with the knowledge that it can be applied to the uplink using the reciprocity principle. In the downlink case, a small gaussian beam is launched from the satellite (indicated with red arrows) that propagates downwards towards the receiver. The beam may be approximated as a plane wave at the top of the atmosphere due to the large divergence. The light is captured by the receiver and focused onto the focal plane, where it can be coupled into an optical fibre. In the uplink case, a gaussian beam is launched through the aperture and propagates upwards through turbulence. The turbulence close to the transmitter in this case can produce significant beam wander and speckle. The light is received by the satellite and similarly coupled into an optical fibre or other detector (not shown in this diagram). Importantly, the far-field reciprocity principle implies that the focal plane of the ground-based telescope and the satellite pupil plane, indicated by dashed blue lines, are reciprocal provided that the distance to the satellite is large enough. This means that we can model the uplink scenario using a much simpler downlink simulation.
In addition it is important to note that the theory in this section is predicated on the fact that the turbulent perturbations to the wavefront are small, i.e. that we are in the weak turbulence regime. This is usually true in the vertical propagation regime but may not hold at very small elevation angles (<10 degrees).

Definitions
We model the turbulent atmosphere as a number of independent turbulent layers at altitudes h i , with thickness δh and layer strengths given by where J i is the turbulent strength of the ith layer in units of m 1/3 and C 2 n (h) is the refractive index structure constant as a function of altitude h, with units of m −2/3 . It is assumed that δh is effectively infinitesimally small compared to the distance between layers such that diffraction within the layer is negligible [37].
Light propagates through the turbulent layers to a receiver at the ground (h = 0). The electric field in the plane of the receiver aperture is denoted where χ and ϕ denote the field log-amplitude and phase respectively. As a result of atmospheric turbulence and the AO system, both χ and ϕ are gaussian random variables with 0 mean that vary in both position r and time t. We assume that χ and ϕ are independent variables, which as we will see will allow for a faster simulation. This assumption holds in the vertical propagation regime [19,38]. The spatial power spectral density (PSD) of these random variables describes the variance power in Fourier space. They are defined as the modulus square of the Fourier transforms of ϕ and χ, averaged over time where f = (f x , f y ) denotes the two dimensional spatial frequency vector and ⟨⟩ denotes the time average. Note that we are using linear spatial frequencies, i.e. the units of f are 1/length [37]. For clarity, the forward and inverse Fourier transforms used here have the explicit form for an arbitrary function g. For a plane wave propagating through a turbulent layer at altitude h i in the small perturbation regime, the PSDs of ϕ and χ, denoted Φ φ and Φ χ , are related to the PSD of refractive index where λ is the wavelength considered, the wave vector k = 2π/λ and f = |f| [37]. We take the commonly used geometric approximation for the phase, i.e. that cos 2 (πλh i f 2 ) ≈ 1 in Eq. (7), since we are simulating a regime where h i is relatively small. Assuming a Von Karman model for the turbulent refractive index fluctuations, the PSD is given by where f m and f 0 are spatial frequency cutoffs defined as for inner and outer scale values l 0 and L 0 [5].

AO corrected phase PSD
The foundation of the Fourier approach is to assume that the AO system acts as a series of spatial filters on the phase power spectrum (Eq. (7)). After correction, the PSD of the residual phase is given by where each term corresponds to the PSD of the fundamental error terms of classical AO. Φ fit φ denotes the DM fitting PSD, describing residual error from the finite sampling of the DM and therefore inability to correct high spatial frequencies. Φ aniso-servo φ (f) denotes the joint anisoplanatism-servo lag PSD, which describes the errors resulting from both an angular separation of WFS target and science target and the temporal separation between measurement of the wavefront and application of the correction. These errors are coupled, therefore the joint PSD is computed. The WFS aliasing PSD Φ alias φ (f) describes the error from the finite sampling of the WFS: high spatial frequencies are seen as low spatial frequencies and therefore erroneously corrected by the system. Finally, the WFS noise PSD Φ noise φ (f) describes the effect of random noise on the WFS signals and how these propagate through the reconstruction.
Analytical expressions for the terms in Eq. (12) are derived in [28,32]. We will not enter into detailed discussion of these derivations here. The general idea is to define a region of smaller spatial frequencies in which the AO system controls the phase. The shape and size of this region depends on the subaperture geometry and correction type (i.e. modal or zonal correction). The aniso-servo, aliasing and noise PSDs are defined only within this region, with the fitting PSD defined in the complementary high spatial frequency region. We show an example of each of these PSDs in Fig. 2.

Focal plane
These focal plane computations follow the work of [19,33], and are presented here for completeness.
The primary value of interest for FSOC is the power ρ coupled into the detector, which most commonly takes the form of a single mode optical fibre in the focal plane. This may be modelled for an instant in time t in the pupil plane as where P(r) represents the receiver pupil function and w(r) the gaussian fibre mode projected into the pupil plane, of the form where A is chosen such that |︁ |︁ ∬ P(r)w(r) dr |︁ |︁ 2 = 1 and w 0 is found by optimising the coupled power in the diffraction limited case. The instantaneous coupling ρ(t) is therefore the instantaneous Fig. 2. Example of the four PSD terms that make up the residual phase PSD. From left to right: fitting, aniso-servo, WFS aliasing, WFS noise. In each case the AO correction region is marked by the red dashed line, which in this case is a square indicating correction on a zonal basis with square geometry. Note that the fitting error PSD (left) is presented at a larger spatial scale to the others, illustrating that is defined only in the high spatial frequency region. Each PSD has a different logarithmic colour scale for clarity, where black represents large values and white small values. electric field Ψ(r, t), with phase and log-amplitude components, integrated over the aperture and weighted by the projection of the fibre mode. For the uplink case, the "focal plane" being computed is actually the reciprocal plane of the satellite in LEO/GEO. The satellite receiver, with aperture diameter usually up to tens of centimetres, is much smaller than the coherence scale of intensity fluctuations and so is effectively point-like. The flux coupled into the satellite receiver therefore depends only on the intensity at the central point of the (reciprocal) focal plane plus an additional constant coupling loss depending on the design of the satellite receiver that we do not compute. This central intensity can be computed by removing the fibre mode from Eq. (13) (i.e. setting w(r) = 1), in which case it is equivalent to the on-axis Fourier transform |F [P(r)Ψ(r)](0)| 2 .
We wish to compute not only the mean value ⟨ρ⟩ but also its distribution, in order to obtain useful quantities for FSOC. To this end, we take a Monte Carlo approach and generate random realisations of Ψ(r) in order to obtain random realisations of ρ.
As previously stated, we make the assumption that χ(r, t) and ϕ(r, t) are independent normally distributed random variables, which allows us to treat their contributions to the coupling separately, i.e. ρ(t) = ρ χ (t)ρ φ (t).
Under weak fluctuation conditions, the variance of the aperture averaged log-amplitude, which we denote σ 2 χ , can be obtained via aperture filtering of the log-amplitude PSD (Eq. (8)) [39,40] where we have summed the contributions from each layer at altitude h i , and the aperture averaged log-amplitudeχ. The spatial filter |F [P(r)w(r)](f)| 2 is easily obtained for an arbitrary aperture via a numerical fast Fourier transform (FFT). The contribution to the coupling is then where randomχ(t) are obtained by drawing samples from a normal distribution with mean µ = 0 and variance σ 2 = σ 2 χ .
The contribution to the coupling from the phase is then given by therefore we generate random instances of ϕ(r) with the correct statistics and perform the integration numerically. These random phase screens can be efficiently generated from the (residual) phase PSD via an inverse FFT of spatially filtered normally distributed random variables [41]. If required, we can also ensure the variance at lower spatial frequencies is correct by adding subharmonics [42]. This makes the simulation similar to full end to end Monte Carlo simulations of AO, but by modelling the AO system in the Fourier domain and avoiding computationally expensive wave-optical propagation we gain a significant advantage in speed.

Temporal correlation
For some applications it is required to obtain temporal sequences of ρ(t), for example to assess the frequency and other statistics of signal fades. To accomplish this, for the log-amplitude we generate a temporal PSD by assuming Taylor frozen flow [43] and integrating the appropriately scaled spatial PSD (Eq. (8)) along the axis orthogonal to the wind direction for each turbulent layer [44]. This PSD is used to temporally filter the randomly generated log-amplitude sequences [19].
For the phase, we have confirmed with end-to-end AO simulation that the Taylor frozen flow approximation applies to all residual phase errors apart from noise. For the fitting, aniso-servo and aliasing we generate only a single random screen per layer and translate them according to the wind speed and direction as in conventional Monte Carlo AO simulation. This is accomplished in Fourier space by applying a spatial shift to the initial screen, i.e.
where ∆t is the simulation timestep and v = (v x , v y ) is the wind speed vector for the layer. Since we generate one screen per layer in this mode the computation time is affected, increasing approximately linearly with the number of layers. Residual phase from noise error is unchanged and uses uncorrelated random samples. When generating temporal sequences one must ensure the phase screens are large enough (or equivalently, the total simulated time short enough) such that the translated screens do not wrap around, resulting in the same phase aberration appearing multiple times and introducing a periodicity to ρ(t). For the FSOC case, we are sampling fades on the order of milliseconds, so it is sufficient to simulate many very short time periods. This has the additional bonus of decreasing the simulation convergence time.

Advantages of technique
The speed advantage over WO simulation provides several benefits. First, parameters of both the AO system (e.g. subaperture size, loop rate) and the atmosphere (e.g. r 0 , turbulence profile) may be scanned quickly to assess how these parameters impact the quality of the optical link. This is useful in the early phases of a specifying a system to rapidly constrain the design.
The increase in number of iterations per unit time also allows better characterisation of the coupled flux distribution, since more random iterations can be computed in the same amount of time. This is particularly important when assessing the likelihood of low probability deep signal fades.
In conventional Monte Carlo AO simulation, any closed loop control system or simulation incorporating temporal effects (loop rate, loop delay, exposure time etc.) must operate using temporal sequences of turbulence. This increases the convergence time of the simulation as each iteration is correlated with previous iterations as opposed to being completely random. In analytical AO simulation there is no such limitation, as random phase screens may be computed from the closed loop residual PSD. This results in much faster simulation convergence in these cases.

Limitations of technique
The most limiting assumption of the Fourier technique is that the statistics of the residual phase are stationary across the pupil. Effectively, an infinite aperture is assumed. In reality the finite size of the DM and WFS mean the phase at the edges of the pupil tends to be poorly controlled, leading to an increase in residual phase variance. This also means that links where the gaussian nature of the beam is important (larger apertures, shorter propagation distances) cannot be simulated accurately since we assume the phase variance is constant across the pupil which is not true in general for gaussian beams [5].
Another significant limiting factor is that the analysis is founded on the assumptions of weak turbulence, meaning FAST will not be accurate in strong turbulence conditions. This limits the simulation to above around 30 degrees elevation angle depending on the turbulence conditions. Many satellite FSOC links will need to operate as low as 10 degrees elevation especially for LEO, and such links may not currently be simulated with FAST.
The modelling of WFS noise is more difficult with an analytical simulation. In end-to-end Monte Carlo simulation the impact of shot noise, read-out noise, etc. on the WFS measurements are straightforwardly obtained since the WFS (and by extension, the detector) is modelled in real space. For the analytical simulation we must use a model to convert between input flux and measurement noise. For a Shack Hartmann WFS several models exist, depending on the centroiding method [45], but must necessarily include some assumptions about the profile of the subaperture spot, so are limited to certain regimes. Additionally, there is currently no method of including the effect of scintillation on the WFS noise in stronger turbulence.

Comparison with end-to-end Monte Carlo wave-optical simulation
Here we present a comparison between results obtained with FAST and full MC-WO AO simulations. These simulations are soapy [21], an end-to-end Monte Carlo AO simulation for astronomy (downlink), and the code used in Osborn et al [11], a Monte Carlo simulation of precompensated ground-space optical links with a simple AO model. Both simulations employ wave-optical propagation between turbulent phase screens. We simulate two scenarios, uplink and downlink communications links with a satellite in GEO. We assume that the transmitter/reciever on the satellite is unresolved, which allows us to make two simplifications. First, in the downlink case the beam at the top of the atmosphere may be approximated by a plane wave. Second, in the uplink case the receiver is a point at the centre of the far-field, therefore the reciprocity theorem applies (see Fig. 1). The general simulation parameters for each case are described in Table 1.
We present coupled flux distributions, defined as the loss in dB relative to the diffraction limited coupling ρ dl : ρ[dB] = 10 log 10 (ρ/ρ dl ).
We also present fade statistics in both cases. The probability of fade is related to the cumulative distribution of ρ, normalised to the average ⟨ρ⟩ [5]. It is simply computed from the random samples P(ρ<I T ) = n(ρ<I T )/n tot , where n(ρ<I T ) is the number of samples with value less than the threshold I T and n tot is the total number of random samples. The threshold I T has units of dB and is computed as I T [dB] = 10 log 10 (ρ/⟨ρ⟩). Table 1. Atmospheric and optical link parameters for the downlink and uplink comparison simulations. Atmospheric parameters (r 0 , θ 0 , τ 0 , Rytov variance) are defined along the 30 degree elevation path. Differences in these parameters between the uplink and downlink are only due to the different wavelengths used; the atmospheric conditions are the same in both cases.

Downlink Uplink
Turbulence profile Hufnagel Valley 5/7 [11] Wind profile Bufton [46] Outer Scale (L 0 ) 25 m MC-WO simulation soapy [21] Osborn et al [11] The mean fade duration is calculated by computing the duration of each fade in the random samples (providing they are temporally correlated), and taking the mean of the result.
For the uplink and downlink scenarios we simulate three AO configurations: uncorrected turbulence (no AO), tip/tilt correction and full AO, with full parameters defined in Table 2.

Downlink
For the downlink case the obtained coupled flux distributions are shown in Fig. 3. We see for the most part good agreement between the MC-WO simulation and FAST for all AO configurations. It is clear for this set of parameters that tip/tilt correction only does not provide much improvement over the uncorrected case, primarily due to the relatively large aperture (D = 1 m). The full AO system however provides a marked improvement. This is continued in Fig. 4 where we compare fade statistics, specifically the probability of fade and mean fade duration. Again we observe good agreement between the MC-WO simulation and FAST. Since the fade threshold I T (Eq. (22)) is computed relative to the mean of the distribution, the uncorrected and tip/tilt lines are almost identical. For this downlink scenario, temporally correlated samples were used in both simulations, allowing the fade duration to be computed which is shown in the lower panel of Fig. 4.

Uplink
For the uplink case we present the coupled flux distributions in Fig. 5. This configuration offers quite different distributions to the downlink case, as the aperture used is smaller (12 cm gaussian beam radius launched on a 48 cm diameter clear aperture) making tip/tilt correction more effective. Reciprocally we may consider the correction of tip/tilt as correction of beam wander on the uplink, which is known to dominate the intensity variability for these beam sizes [5]. For uplink in this scenario there is only a slight advantage for full AO. We observe good agreement in all configurations between the simulations.
In Fig. 6 we show the fade probability for the three AO configurations. This shows good agreement for the uncorrected and full AO cases, with an overestimation of the fade probability for deep fades in the tip/tilt case. This is most likely because the Fourier model is poorer representing low-order modes leading to a discrepancy in this case. This could likely be alleviated by a change of pixel scale or by including a greater number of subharmonics in the phase screen generation.

Speed comparison
In Table 3 we show some timing information for the simulations used above. All simulations are benchmarked on a server running an Intel Xeon Gold 6244 CPU with 32 cores. For each simulation we measure the time taken to initialise the simulation (setup time) and the time per iteration once it is running. We can see that FAST offers approximately an order of magnitude speed increase compared to soapy in the downlink case. It should be noted that soapy is highly optimised, and that many of these optimisations are also available to FAST but have yet to be implemented.
For the uplink, FAST is almost 200 times faster than the simulations presented in Osborn et al [11]. These uplink simulations are particularly slow due to the sampling requirements of gaussian beam propagation, i.e. large phase screens are required. Generally since FAST does not perform numerical optical propagation, sampling requirements are not as strict leading to a significant speed increase.

Conclusion
We have presented a semi-analytical simulation, FAST, that is capable of modelling the performance of an AO corrected bidirectional ground-space optical link at speeds between 10 and 200 times faster than conventional MC-WO simulations. This is accomplished firstly by using a Fourier domain AO model, which allows us to analytically compute the PSD of the residual phase. This PSD is used to create phase screens which allow the rapid Monte Carlo characterisation of the optical link.
We also avoid wave-optical propagation, the cost of which is the simulation functions only in weak turbulence conditions. However, despite the various limitations and assumptions of FAST, we have shown that it gives consistent results with MC-WO simulation for both uplink pre-compensated and downlink post-compensated beams for Rytov variances as high as 0.46. The scenario simulated here is representative of a GEO feeder link and is therefore of direct relevance to FSOC.
The simulation occupies the middle ground between theoretical formulae with limited application to the AO-corrected case, and full end-to-end MC-WO simulation. Since a Fourier domain AO model is used, many parameters (e.g. DM pitch, loop rate, WFS noise) can be modified and their effect on the optical link quantified using FAST. At the same time, the simplified AO and propagation models allow these parameter spaces to be rapidly spanned due to the short simulation times. We therefore believe that FAST will be a useful tool for the growing AO/FSOC community, and the code is available at https://github.com/ojdf/fast.