Standard Deviation of Fiber-Coupling Efficiency for Free-Space Optical Communication Through Atmospheric Turbulence

In fiber-based optical receivers such as coherent receiver, the optical signal reaching the receiver's aperture is coupled to the optical fiber through a set of lenses. Due to atmospheric turbulence, however, the coupling efficiency, which is defined by the ratio of the fiber-coupled optical power to the optical power reaching the aperture, fluctuates over time, and thus should be treated as a random variable. Previous theoretical works predict the average coupling efficiency accurately, but fail to estimate the standard deviation, especially when the turbulence is strong. In this study, we propose new theoretical formulas for the standard deviation of coupling efficiency applicable to a wide range of turbulence conditions. For this purpose, we derive a new approximated formula for the second moment of coupling efficiency using the second-order Taylor expansion. Next, a new formula for the second moment of coupled optical power is developed by introducing the spatial filter to the fourth-order coherence function. The computer simulation is carried out to evaluate the accuracy of the formulas. The results show that our formulas agree with the simulation results over a wide range of turbulence strength. Also, the new formulas provide better accuracy than the previously reported ones.


I. INTRODUCTION
D UE to the tremendous increase in demand for wireless data traffic, free-space optical communications (FSOCs) have attracted a great deal of attention as a major alternative technology to conventional radio-frequency counterpart [1]. Since the FSOC systems operate at very high frequencies (e.g., 350 THz), they have several benefits including ultra-wide bandwidth, small antenna size, high antenna gain, low power consumption, and unlicensed spectrum. However, high carrier frequency of optical signals makes them experience severe amplitude and phase distortions when they propagate through atmospheric channel. It is well known that a laser beam propagating through atmospheric turbulence suffers from scintillation, beam wandering, angle-of-arrival fluctuation, and beam spreading. Thus, modeling the beam propagation through atmospheric turbulence has long been a major subject of research for the performance estimation of FSOC systems [2]. Recently, there have been substantial efforts to improve the performance of FSOC systems by utilizing highly sensitive receivers such as optically pre-amplified receivers or coherent receivers [3], [4], [5], [6]. This is because the performance of most FSOC systems is determined by the received optical power as long as the transmitter and receiver provide sufficient bandwidth to generate and detect the signal, respectively. In these receivers, the optical signal reaching the receiver's aperture is first focused by a set of optical lenses and then coupled to the single-mode fiber (SMF). However, since the amplitude and phase of optical signal are distorted by atmospheric turbulence, the coupling efficiency of these fiber-based receivers fluctuates over time, and as a result, makes the coupling efficiency a random variable. To estimate the performance of FSOC systems employing the fiber-based receivers theoretically, it is essential to understand the statistical characteristics of the coupling efficiency.
In [7], Dikmelik and Davidson derived a formula for the average coupling efficiency using the second-order coherence function of incident light and the back-propagated mode field distribution. Subsequently, theoretical analyses on the average coupling efficiency under various conditions have followed, including for non-Kolmogorov turbulence [8], angle-of-arrival fluctuations [9], multi-mode fiber [10], aberration and positional errors of the lens [11].
Together with the average, the standard deviation of the coupling efficiency should be identified to estimate the performance of FSOC systems using fiber-based receivers. Due to random nature of atmospheric turbulence, the coupling efficiency varies widely over time. The link outage occurs when the optical power coupled to the fiber is lower than the receiver sensitivity. Bit-error ratio (BER) performance also deteriorates when the coupling efficiency becomes low instantaneously. Thus, to predict the link outage probability and BER performance of FSOC systems, it is crucial to obtain precise expressions about the higher-order statistics, such as the standard deviation of coupling efficiency. For instance, if the accurate values of the average and standard deviation of coupling efficiency are given, we can estimate the performance of FSOC systems theoretically through the probability distribution function of the detected optical power at the receiver, without exhaustive and time-consuming measurements and/or computer simulations. Unlike the average, the standard deviation of coupling efficiency is extremely difficult to calculate since the fourth-order coherence function is required. The fourth-order coherence function of the light passing through atmospheric turbulence is immensely complex in comparison with the second-order coherence function, and thus it is difficult to develop a consistent theory covering a wide range of turbulence conditions. For these reasons, previous theoretical studies on the standard deviation of coupling efficiency are rare and their limitations are clear. In [12], the fourth-order coherence function was approximated as an expression of the second-order coherence functions with the Gaussian Schell-model. In [13], Ma et al. derived the probability distribution of coupling efficiency as a modified Rician distribution using the statistical optics. However, the results of these models do not agree well with the simulation results. In particular, due to the assumptions used to derive the formulas, the standard deviation of coupling efficiency converges rapidly to its average as the turbulence becomes strong.
On the other hand, X. Ke et al. numerically calculated the relative fluctuations of the optical power coupled to fiber for the Maksutov-Cassegrain telescope system, using the fourth-order coherence function derived from the Rytov perturbation theory [14]. However, the Rytov theory for the fourth-order coherence function is valid only under weak turbulence condition, and cannot be extended to the strong turbulence region without modification. Furthermore, they did not compare these numerical results with simulation nor experimental results. Thus, the scope of application for this theory is not identified.
In this article, we propose new theoretical formulas for the standard deviation of the coupling efficiency applicable to a wide range of turbulence conditions. For this purpose, we make an approximation of the second moment of coupling efficiency using the second-order Taylor expansion. We then derive a formula for the second moment of coupled optical power by introducing the spatial filter function to the fourth-order coherence function. This filter function serves to prevent the second moment of coupled optical power from diverging to infinity under strong turbulence conditions. To validate the developed theory, we carry out computer simulation by solving the stochastic wave equation and compare the results with the theory. The split-step Fourier method is employed to solve the wave propagation through phase screens and free space. The results show that our theory agrees with the simulation results over a wide range of turbulence conditions.

