Atmospheric effects on continuous-variable quantum key distribution

Compared to fiber continuous-variable quantum key distribution (CVQKD), atmospheric link offers the possibility of a broader geographical coverage and more flexible transmission.However, there are many negative features of the atmospheric channel that will reduce the achievable secret key rate, such as beam extinction and a variety of turbulence effects. Here we show how these factors affect performance of CVQKD, by considering our newly derived key rate formulas for fading channels,which involves detection imperfections, thus form a transmission model for CVQKD. This model can help evaluate the feasibility of experiment scheme in practical applications.We found that performance deterioration of horizontal link within the boundary layer is primarily caused by transmittance fluctuations (including beam wandering, broadening, deformation, and scintillation),while transmittance change due to pulse broadening under weak turbulence is negligible. Besides,we also found that communication interruptions can also cause a perceptible key rate reduction when the transmission distance is longer, while phase excess noise due to arrival time fluctuations requires new compensation techniques to reduce it to a negligible level. Furthermore, it is found that performing homodyne detection enables longer transmission distances, whereas heterodyne allows higher achievable key rate over short distances.


Introduction
Nowadays quantum key distribution (QKD) [1] through atmospheric turbulence channel over long distance has been realized [2][3][4], and satellite-to-ground discrete variable QKD [1] over 1200 km has been verified recently [5]. However, systems using single-photon detectors suffer from background noise [6], while homodyne or heterodyne detection with bright local oscillator (LO) acting as a filter can reduce the background noise [7]. Experiments measuring Stokes operators [8][9][10][11][12][13] have shown the filtering effect of LO. Besides, quantum-limited coherent measurements between a geostationary Earth orbit satellite and a ground station has been conducted [14]. Nevertheless, there is still no complete Gaussian-modulated coherent state (GMCS) CVQKD [15,16] experiments being reported. Therefore, for future experiments and applications, it is quite necessary to analyze the atmospheric effects on GMCS CVQKD.
Recently, an elliptical beam model [17,18] considering beam wandering, broadening and deformation has been established for quantum light through the atmospheric channel. The states of entanglement-based CVQKD through fading channels have been deduced, and the secret key rate without considering detection efficiency and noise has also been calculated [19]. However, the detection efficiency and noise have significant impacts on the final achievable key rate, and the atmospheric effects on signal transmission are not only the three effects included in the elliptical beam model but also many other effects such as arrival time fluctuations, temporal pulse broadening, angle-of-arrival fluctuations and scintillation [20]. Therefore, a comprehensive transmission model and corresponding performance analysis of CVQKD in the atmosphere are necessary.
In this report we consider GMCS CV-QKD horizontal link on the surface of the Earth. We deduce a new achievable secret key rate for the atmospheric channel of CVQKD scheme with imperfect homodyne and

Secret key rate through atmospheric channels
To investigate CVQKD in the atmospheric channel, we first analyze the secret key rate through fading (fluctuating) channels. The description of entanglement-based GMCS CVQKD over the fading channel is shown in figure 1. Alice and Bob share an entangled state generated by the EPR source with variance V. One mode of the entangled state, B 0 is transmitted to Bob through a fading channel characterized by a distribution of transmittance T, and Bob performs homodyne or heterodyne detection to measure field quadratures. The imperfection of the detector is described by detection efficiency η and electronic noise el u contained in variance ν.
In the asymptotic regime, the secret key rate K is given as [23] K I , 1 where β is the reconciliation efficiency, I AB is the Shannon mutual information of Alice and Bob, and χ BE is the Holevo quantity, which can be expressed as [24] S mpm S d , 2 where m B represents the measurement of Bob, p m B ( ) represents the probability density of the measurement, m E B r represents the eavesdropper's state conditional on Bob's measurement result, and S (·) represents the Von Neumann entropy. To calculate I AB and χ BE , we first need to find the covariance matrix after fluctuating channels. The covariance matrix of a two-mode squeezed vacuum state generated by the EPR source is given as )is the Pauli matrix. After a channel with excess noise ε and random variable transmittance T, the covariance matrix becomes [19] It can be seen from equation (4) that the influence of the fading channel is primarily reflected in T á ñ and T á ñ. Thus, considering the detection efficiency η and electronic noise el u , we can obtain the mutual information of Alice and Bob for homodyne detection . Nevertheless, the five symplectic eigenvalues in [24] can not be used directly. Thus, we deduce the symplectic eigenvalues of both homodyne and heterodyne detection for fading channels (for details see appendix A). It is noteworthy that the results presented in appendix A are not only applicable to the atmospheric turbulence channel but also to other channels whose transmittance changes randomly, such as underwater channels.
However, for the atmospheric channel which may cause angle-of-arrival fluctuations, equation (1) should be revised to where P stands for interruption probability due to angle-of-arrival fluctuations. The details of angle-of-arrival fluctuations for equation (8) will be demonstrated in section 3.4.