II. THEORETICAL ANALYSIS
As the optical signal passes through the atmospheric channel, the phase of the light is distorted by the random spatial distribution of refractive index of the channel induced by atmospheric turbulence. In fiber-based optical receivers, the optical signal reaching the aperture is coupled to the optical fiber through a coupling lens, as shown in Fig. 1. Since the SMF supports only the HE 11 mode, the modal component that matches to this fundamental mode can be coupled to the SMF.
The optical power coupled to the fiber (P c ) can be calculated by the overlap integral between the incident light and the back-propagated mode-field distribution of SMF [7]. Thus, we have (1) where r i denotes the position vector on the aperture , and E represents the field distribution of the incident optical signal. Also, M is the normalized back-propagated mode-field distribution of SMF at the aperture plane, which is given by where w 0 is the mode-field radius of the fundamental mode of SMF, λ is the wavelength, and f is the focal length of the coupling lens. The total optical power at the aperture plane (P a ) is simply obtained by P a = E (r)E * (r)d 2 r. Then, the coupling efficiency (η) is defined by [7]

A. Average Coupling Efficiency
The average coupling efficiency is written by [7] If we assume a plane wave at the transmitter and its intensity is normalized to 1, P a can be written as the area of the aperture, i.e., S = πD 2 /4, where D is the diameter of the circular aperture. Using (1), we can write the average coupling efficiency as is the second-order coherence function of incident light. From the Rytov theory, Γ 2 can be expressed for a normalized plane wave as [2] where k is the wave number, L is the propagation distance, and is the modified atmospheric spectrum for the refractive-index fluctuation, which is given by [2] Φ n (κ) where l 0 and L 0 are the inner and outer scales of turbulence, and C 2 n is the refractive index structure constant. It should be noted that the coupling efficiency depends upon the coupling geometry through a design parameter a = D 2 πw 0 λf [7]. Thus, we set the focal length of the coupling lens to satisfy a = 1.121 throughout our work. In this condition, we have a maximum average coupling efficiency of 0.8145 in the absence of atmospheric turbulence.

B. Standard Deviation of Coupling Efficiency 1) Second-Order Taylor Approximation:
To obtain the standard deviation of the coupling efficiency, the second moment of the coupling efficiency is needed together with the average value. In [12], the second moment of coupling efficiency was approximated as η 2 ≈ P 2 c / P 2 a . This approximation was based on two assumptions. First, P 2 c and P 2 a are statistically independent. Second, the mean of the reciprocal of P 2 a . can be approximated as the reciprocal of its mean. Under these assumptions, the standard deviation of η was approximated as However, this previously reported approximation exhibits significant errors compared to simulation results, especially when atmospheric turbulence is strong. In this work, we propose a new approximation of the standard deviation of η which improves its accuracy across a wide range of turbulence range. Using the second-order Taylor expansion with expansion point ( P c , P a ) and (4), we can express the variance of η as where σ P c and σ P a are the standard deviations of the coupled power and total power at the aperture, respectively. Also, cov[P a , P c ] is the covariance between P a and P c . If the strength of the turbulence is strong enough so that P c << P a , (9) can be simplified as We will assess the accuracy of this newly derived formula about the standard deviation of η in Section III. It is worth noting that by using this new formula, the second moment of coupling efficiency can be expressed as η 2 ≈ P 2 c / P a 2 .
2) Spatial Filter Function: To estimate the standard deviation of η (i.e., σ η ), it is necessary to obtain the second moment of P c , P 2 c . It is given by where Γ 4 (r 1 , r 2 , r 3 , is the fourth-order coherence function of incident light. Under weak turbulence conditions, Γ 4 can be expressed for a normalized plane wave as [2] Γ 4 (r 1 , r 2 , r 3 , r 4 where The fourth-order coherence function can also be written as Here, where r 0 denotes the Fried's parameter which is given by r 0 = (0.42C 2 n k 2 L) −3/5 .
A grave problem of this formula is that it cannot be applied to the strong turbulence region. This is because as turbulence becomes stronger, the productive term r −5/3 0 included in the exponential function in (15) increases rapidly, whereas F is independent of r 0 . Fig. 2 shows the values of P 2 c obtained by using (11) and (15) as a function of D/r 0 . The ratio D/r 0 is used to indicate the strength of atmospheric turbulence [15], [16]. Most previously reported theoretical works on the coupling efficiency also employ this ratio [12], [13], [14]. The figure shows clearly that the value of P 2 c diverges to infinity as the strength of turbulence gets higher. Also, due to the productive term (1 − k κ 2 L sin κ 2 L k ) contained in the integrand of (16), this tendency becomes more pronounced as the propagation distance L increases. This problem is originated from the fact that (12) is derived from the perturbation theory assuming weak turbulence.
To solve this problem, we propose to utilize the extended Rytov theory and the spatial filter function proposed by Andrews et al. [2]. The extended Rytov theory is a modification of the Rytov approximation that allows us to extend the results of P 2 c obtained numerically by using (11) and (15). D = 10 cm, λ = 1550 nm, L 0 = 2 m, and l 0 = 10 mm. The intensity of initial plane wave is normalized to 1.
weak turbulence region to moderate-to-strong turbulence. In the extended Rytov theory, the light passing through turbulence is distorted by eddies of two different scales. Large-scale eddies mainly distort the phase and wavefront of light, and consequently they produce the refractive effects such as beam wandering and defocusing. On the other hand, small-scale eddies bring about the diffractive effects of light, which are related to beam spreading and/or intensity fluctuations. In order to reflect these effects, the eddies should satisfy certain conditions: small-scale eddies should be smaller than the coherence radius ρ 0 = (1.46C 2 n k 2 L) −3/5 ; large-scale eddies should be larger than the scattering disk L/kρ 0 . As a result, in the moderateto-strong turbulence region, the influence of eddies having the size between these two scales is suppressed. This phenomenon can be expressed by applying a filter to the spatial power spectrum of refractive index. In [2], Andrews et al. proposed the effective atmospheric spectrum Φ n,e (κ) to analyze the scintillation index of beam under the strong turbulence condition: where G is the spatial filter function .
(18) Here, f (κl 0 ) and g(κL 0 ) represent inner and outer scale filter modifications, which are defined by respectively. In these equations, κ X is the large-scale spatial frequency cutoff and κ Y is the small-scale spatial frequency cutoff. Those cutoff frequencies satisfy the following conditions By using the spatial filter G, we can make only the effects of the low spatial frequencies (κ < κ X ) and the high spatial frequencies (κ > κ Y ) contribute to the theory.
In this work, we make the following assumptions to suppress the divergence of Γ 4 under moderate-to-strong turbulence: 1) The spatial filter function G of (18) is included in F.
2) The second-order coherence function Γ 2 in (15) is given by (6). These assumptions are based on the observation that the values of P 2 c obtained by using (15) agree well with the simulation results when the effect of the term F is completely suppressed (See Appendix). Based on the assumptions, the fourth-order coherence function Γ 4 included in (11) is now replaced by where F m includes the spatial filter function G given by (18) and is given by r 1 , r 2 , r 3 , r 4