Atmospheric channel effects on CVQKD
The secret key rate demonstrated in equation (8) indicates that parameters (T, ε and P) closely related to atmospheric effects should be analyzed in depth so that we can approximately assess the performance of atmospheric CVQKD by developing the method proposed in [19]. In this section, some well-developed models in atmospheric optical communications will be employed to accomplish the assessment of performance, in which we will build links between the models and atmospheric CVQKD, and show how much influence would be. Besides, a new phase excess noise caused by pulse arrival fluctuations will be derived. Atmospheric channel effects primarily include beam extinction and turbulence effects. On the one hand, extinction is caused by absorption and scattering by molecules and aerosol which leads to attenuation of light intensity. On the other hand, random variations in the refractive index of atmosphere may cause pulse temporal broadening, transmittance fluctuations (signal fading), communication interruptions, and extra phase excess noise.
In this paper, we consider that (temporal and spatial) Gaussian beam is transmitted horizontally on the surface of the Earth, as depicted in figure 2. Time and polarization multiplexed LO and signal pulses are generated by the GMCS module and collimated by the telescope. Subsequently, the pulses that undergo turbulence and extinction, are collected by a telescope and measured by coherent detector, i.e., homodyne or heterodyne. It is noteworthy that LO can also be generated in Bob's terminal [25,26]. This 'local' LO scheme is capable of avoiding loopholes due to sending LO through the channel, but not mature relative to the scheme of simultaneous transmission of LO and signal, which has been developed over 15 years since the first experiment [16]. In order to integrate with optical fiber systems, the wavelength of laser is chosen as 1550 nm which is also in the atmospheric transmission window. The Rytov variance is employed to describe the strength of turbulence which is given by [22] C k L 1. where k=2π/λ is the optical wave number, λ is the wavelength of light, L is the horizontal propagation distance, and C n 2 is the index of refraction structure parameter. Many models of C n 2 have been proposed over the years [27]. However, since the link is assumed to be located within the boundary layer, it is rather reasonable to assume C n 2 (median) to be constant as shown in table 1, which are based on results of long-term radiosonde measurement conducted in Hefei, Anhui, China [28]. Now let us estimate the impact of atmospheric effects on the performance of CVQKD.

Transmittance: beam extinction
For CVQKD, beam extinction means that the transmittance associated with wavelength and propagation path decreases as transmission distance increases. For horizontal paths, the transmittance is given by [21] T L e , 10 L ext where the total extinction coefficient α(λ) comprises the aerosol scattering, aerosol absorption, molecular scattering, and molecular absorption terms: There are some models that can be used to estimate the transmittance of a particular environmental conditions for CVQKD, such as LOWTRAN, MODTRAN and FASCODE [21]. It is assumed that the horizontal link is located in the mid-latitude countryside and has a visibility of 23 km, then the extinction coefficients can be estimated by LBLRTM [28]. Since C n 2 in spring and summer are close to each other, here we consider the worse one, summer. Besides, the strongest turbulence occurs in winter, thus the case of winter should be considered. The extinction coefficients are listed in table 2.