III. VERIFICATION OF THEORY THROUGH SIMULATION
In this section, we evaluate the accuracy of our formulas by comparing the theoretical values with the simulation results. Computer simulations are widely used to analyze the characteristics of the atmospheric channel [17], [18], [19]. We solve the stochastic wave equation using the split-step Fourier method. Phase screens are first generated by using the fast Fourier transform (FFT) method and periodically placed along the propagation path. For the atmospheric channel with a propagation distance L, the power spectrum of a phase screen is given by where κ r = κ 2 x + κ 2 y is the angular spatial frequency. The phase screens are implemented by generating two-dimensional random phase distribution having a power spectrum of (24). In the FFT method, the phase at the position (x, y) of the phase screen is obtained by sampling as follows where h(κ x , κ y ) is a Hermitian complex Gaussian white noise process with a zero mean and a unit variance. The size of the phase screen is A×A. In the split-step Fourier method, the propagation link from the transmitter to the receiver is divided into N PS sections. Each section has a length of Δz = L/N P S and a phase screen generated by the FFT method is located at the center of each section. These phase screens have a Fried's parameter of r 0,Δz = (0.42C 2 n k 2 Δz) −3/5 . The sections between the phase screens reflect the diffraction effect of light in free space through the Huygens-Fresnel integral. By solving the wave equation at every section and repeating it to the end of propagation distance, we can calculate the beam propagation through atmospheric turbulence.
Reference [18] shows the following criteria should be satisfied for the simulation using phase screens: (1) For the phase screens to reflect the effects of the outer and inner scales, the size of the phase screens should be sufficiently larger than the outer scale, and the pixels of each phase screen should be sufficiently smaller than the inner scale. (2) To satisfy the Nyquist sampling theorem, the phase difference between neighboring pixels of the phase screen should be less than π. (3) To prevent the aliasing in the propagation of the beam through the Fourier transform, the distance between the phase screens should satisfy Δz < k π AΔx , where Δx ( = Δy ) is the pixel size of phase screens. (4) To ignore the edge effect of the phase screens, the size of the phase screens must satisfy A > 2θ max Δz for the maximum scattering angle, θ max = max[ 1 k (∂φ/∂x) 2 + (∂φ/∂y) 2 ]. We set the simulation parameters in our work to satisfy these criteria. We also make the resolutions of the phase screens higher than the above criteria for accurate numerical calculation of the overlap integral of (1). Throughout the simulation, the wavelength, the mode-field radius of SMF, and the design parameter are set to be 1550 nm, 5 μm, and 1.121, respectively.
We first evaluate the accuracy of the second-order Taylor approximation in (10). For this purpose, we define the error between the theoretical and simulation values as where σ η is the standard deviation of instantaneous coupling efficiency and σ η denotes the estimated standard deviation formulated in (8) or (10). Fig. 3 shows the errors of σ η as a function of D/r 0 under a couple of conditions. In order to isolate the results of the Taylor approximation from the introduction of spatial filter function, we utilize the samples of P a and P c obtained from simulation to calculate σ η . In this way, we can estimate the standard deviation formulated in (8) or (10) without using the fourth-order coherence function. The results show that when the turbulence is stronger than a certain level (e.g., D/r 0 > 10), the errors of the new approximation [i.e., (10)] are smaller than those of the previously reported one [i.e., (8)]. In this range, the errors of (10) falls mostly within 10%. Thus, we can say that (10) is a better approximation than (8) under strong turbulence conditions. However, our new approximation of (10) is slightly inaccurate compared to (8) when the turbulence is weak. For example, when D, r 0 and L are 10 cm, 2 cm, and 5 km, respectively, the error is calculated to be 18% for our new formula, but it is 2% for the previously reported approximation. Nevertheless, the new approximation exhibits decent accuracy over a wide range of turbulence strength.
Next, we evaluate the efficacy of the introduction of spatial filter function to the fourth-order coherence function for the estimation of σ η . In FSOC links using fiber-based receivers, the aperture diameter D and the propagation distance L are important parameters. In addition, the outer scale L 0 and the inner scale l 0 are the parameters that specify the turbulence characteristics. Therefore, we investigate the behavior of the errors for the proposed formula (22) as a function of D/r 0 by varying D, L, L 0 and l 0 . The numerical values of theoretical formulas are obtained from the quasi-Monte Carlo integration of 20 million iterations using the Sobol sequence [20]. Fig. 4 shows the errors of the standard deviation of coupling efficiency obtained by theoretical formulas for various parameters. We can see that when the aperture diameter is large enough  (e.g., D ≥ 6 cm), the errors of the proposed formula are less than 20%, except under extremely strong turbulence conditions (e.g., D/r 0 > 40). Although the errors of new formulas are relatively large when the aperture diameter is small, they are much lower than the errors obtained without using the spatial filter under the entire range of turbulence strength we tested. Meanwhile, the results show that the effects of L, L 0 and l 0 on the accuracy of our formulas are not pronounced. Overall, our new formulas estimate the standard deviation of coupling efficiency even when D/r 0 is as large as 40. Therefore, we conclude that our new formula provide accurate estimation of the standard deviation of fiber coupling efficiency over a wide range of turbulence conditions.

IV. CONCLUSION
In this article, we propose new formulas for the standard deviation of fiber coupling efficiency applicable to a wide range of turbulence conditions. We utilize the second-order Tayler expansion to obtain the approximation of the second momentum of the coupling efficiency. We also introduce the spatial filter to the fourth-order coherence function to avoid this coherence function from diverging to infinity as the turbulence gets stronger. Computer simulations are carried out to validate the newly developed formulas. The results show that the formulas provide better accuracy than the previously reported ones over a wide range of turbulence conditions. In particular, the standard deviations of the fiber coupling efficiency obtained from our proposed formulas agree well with the simulation results even if D/r 0 is as large as 40. Thus, we believe that the new formulas could be used to theoretically predict the performance of high-speed free-space optical communication systems utilizing optically pre-amplified receivers or coherent optical receivers.

APPENDIX
The purpose of introducing the spatial filter G in Section II is to avoid the values of Γ 4 formulated by (15) from diverging to infinity. When the propagation distance is very short, the value of F converges to zero, and thus P 2 c can be approximated as P 2 c ≈ Γ 2 (r 1 , r 2 )Γ 2 (r 2 , r 3 )Γ 2 (r 3 , r 4 )Γ 2 (r 4 , r 1 ) Γ 2 (r 1 , r 3 )Γ 2 (r 2 , r 4 ) ×M * (r 1 ) M (r 2 ) M * (r 3 ) M (r 4 ) d 2 r 1 d 2 r 2 d 2 r 3 d 2 r 4 (27) Fig. 5 shows the values of P 2 c obtained from the simulation and (27) when the propagation distance is 10 meters. It shows clearly that the approximation of (27) agrees very well with the simulation over a wide range of D/r 0 . This exemplary comparison indicates that the theoretical formula for P 2 c should converge to (27) as L goes to zero. Actually, this approximation agrees very well with the simulation under short propagation distance conditions, regardless of parameters such as D, L 0 , and l 0 . Therefore, it is reasonable to assume that the spatial filter function applies only to F [or inside the exponent of the exponential function in (15)], not to the spectrum included in the second-order coherence function. Fig. 6 shows the comparison between the theory and simulation when the propagation distance is 1 km. Other parameters are listed as follows: D = 10 cm, λ = 1550 nm, L 0 = 2 m, and l 0 = 10 mm. By introducing G and thus modifying F as formulated by (23), the theoretical values of P 2 c do not diverge to infinity and match very well with the simulation results. Here, the discrepancy between the simulation results and the theory using the spatial filter function varies depending on the parameters, such as D, L, L 0 , and l 0 . The results depicted in Fig. 4 in Section III show the effects of these parameters.