Transmittance: temporal pulse broadening
In this section we will study the transmittance change due to temporal pulse broadening under weak turbulence, regardless of the beam extinction.
Temporal pulse broadening is primarily caused by two mechanisms [29]: first, the difference in arrival time of each individual pulse when only single scattering is affecting the pulse, i.e., pulse arrival time fluctuations (pulse wandering) which is also responsible for extra excess noise and will be further discussed in section 3.5, second, the pulse spreading brought by multiple scattering process of each pulse. The combination of these two mechanisms causes temporal pulse broadening, which can be described by the averaged pulse intensity I L t r, ; á ñ ( ), or referred to as mean irradiance. Without loss of generality, here we consider the temporal pulse broadening of Gaussian pulse [20], whose intensity has a shape of I t t T exp is the pulse half-width. Here, R dut and f rep are the duty ratio and pulse recurrence frequency (PRF), respectively. The temporal pulse broadening of Gaussian pulse is demonstrated in figure 3. The temporal width of pulse is broadened by the atmosphere, and the pulse intensity is also attenuated. The free-space irradiance of a collimated beam under the near-field (Ω?1) and far-field (Ω=1) approximations is given by [30] respectively, where c is the light speed in free space, ω is the angular frequency of the light, W 0 is the beam-spot radius at the transmitter and kW L 2 0 2 W = is the Fresnel parameter. However, since the amount of spreading and arrival time of each pulse are different, equations (13) and (14) are not able to be applied to characterize broadening in turbulence. This is exactly why ) is required. Under near-field assumption the mean irradiance in weak turbulence is given by [31] It can be seen from equations (15)-(17) that the turbulence-induced temporal pulse broadening in both nearfield and far-field approximations can be characterized by T 1 . Here, we define the pulse broadening ratio as which only in the femtosecond order will have a significant impact, as indicated by figure 4. The outer scale on the ground is assumed to be 0.4 m [33]. The result is for winter whose turbulence strength is strongest. Now we consider the transmittance change due to temporal pulse broadening. Comparing equations (13) and (15) shows that pulse broadening will result in a T 0 /T 1 -fold attenuation of the average light intensity of the received signal, as shown in figure 3. The same result can be found by comparing equation (14) with equation (17). Thus, the mean value of pulse broadening introduced transmittance can be expressed as The mean transmittance in winter is demonstrated in figure 5, T bro á ñ varies quickly from the femtosecond level to the picosecond level, but after the picosecond level, T bro á ñ is approximately equal to one. Therefore, T bro á ñ is also approximately equal to one (for details see appendix B). In other words, the transmittance introduced by pulse broadening is actually negligible in the regime of weak turbulence. However, in the regime of strong turbulence, the analysis of equation (18) requires a large amount of numerical calculations [34]. In [34], the results also show that broadening is only perceptible on the order of femtosecond in strong turbulence, i.e., the broadening-induced transmittance approaches one, thus negligible. Therefore, there is no need to consider the transmittance change caused by pulse broadening in the following performance analysis in section 4.

Transmittance: beam wandering, broadening, deformation and scintillation
In this subsection we will concentrate on transmittance fluctuations (signal fading) which is primarily caused by beam wandering, beam broadening, beam deformation, and scintillation. The elliptical beam model [17] well describes beam wandering, broadening and deformation in weak and strong turbulence, as shown in figure 6. However, note that the moderate-to-strong transition regime of this model has not been clarified yet, the corresponding performance analysis of CVQKD in this regime may need a better transmittance model. Weak, moderate and strong turbulence correspond to   Here we define ζ=f−j 0 , the transmittance is then given approximately by [17] where a is the receiving aperture radius, and specific expressions of the other parameters are shown in appendix C. The transmittance T ell is a function of five real parameters, x y , , , , In the isotropic turbulence case, f can be seen as a uniformly distributed random variable, having no correlations with other parameters. Considering Based on equation (19), appendix C and table D1, the probability distribution of T ell can be estimated by Monte Carlo simulations. The transmittance probability density function (PDF) of T ell in summer on the ground is shown in figure 7, with extinction involved. The receiving aperture radius a and the initial beam-spot radius W 0 in figure 7 are assumed to be 110 mm and 80 mm, respectively. Now we consider transmittance fluctuations due to scintillation which is not incorporated in the elliptical beam model. The total transmittance is considered as multiplication of equation (19) and transmittance due to  scintillations, which approximately gives a lower bound for atmospheric CVQKD, so as to make a conservative estimation of performance of CVQKD. Inherently, turbulence effects contained in the elliptical model and scintillations should be combined together, and this will be further investigated in the next step of our work.
The scintillation effect is illustrated by the light intensity spatial distribution of the beam cross-section (see the schematic diagram in figure 8). A great deal of turbulent vortices contained in the cross section independently scatters and diffracts the portion of the light impinging thereon, such that the intensity of light at each spatial point (irradiance) in the cross section varies randomly.
Over the years, many irradiance PDF models have been proposed to characterize the randomly fading irradiance signal, such as lognormal distribution, K distribution, I−K distribution, lognormal-Rician distribution, and gamma-gamma distribution [22].
where r is a transverse vector, á ñ · is an ensemble average, L r, ) is the scintillation index (for details see appendix E), and I L r, á ñ ( ) is the (normalized) mean irradiance (for details see appendix G). Considering largescale and small-scale effects, the (normalized) irradiance in strong fluctuation can be described by gammagamma distribution [22] p I I K I 2 2 , 21 where Γ(·) is the gamma function, K α−β is the modified Bessel function of the second kind, α is the effective number of large scale cells of the scattering process, and β is the effective number of small scale cells. Both α and β are related to the scintillation index, and detailed in appendix F. Nevertheless, the distribution of irradiance only describes the intensity fluctuations at a certain spatial point. Hence, the irradiance should be integrated within the plane of the receiving aperture  : Since I L r, ( )is a random variable, it is quite difficult to calculate equation (23) directly. Here, we still apply Monte Calro simulations to estimate T sci . For simplicity, only the scintillation and beam broadening reflected in I L r, á ñ ( ) are considered, regardless of beam wandering and deformation. In this case, the scintillation-induced transmittance fluctuations is quite small, and the reduction in transmittance is primarily caused by the beam broadening, as shown in figure 9. Since the scintillation-induced transmittance fluctuation is too small, the PDF in figure 9 looks just like a line. Therefore, the PDF of transmittance at L=10 km in the inset of figure 9 is presented as an example. The PDF shape of L=5 and 15 km are the same as the shape of L=2 km.

Interruption probability
Due to the high directivity of laser transmissions, there exists the possibility of communication interruption when there is a large angle-of-arrival fluctuation. The direct reflection of angle-of-arrival fluctuations on the receiving aperture plane is image jitter on a focal plane. When the focus is not within the receiving fiber core, communication is interrupted at this time (see figure 10).
The interruption phenomenon is closely related to beam wandering comprised in the elliptical model. However, the major concern of the elliptical beam model is the total energy collected by RL, i.e., truncation of beam spot. It is still possible that part of the beam spot is within RL while the focus displaced out of FC, as illustrated in figure 11. At this point, if only the elliptical beam model is involved, the signal transmission will still be considered as successful, but in fact not. Thus, the consideration of communication interruption is necessary. Due to the close relationship between beam wandering and interruption, the statistics of them should be the same, that is, β a is Gaussian distributed.
Assuming the mean value of arriving angle 0 a b á ñ = and β a is small enough so that sin a a b b @ , the variance of β a can be written as where f is the focal length of the collecting lens.
Since β a is normally distributed, the focus displacement is also normally distributed. Thus, the communication interruption probability can be expressed as where d cor is the diameter of the fiber core in meters. A typical single-mode optical fiber has a core diameter from 8.3 to 10.5 μm. Here, as an illustration, we assume that the core diameter is 9 μm, f=220 mm, a=110 mm, and W 0 =80 mm. The interruption probability is shown in figure 12.

Excess noise
In this subsection we will focus on pulse arrival time fluctuations observed by a fixed observer (see figure 13). This effect may bring extra excess noise by causing phase fluctuations. The pulse arrival time between the LO and signal Δt is now a random variable. To clarify Δt, we first investigate the arrival time of a single pulse t a . The mean value and variance of t a are given by [31] t , , Lt t r r , ; , ; d 28 is the nth moment with the complex envelope v L t r, ; 0 ( ). Under weak turbulence, near-field and far-field approximations, the mean value and on-axis variance of arrival time is given by [20]   where T 1 is given in equation (16). The on-axis variance of strong turbulence is deduced by [43]. Now we define the random variable Δt as where ω is the angular frequency of light mentioned in equation (14). When the phase fluctuation is small enough, the excess noise can be expressed as [25] V , 3 4 where V A =V−1 is the modulation variance of Alice. Figure 14 shows the excess noise varies with distance L. Here we consider that f PRF =100 MHz, R dut =10%, 1 10 ta 13 r = --, and V A =2 in shot noise unit (SNU). With such a high correlation coefficient and weak turbulence condition, phase excess noise can eventually reach an acceptable level, which is hard to achieve in practice. The phase compensation method for fiber CVQKD [35] can compensate small phase fluctuations,  but it is not applicable to atmosphere CVQKD whose phase fluctuations is very large. Therefore, we hope that an effective phase compensation method for atmospheric CVQKD will be proposed in the future. It is also worth noting that the turbulence-induced phase excess noise may be more readily to be decreased, if a non-pilot aided 'local' LO scheme, which does not require reference pulses or pilot, e.g. the experiment [14], is successfully applied in GMCS CVQKD. This is because only the arrival time of the signal is fluctuant, whereas the LO is not transmitted through the channel. Unfortunately, in this scheme, part of signal would be split off to lock phase. This operation would bring extra noise and decrease total detection efficiency in a complete experiment. Therefore, it would be a trade-off between the simultaneous transmission and 'local' LO scheme.

Performance analysis
In this section, we will combine the results of sections 2 and 3 to analyze the achievable final key rate.
The Monte Carlo method is applied to estimate the secret key rate in equation (8). Since the excess noise caused by phase fluctuations can not yet be accurately compensated, it is quite difficult to estimate the practical excess noise after using different effective compensation methods. Thus, here we do not consider the phase excess noise that changes with atmospheric conditions for the time being, but we still examine the achievable key rate under different fixed excess noise level, namely ε=0.01 and ε=0.03 in SNU. And as explained at the end of section 3.2, the temporal pulse broadening will not be considered in this section. The extinction coefficients used are listed in table 2 and all other parameters needed in performance analysis are presented in table 3.
The secret key rate with ε=0.01 SNU is demonstrated in figure 15(a). As we can see, the transmission distance of system using homodyne detection is longer than heterodyne detection, but not much. Now we increase the excess noise to ε=0.03 SNU, as depicted in figure 15(b). Compared with figure 15(a), the achievable transmission distance of is obviously shorter.
The performance analysis conducted in this section indicates several key points. First, homodyne detection should be applied to practical systems if farther transmission distance is demanded, otherwise, heterodyne provides higher achievable key rate at short distance. Here, especially, the heterodyne case deserves attentions, since only the no-switching protocol [36] had been proven to against general attacks in a realistic finite-size regime [37,38]. Second, we found that the transmittance fluctuations are destructive to the key rate. Accordingly, in practical experiments, the main effort should be devoted to inhibiting the effects of beam wandering, broadening and deformation. Third, since the impact of excess noise is quite significant, and the phase excess noise would be much more than 0.03, effective methods of controlling phase excess noise are needed to increase the secret key rate. It is noteworthy that our performance analysis was conducted on the assumption that there is no enhancement technique, such as adaptive optics [14,39] and post selection [19]. Adaptive optics are shown to be cable of relieving signal fading [40], correcting wavefront and improving fiber coupling efficiency [41]. The incoming result of post selection is, actually, increasing T á ñ by discarding data when instantaneous transmittance is too low. These methods or techniques can directly or indirectly weaken atmospheric effects. Therefore, the performance of practical systems may be better than the results shown in this paper. We expect that phase compensation techniques for atmospheric CVQKD can also be proposed so that the phase excess noise can be reduced to a negligible level.

Conclusion
We analyzed atmospheric effects on the horizontal CVQKD links on the Earth's surface, thus establishing a transmission model, which can help the performance assessment of practical CVQKD systems. The newly derived key rate formulas for fading channels with detection efficiency and noise considered shows that there are three main parameters that affect the final key rate: the transmittance, interruption probability and excess noise. The transmittance change caused by temporal pulse broadening under weak turbulence regime is found to be negligible. Transmittance fluctuations caused by beam wandering, broadening, deformation, and scintillation make the final key rate deteriorated rapidly. Angle-of-arrival fluctuations may cause communication interruptions which leads to a more obvious decline in the key rate over long-distance transmission. The phase excess noise caused by pulse arrival time fluctuation is found to be quite large. We found, in fading channels, systems employing homodyne detection can transmit far more distances than heterodyne detection, while employment of heterodyne detection at short-range transmission has a higher key rate than homodyne detection.
The symplectic eigenvalues λ 1,2 can be calculated for both homodyne and heterodyne detection by ) is the variance of T . Then, λ 3,4,5 are the symplectic eigenvalues of covariance matrix m AFG B g , which can be expressed as  Appendix B. The estimation of T bro á ñ In this appendix we estimate the value of T bro á ñ when T 1 bro á ñ @ .
Although both T bro and T bro are random variables, the relationship between them can be determined, as shown in figure B1. Since T 1 bro á ñ @ and T bro 1, only the area around T bro ≅1 needs to be considered. We can see from the inset in figure B1 that when T 1 bro @ , the relationship between T bro and T bro is approximately linear, namely,

Appendix D. Mean values and covariance matrix elements of v
The mean values and covariance matrix elements of v is shown in These mean values and elements are given for horizontal links. Figure B1. The relationship between T bro and T bro . The inset shows the details of the relationship from T bro =0.9 to 1. r 0 , , , where F a b x ; ; 1 1 ( )is the confluent hypergeometric function, and the longitudinal component is given as