Satellite-mediated quantum atmospheric links

TL;DR

This paper analyzes how Earth's atmosphere affects satellite-based quantum communication links, focusing on signal distortions, transmittance fluctuations, and security implications for quantum protocols.

Contribution

It provides a detailed geometric and statistical model of atmospheric effects on satellite quantum links, aiding secure quantum communication design.

Findings

Transmittance fluctuations saturate at large zenith angles.

Atmospheric effects strongly depend on zenith angle.

Derived transmittance probability distribution for security analysis.

Abstract

The establishment of quantum communication links over a global scale is enabled by satellite nodes. We examine the influence of Earth's atmosphere on the performance of quantum optical communication channels with emphasis on the downlink scenario. We derive the geometrical path length between a moving low Earth orbit satellite and an optical ground station as a function of the ground observer's zenith angle, his geographical latitude, and the meridian inclination angle of the satellite. We show that the signal distortions due to regular atmospheric refraction, atmospheric absorption, and turbulence have a strong dependence on the zenith angle. The observed saturation of transmittance fluctuations for large zenith angles is explained. The probability distribution of the transmittance for slant propagation paths is derived, which enables us to perform the security analysis of decoy state protocols implemented via satellite-mediated links.

Tables

Table 1. Table 1: Atmospheric, optical beam, and geographical parameters used in the simulations. Zenith-dependent values are given at Z a = 0 ∘ subscript 𝑍 𝑎 superscript 0 Z_{a}=0^{\circ} .

Parameter	Notation	Value
Radius of spatial coherence	${\rho }_0$ (cm)	$13$
Reference structure parameter	$C_{n,0}^2$ (m-2/3)	${10}^{-17}$
Wavelength	$\lambda$ (nm)	$840$
Initial beam-spot radius	$W_0$ (cm)	$2$
Receiver’s aperture radius	$a$ (m)	$0.5$
Wave-front radius	$F$ (m)	${10}^5$
Observer’s coordinates		${48}^{\circ }$N, ${11.5}^{\circ }$ E

Table 2. Table 2: Overview of conducted measurement campaigns.

Year	Satellites/	Wavelength (nm)	No. Measurements
	(Laser terminal)
2006	OICETS (LUCE)	847	5
2009	OICETS (LUCE)	847	4
2015	ISS (OPALS)	1550	1
2016	ISS (OPALS)	1550	1
2016	Socrates (SOTA)	1550	1

Table 3. Table 3: Parameters for a QKD system.

Parameter	Notation	Value
Detector efficiency	${\eta }_{\det }$	$60\%$
Extinction due to optics	${\chi }_{\mathrm{opt}}$	$84\%$
Background yield (dark count)	$Y_0^{\mathrm{DC}}$	$5.89\times {10}^{-7}$
Erroneous detector prob.	$e_{\det }$	$1\%$
Background error rate	$e_0$	$50\%$
Failure probability	$\epsilon$	${10}^{-5}$
Mean intensity of signal	${\mu }_s$	$0.8$
weak decoy state	${\mu }_d$	$0.1$
Number of sent bits	$N$	${10}^{11}$
Pulse repetition rate	$r_N$ (MHz)	$150$
Generation prob. signal bits	$p_s$	$65\%$
weak decoy bits	$p_d$	$25\%$
Generation prob. of x-basis bits	$p_a^x$, $\left(a=s,d\right)$	$60\%$
Error correction efficiency	$f\left(\text{QBER}\right)$	$1.16$
Tracking precision	${\sigma }_{\mathrm{tr}}$ ($\mu$rad)	$1$

Table 4. Table 4: The reference altitudes and values and gradients of the linearly segmented temperature-height and pressure-height profiles from the Earth’s surface up to altitude of 85 km.

Sub-	Height	Temperature	Temperature	Pressure
script		gradient
$b$	$H_b$, (km)	${\lambda }_b$, (K/km)	$T_b$, (K)	$P_b$, (mb)
0	0	-6.5	288	1 013
1	11	0.0	217	226
2	20	+1.0	217	54.7
3	32	+2.8	229	8.68
4	47	0.0	271	1.11
5	51	-2.8	271	0.67
6	71	-2.0	215	0.04
7	84.8	-	188	0.004

Table 5. Table 5: The defined reference levels, gradients, and values of the linearly segmented refractive index-height profiles from surface to 85 km.

Sub-	Height	Refractive index	Refractive index
script		gradient
$i$	$H_i$, (km)	${\left(dn/dh\right)}_i\times {10}^{-6}$, (km-1)	$\left(n_i-1\right)\times {10}^8$
0	0	-	27 340
1	5	25.68	14 660
2	7	17.58	11 142
3	11	12.50	6 141
4	15	7.183	3 268
5	20	3.565	1 485
6	32	1.042	235
7	47	0.134	34
8	51	0.034	21
9	71	0.010	1
10	84.8	0.001	0.1

Equations

$$L(Z)=\sqrt{H^2+2H{R}_{\oplus }+R_{\oplus }^2{\cos }^{2}Z}-R_{\oplus }\cos Z,$$

$$Z_{\min }^{\Delta \iota }=\arccos \left[\sqrt{1-{\cos }^2\Psi {\sin }^2\Delta \iota }\right],$$

$$\Delta \iota =n\frac{T_{\mathrm{sat}}{v}_{\oplus }}{R_{\oplus }},$$

$$Z_a=\arcsin \left(\frac{1}{n_0}\sin Z\right)$$

$$L_i$$

$$\displaystyle-2(R_{\oplus}+H_{i-1})(R_{\oplus}+H_{i})\cos[\Phi(Z,r_{i})]\Bigr{\}}^{1/2}.$$

$$L_{11}$$

$$\displaystyle-2(R_{\oplus}+H_{10})(R_{\oplus}+H)\cos[\Theta(Z,r_{10})]\Bigr{\}}^{1/2}.$$

$${\epsilon }_r(Z_a)=\frac{1}{L(Z_a)}\sum _{i=1}^{11}{L}_i(Z_a),$$

$${\epsilon }_r(Z_a)$$

$$-4.066061\times 10^{-5}|Z_a{|}^3+3.813573\times 10^{-6}{Z}_a^4$$

$$-1.920844\times 10^{-7}|Z_a{|}^5+5.710429\times 10^{-9}{Z}_a^6$$

$$-1.032821\times 10^{-10}|Z_a{|}^7+1.117105\times 10^{-12}{Z}_a^8$$

$$-6.644358\times 10^{-15}|Z_a{|}^9+1.672433\times 10^{-17}{Z}_a^{10}.$$

$$L_r(Z_a)={\epsilon }_r(Z_a)L(Z_a),$$

$$\displaystyle\chi_{\mathrm{ext}}=\exp\Bigl{[}-\int\limits_{0}^{L_{r}(Z_{a})}{\rm{d}}L^{\prime}\beta_{\mathrm{ext}}(L^{\prime})\Bigr{]},$$

$${\beta }_{\mathrm{ext}}(L^{\prime })={\beta }_{\mathrm{ext}}^0\frac{N(L^{\prime })}{N_0}={\beta }_{\mathrm{ext}}^0\exp \left[-L^{\prime }/({\mathcal{H}}_0\sec Z_a)\right],$$

$$\displaystyle\chi_{\mathrm{ext}}=\exp\Bigl{[}-\beta_{\mathrm{ext}}^{0}\frac{\mathcal{H}_{0}\sec Z_{a}}{1000}\left(1-e^{-L_{r}(Z_{a})/(\mathcal{H}_{0}\sec Z_{0})}\right)\Bigr{]}.$$

$$n_i(\mathbf{\rho },t)=n_i+\delta n(\mathbf{\rho },t),$$

$${\Phi }_n(\mathbf{k})=\frac{1}{(2\pi {)}^3}\underset{{\mathbb{R}}^3}{\int }{\mathrm{d}}^3\mathbf{\rho }{B}_n(\mathbf{\rho })e^{-i\mathbf{k}\cdot \mathbf{\rho }},$$

$$B_n({\mathbf{\rho }}_1-{\mathbf{\rho }}_2)$$

$$=⟨\delta n({\mathbf{\rho }}_1)\delta n({\mathbf{\rho }}_2)⟩.$$

$${\Phi }_n(\mathbf{\kappa };z)=\frac{1}{\pi }\underset{0}{\overset{z}{\int }}\mathrm{d}{z}^{\prime }{F}_n(\mathbf{\kappa };z^{\prime }),$$

$$F_n(\mathbf{\kappa };z)=\underset{{\mathbb{R}}^2}{\int }{\mathrm{d}}^2\mathbf{r}{B}_I(\mathbf{r},z)e^{-i\mathbf{\kappa }\cdot \mathbf{r}},$$

$${\Phi }_n(\mathbf{\kappa },z)=0.033C_n^2(z)|\mathbf{\kappa }{|}^{-\frac{11}{3}}.$$

$${\mathcal{D}}_S({\mathbf{r}}_1,{\mathbf{r}}_2,{\mathbf{r}}_1^{\prime },{\mathbf{r}}_2^{\prime };z)=⟨{\left[S({\mathbf{r}}_1,{\mathbf{r}}_1^{\prime };z)-S({\mathbf{r}}_2,{\mathbf{r}}_2^{\prime };z)\right]}^2⟩,$$

$$S^{\mathrm{UL}}(\mathbf{r},{\mathbf{r}}^{\prime };z)=\frac{k}{2}\underset{0}{\overset{z}{\int }}\mathrm{d}{z}^{\prime }\delta n\left(\mathbf{r}\frac{z^{\prime }}{z}+{\mathbf{r}}^{\prime }\frac{z-z^{\prime }}{z},z^{\prime }\right)$$

$$S^{\mathrm{DL}}(\mathbf{r},{\mathbf{r}}^{\prime };z)=\frac{k}{2}\underset{0}{\overset{z}{\int }}\mathrm{d}{z}^{\prime }\delta n\left(\mathbf{r}\frac{z-z^{\prime }}{z}+{\mathbf{r}}^{\prime }\frac{z^{\prime }}{z},L-z^{\prime }\right)$$

$$⟨\delta n(\mathbf{r};z)\delta n({\mathbf{r}}^{\prime };z^{\prime })⟩$$

$$\times {\Phi }_n(\mathbf{\kappa },z)\exp \left[i\mathbf{\kappa }\cdot (\mathbf{r}-{\mathbf{r}}^{\prime })\right],$$

$${\mathcal{D}}_S({\mathbf{r}}_1,{\mathbf{r}}_2,{\mathbf{r}}_1^{\prime },{\mathbf{r}}_2^{\prime };z)={\mathcal{D}}_S({\mathbf{r}}_1-{\mathbf{r}}_2,{\mathbf{r}}_1^{\prime }-{\mathbf{r}}_2^{\prime };z)$$

$$=\pi k^2\underset{0}{\overset{z}{\int }}\mathrm{d}{z}^{\prime }\underset{{\mathbb{R}}^2}{\int }{\mathrm{d}}^2\mathbf{\kappa }{\Phi }_n(\mathbf{\kappa },z-z^{\prime })$$

$$\displaystyle\times\Biggl{\{}1-\exp\Bigl{(}i\boldsymbol{\kappa}{\cdot}\Bigl{[}(\boldsymbol{r}_{1}{-}\boldsymbol{r}_{2})\frac{z-z^{\prime}}{z}+(\boldsymbol{r}_{1}^{\prime}{-}\boldsymbol{r}_{2}^{\prime})\frac{z^{\prime}}{z}\Bigl{]}\Bigr{)}\Biggr{\}}.$$

$${\mathcal{D}}_S$$

$$=2{\rho }_0^{-\frac{5}{3}}\underset{0}{\overset{1}{\int }}\mathrm{d}\xi \frac{C_n^2([1-\xi ]L_r)}{C_{n,0}^2}{\left|\mathbf{r}(1-\xi )+{\mathbf{r}}^{\prime }\xi \right|}^{\frac{5}{3}},$$

$${\rho }_0\approx (1.5C_{n,0}^2{k}^2{L}_{\mathrm{turb}}{)}^{-3/5}$$

Full text

††preprint: PHYSICAL REVIEW A 99, 053830 (2019)

I Introduction

The recent success in practical realization of quantum state transfer between satellites and ground stations Vallone2015 ; Vallone2016 ; Yin2017 ; Yin2017a ; Guenthner2017 ; Liao2017 ; Liao ; Takenaka ; Liao2018 ; Bedington2017 ; Calderaro2018 has established the cornerstone of future global quantum communication. This achievement is especially impressive in the view of how fast the quantum communication over free-space quantum links reached its maturity. During only three decades the free-space quantum communication reached distances of 7600 km Liao2018 , while starting from proof-of-principle demonstration over a 30-cm quantum channel Bennett1989 .

An acceleration of the development of satellite-based quantum technologies is caused by the need to establish global secure communication networks. The present classical public key cryptography is based on mathematical problems that admit no efficient solutions with currently available technologies. However, its security is vulnerable to quantum computer hacking attacks, in particular, those based on Shor’s factoring algorithm Shor1994 . Therefore, the future realization of quantum computers threatens the currently used classical cryptographic protocols to become insecure and hence useless.

Conversely, quantum key distribution (QKD) establishes unconditional secure cryptographic keys between two distinct parties Bennett1984 ; Scarani2009 . The security of QKD is based on some fundamental principles of quantum physics such as the no-cloning theorem Wooters1982 , Bell correlations Ekert1991 , uncertainty relations for the most of continuous-variable protocols Hillery ; Cerf ; Gottesman ; Grosshans2002 ; Grosshans2003 ; Madsen ; Braunstein , etc. Ironically, while providing the security of QKD, the no-cloning theorem at the same time puts limitations on possible communication distances due to inability to amplify a quantum signal. One way to bring quantum communication to global scale is the use of satellites. Satellite-mediated QKD networks could span large communication distances by linking widely separated ground stations Liao2018 ; Bedington2017 .

The setup of the satellite-mediated quantum links is a demanding task. The following problems must be overcome for successfully establishing a feasible quantum communication with satellites: the relative motion of the communication parties Bonato2006 ; Zhang2014 ; Pramanik2017 , the influence of gravity Bruschi2014 ; Bruschi2018 , the clock synchronization problem Robert2016 ; Wang2016 , acquisition, tracking and pointing issues with moving platforms Nauerth2013 ; Wang2013 ; Pugh2017 , the influence of background noise Erlong2005 ; Bonato2009 ; Tomaello2011 , to name just a few. For low Earth orbiting (LEO) satellites, the communication time is limited to a few minutes Bacco2013 ; Casado2017 and this puts additional tight bounds on communication security. Finally, Earth’s atmosphere contributes significantly to the loss budget and to the deterioration of the quantum signal due to diffraction, scattering, absorption, and atmospheric turbulence Tatarskii2016 ; Andrews2005 ; Bourgoin2013 .

The common feature of the aforementioned issues is the strong dependence of disturbing factors on the relative position of the moving satellite and the ground station. For the observer at the ground station the instantaneous position of the satellite is given by the zenith angle, i.e., by the angle between the satellite and the vertical direction. The dependence of atmospheric disturbances on the zenith angle has attracted the attention of scientists since ancient times. Already Ptolemy (100-175 A.D.) noticed that the atmospheric refraction displaces the apparent position of celestial bodies toward the zenith Lehn2005 . Such a regular refraction is caused by the altitude variation of the atmospheric refractive index. It influences the precision of satellite measurements and is responsible for apparent elongation of celestial objects, mirages, green flashes on setting sun, to name just a few Young2006 . Other types of disturbances, such as intensity scintillation Andrews2000 ; Dravins1997 and image dancing Tatarskii2016 ; Chiba1971 , arise from turbulent fluctuations of the atmospheric refractive index. The strength of these fluctuations varies with altitude and hence its influence on optical light propagation is also a function of the zenith angle. Thus, the examination of the zenith angle dependence of relevant characteristics is of immense importance for the establishing of robust and reliable optical quantum communication links via the satellites.

In this article we address several aspects of ground-to-space communication connected with Earth’s atmosphere. We focus on quantum links between optical ground stations (OGS) and LEO satellites. Orbital periods of LEO satellites are of 120 min or less, leading to short passage times of the OGS, with which communication is aimed to be established. Depending on the inclination angle of the satellite orbit, the expected communication time varies in the range of 8-16 min. This communication time window is actually narrower due to the time needed for satellite acquisition. Additionally, satellites with trajectories of small inclination angle could lose the contact with the OGS because of natural barriers, while OGS telescopes could disconnect communication with satellites at small zenith angles due to mechanical construction of telescope mounts. Under these circumstances, it is desirable to start the acquisition procedure when the satellite appears on the horizon. In this case, quantum communication could be established already at small satellite elevation angles (large zenith angles). Since the disturbing influence of the atmosphere grows with the growth of the zenith angle, the zenith angle dependence is crucial for the strict analysis of quantum satellite communication links.

In the following we provide a systematic analysis of quantum communication channels with the inclusion of disturbance effects due to geometrical elongation of communication links, atmospheric regular refraction, extinction, and turbulence. In our consideration we focus on the zenith angle dependence of the associated signal losses. Additionally, the finite communication time puts severe restrictions on secure quantum key length in satellite-mediated communication. Careful estimation of the loss budget of quantum channels is needed to estimate the lower bound of the key length and its rate. In this article we perform the security analysis of a decoy state protocol that is commonly used in practice.

This article is structured as follows. In Sec. II we review some basic communication scenarios with satellites and in Sec. III.1 we obtain the geometrical distance between a ground station and a satellite. We focus on atmospheric refraction effects in Sec. III.2. The regular extinction of the optical signal due to scattering and absorption on air particles is discussed in Sec. III.3. In Sec. IV, the dependence of optical field correlation functions on the zenith angle for the optical field propagating in turbulence is studied. These functions allow us to calculate the scintillation index, the mean beam-spot radius, and the beam wandering variance. In Sec. V, we derive the probability distribution of the transmittance that characterizes satellite-mediated atmospheric quantum channels. An application of the developed theory to the security analysis of decoy-state protocols is discussed in Sec. VI. A summary and conclusions are given in Sec. VII.

II Satellite-mediated communication links

A low Earth orbit (LEO) is an orbit around Earth with an altitude above Earth’s surface $H$ ranging from 160 to 2000 km and an orbital period between about 80 and 130 min. The fast movement of the satellite puts certain restrictions on the communication performance with the ground stations and limits the communication window to several minutes. The duration of the communication session depends on OGS’ geographical position, the duration of acquisition procedure, the presence of natural obstacles along the optical path, etc. The signal-to-noise ratio for satellite-mediated optical link is connected with the length of the optical path between the communication parties and hence it is related with the shape of the satellite orbit and its relative position to the ground-station communication party. Moreover, atmospheric refraction bends the optical ray paths, while the atmospheric turbulence corresponds to various random refraction and diffraction phenomena. In this section we remind the reader of some basic geometrical concepts of satellite-mediated optical links and leave the detailed analysis of some aforementioned aspects to the following sections.

Figure 1 shows the generic communication experiment with classical or quantum light. Two ground stations OGS1 and OGS2 establish the communication links with satellite $S$ shown in successive moments of time. The instantaneous positions of the satellite at times $t_{1}$, $t_{2}$, and $t_{3}$ are denoted as $S_{1}$, $S_{2}$, and $S_{3}$, respectively. The satellite has a polar orbit which passes over the Earth’s polar regions from south to north. Moreover, since orbits of most LEO satellites have vanishingly small ellipticity, we will consider circular orbits in the following. Initially at time $t_{1}$, the orbit, denoted as (i), passes through the meridian plane of OGS2 but has some inclination to meridian plane of OGS1. At later time $t_{2}$, the satellite approaches the horizon, makes one full orbit revolution and appears at the point $S_{3}$. For a ground observer111In the following we use the word ”observer” in general referring to the communication party located in the optical ground station. the apparent satellite trajectory changes from (i) to (ii) during this single revolution period due to the Earth rotation. In the following we will also assume that the observer’s parallax due to Earth’s rotation during one communication session can be neglected. This is justified if the flight time over the observer horizon is much smaller than the satellite orbiting period222Strictly speaking, for satellites with small altitudes above the ground and for observers located at small geographical latitudes the corrections due to the aforementioned parallax effects should be included..

The satellite zenith trajectory that passes through the observer’s meridian plane can be considered as the most favorable for establishing a communication link and is called as the ”best pass” by some authors Bourgoin2013 . Indeed, in this case the link length is smaller compared to the link’s lengths for inclined satellites. However, in many practical cases the mechanical mount of the OGS receiver or sender telescopes prevents one to use the whole advantage of the link that can lead to the interruption of the communication process near small zenith angles of the observer Casado2017 . In such cases the trade-off exists between the use of optimal propagation path and the duration of communication session. In this context the analysis of communication with the satellites moving along the inclined orbits is important.

Depending on the location of the source we distinguish downlink [paths (a)-(c) in Fig. 1] and uplink [path (d)] communication scenarios. A downlink geometry is favorable for optical signal transmission since the optical beam starts to propagate in vacuum until it enters the atmosphere. As a consequence the transmitted beam shows smaller diffraction-induced broadening and small beam wandering due to refraction on turbulent atmospheric inhomogeneities. This is the reason why the majority of quantum optical experiments are being performed in downlink configuration Guenthner2017 ; Takenaka ; Casado2017 ; Liao2018 . Under comparable atmospheric conditions an uplink communication shows inferior performance due to the influence of atmospheric turbulence already on early stages of optical signal transmission. Theoretical studies Bourgoin2013 reported the reduction in received key bits being less than one order of magnitude compared to downlink scenario. However, an uplink for quantum communication purposes poses the following advantages: simple design of satellite missions or use of already launched satellites Vallone2015 , relaxed requirements on data storage and processing equipment, variability of quantum light sources and their accessibility for maintenance and repair Pugh2017 . In this article, we will focus primarily on downlink communication scenarios.

III Regular losses of satellite-mediated link

An optical Gaussian beam undergoes diffraction-induced broadening while propagating in free space. The amount of broadening depends on propagation path length and whether beam is focused, collimated, or divergent. Since the receiver collects the incoming light with the aid of finite-aperture device such as telescope, only a fraction of the signal reaches the detector. For fixed aperture size, the signal losses due to truncation on the aperture will increase with the increase of the propagation path. Moreover, absorption and scattering on atmospheric gases and aerosols leads to the degradation of the signal-to-noise ratio as well.

In this section we consider the effects leading to the regular diffraction-induced and extinction losses in satellite-mediated links. By the word ”regular” we refer to effects that occur in a systematical manner. Among such effects we distinguish purely geometric optical path elongation due to increase of relative positions of communication parties. The regular losses in this case depend on the geographical location of the OGS, on the type of the satellite orbit, and on such characteristics as the satellite altitude, orbit inclination angle, and the satellite declination angle to the Equator plane. Additionally, regular atmospheric refraction bends the light rays, increasing the optical propagation path and contributing to loss budget. This effect is especially pronounced if the satellite is positioned close to the observer horizon. Finally, we discuss the signal loss due to atmospheric absorption and scattering.

III.1 Slant range

We consider the communication scenario with an orbiting satellite which has been described in the previous section and derive the length of the line segment connecting the ground observer and the satellite. This purely geometric length referred as the slant range does not account for any elongation effects due to atmosphere and depends on the geographical location of the observer and on the parameters of the satellite orbit.

In the following we consider a perfect polar satellite orbit that passes through both north and south poles. For simplicity, we assume that the orbit is perfectly circular with the radius $R=R_{\oplus}+H$, where $R_{\oplus}=6\,371$ km being the Earth radius and $H$ being the satellite altitude above the ground. The slant range between the OGS and the satellite is obtained then from simple trigonometric considerations, [see Fig. 2 (a) and Appendix A]

[TABLE]

where $Z$ is the zenith angle lying between the vertical direction of the observer and the direction pointing on the satellite. In Appendix A we show that for the satellite orbits inclined to the observer meridian plane on angle $\Delta\iota$ the zenith angle varies in the range $Z\in[Z_{\mathrm{min}}^{\Delta\iota},\pi/2]$ with

[TABLE]

where $\Psi$ is the geographical latitude of the observer. Consequently, the slant range (1) varies from the minimal value $L(Z_{\mathrm{min}}^{\Delta\iota})$ at zenith to the maximal value $\sqrt{(R_{\oplus}+H)^{2}-R_{\oplus}^{2}}$ at the observer’s horizon.

In order to relate the orbit inclination angle $\Delta\iota$ with the characteristics of the relative motion of the satellite and the observer, we restrict our attention to the communication scenario shown in Fig. 1 for the $\text{OGS}_{2}$. If the satellite is initially at zenith of the observer (indicated as $S_{1}$ in Fig. 1) and moves towards the observer’s horizon along the zenith trajectory (i), after one satellite orbiting period $T_{\mathrm{sat}}$ it reappears at $S_{3}$ and moves along the inclined trajectory (ii). The trajectory (ii) is then inclined to the observer’s meridian plane due to the Earth rotation. After $n$th satellite revolution the inclination angle $\Delta\iota$ between the satellite orbit plane and the observer meridian is given by

[TABLE]

where $v_{\oplus}=1669.8$ km/h is the speed of Earth’s rotation at the Equator. We refer the orbit with zero inclination, $\Delta\iota=0$, as to the zenith orbit. In this case, the orbit plane coincides with the observer meridian. In our convention, $\Delta\iota$ is positive (negative) if the satellite flies westwards (eastwards) of the observer meridian plane.

In order to get some impression of the communication geometry with the satellite with the inclined orbit, we consider two satellite positions $S$ and $S^{\prime}$ at successive times and the observer at $O$ as shown in Fig. 2 (b). These positions are chosen in such way that the slant ranges $OS$ and $OS^{\prime}$ are equal as well as the corresponding zenith angles, i.e., $\sphericalangle SOB=\sphericalangle S^{\prime}OB=Z$. Let us choose the plane $ABOC$ as the observer meridian. Then, the instantaneous position of the satellite is given by both its zenith and azimuth angles. The latter is the angle between the projected segment that connects the observer with the satellite and the reference vector in the meridian plane pointed toward the north pole. Clearly, for the inclined orbit at geographical latitudes different from the polar or equatorial ones, the positions $S$ and $S^{\prime}$ have the same zenith angle but different azimuth angles.

One can imagine that the points $S$ and $S^{\prime}$ lie on the cross sections of the satellite trajectory with the right circular cone whose height is aligned along the observer’s zenith and the apex coincides with $O$ [cf. the segment of the cone $OASS^{\prime}$ in Fig. 2 (b)]. For inclined orbits that pass above the observer’s horizon, each cone of angle $2Z$ crosses the trajectory in two points if $Z\in(Z_{\mathrm{min}}^{\Delta\iota},\pi/2)$ and in one point if $Z=Z_{\mathrm{min}}^{\Delta\iota}$. Due to this symmetry we are able to characterize the length of geometrical link between the observer and the satellite by Eq. (1) with $Z\in[Z_{\mathrm{min}}^{\Delta\iota},\pi/2]$ and ignore the detailed information on the corresponding azimuth values. The latter assumption is well justified if the detailed position of satellite is determined by a tracking system that aligns automatically the sender and the receiver telescopes. In this case, the observer coordinate system is associated with the plane $OSC$ [cf. Fig. 2 (a)] that rotates around the axis $OC$. Therefore, for optical communication with the automatical azimuth angle tracking the relevant link information is incorporated in the zenith angle dependence of the relevant quantities. It is worth to note that the cross section of the cone of given $Z$ with other satellite trajectory determines the slant ranges of the same length. For example, link $AO$ on Fig. 2 (b) for the zenith orbit $AB$ is equivalent to the links $SO$ and $S^{\prime}O$ which results in equal propagation properties of light along these links.

III.2 Regular refraction

In this subsection we calculate the elongation of the slant range due to atmospheric refraction. Atmospheric refraction phenomena are based on the fact that Earth’s atmosphere has an optical refractive index that is different from its value in vacuum. Furthermore, the value of the refractive index varies with the altitude, geographical location, and meteorological conditions and hence refraction depends on space and time variables. We refer to refraction phenomena that systematically occur in the atmosphere as to regular refraction Mahan1962 ; Wittmann1997 ; Young2006 ; Jiang2014 . Regular refraction changes with altitude in a theoretically predictable fashion. Time variation of regular atmospheric refraction has rather seasonal behavior even for large zenith angles  Sampson2003 and can be ignored.

Due to a spatial variability of the refractive index with altitude, the light coming from a distant source and reaching the ground observer propagates along a curved path rather than a straight line. As a consequence, the signal from distant objects arrives under the apparent zenith angle $Z_{a}$ rather than under the true zenith angle $Z$. Both zenith angles are related as

[TABLE]

where $n_{0}=1.00027$ is the air refractive index near the ground. We consider the effect of regular refraction on elongation of the slant range given $L(Z)$.

Earth’s atmosphere can be viewed as a spherically stratified medium with specific distribution of refractive index values within each strata. In order to obtain this distribution we use the so-called standard atmosphere model StandardAtmosphere . The standard atmosphere is an idealized steady-state representation of Earth’s atmosphere that gives values of atmospheric pressure, temperature, and other parameters for altitudes up to 1000 km. The altitude-dependent values of the air refractive index can be found using the distributions of temperature and pressure according to the Edlén equation Edlen1966 ; Birch1994 . We distinguish 10 atmospheric layers above the ground and within each layer we approximate the latitude dependence of refractive index in linear manner (for details see Appendix B) If we denote the latitude of $i$-th layer upper bound as $H_{i}$, the linear path within the layer is determined as

[TABLE]

The rest of the total optical path lies in the vacuum and is given by

[TABLE]

The angles $\Phi$ and $\Theta$ determine the relative position of incoming and refracted light rays within each layer and are complex functions of observer’s zenith angle $Z$ and the so-called refraction integral $r_{i}$ Young2006 . The explicit expressions for the lengths (5) and (6) are given in Appendix C [cf. Eqs. (116) and (118)].

It is worth to note that it was assumed that the paths segments $L_{i}$ in the form (5) are obtained with the assumption that they are linear within each atmospheric layer. This is approximately true for the standard atmosphere model, which gives the minimal ray curvature of $4.4R_{\oplus}$ near the ground at $Z_{a}{=}90^{\circ}$ (cf. Ref. Gaifillia2016 for empirical formulas). This allows us to neglect the ray curvature effects in the following. This is of course an idealization and for specific daytime and meteorological conditions the accounting for ray curvature effects turns out to be important Hirt2010 .

Figure 3 shows the path elongation factor due to atmospheric refraction,

[TABLE]

where $L(Z_{a})$ is given by Eqs. (1) and (4) and $L_{i}(Z_{a})$ are defined by Eqs. (5) and (6). Alternatively, Eq. (7) can be written as a function of the true zenith angle, i.e. $\varepsilon_{r}(Z)$. The behavior of $\varepsilon_{r}$ as a function of the apparent zenith angle $Z_{a}$ starts to deviate from $\varepsilon_{r}(Z)$ while approaching the horizon due to larger optical density along the propagation path near the horizon. The elongation factor diminishes with the growth of orbit altitude. This is due to the increase of the propagation path in the vacuum with the increase of satellite altitude. For further convenience we give a polynomial fit to the elongation factor as a function of the apparent zenith angle (in degrees) for the orbit with $H=780$ km:

[TABLE]

Finally, the slant range that accounts the elongation due to atmospheric refraction reads as

[TABLE]

where $L$ is given by Eq. (1).

III.3 Regular extinction

Another source of losses that can be considered as regular or deterministic losses are associated with molecular and aerosol scattering. For the horizontal atmospheric links, the extinction factor is $\chi_{\mathrm{ext}}=\exp[-\beta_{\mathrm{ext}}(h)L]$, where $\beta_{\mathrm{ext}}(h)$ is the extinction coefficient due to molecular absorption and scattering at given height $h$ above sea level and $L$ is the distance between the communication parties. Clearly, for elevated links this formula should be modified in order to include the variation of the extinction coefficient with the height. If the slant range to the satellite is given by Eq. (9) we can write

[TABLE]

where the observer is assumed to be located at sea level and $L^{\prime}=h\sec Z_{a}$. Since the value of the extinction coefficient depends on the number density $N(h)$ of air constituents at given height $h$ its altitude dependence can be written as

[TABLE]

where $\beta_{\mathrm{ext}}^{0}$ is the extinction coefficient at sea level and $N_{0}$ is the corresponding number density of air constituents. Here, we have used the altitude dependence of $N$ derived within the standard atmosphere model [cf. Appendix B, Eq. (97)] with the scale parameter $\mathcal{H}_{0}=6600$ m. Substituting Eq. (11) in Eq. (10) we finally derive

[TABLE]

Here we have adopted the conventional units of the extinction coefficient to be given in km*-1*. Molecular or Rayleigh scattering contributes $2.544\times 10^{-3}$ km*-1* to this coefficient for optical wavelength $\lambda=800$ nm  Elterman1968 . The aerosol distribution has a more complex dependence on altitude than the model (97). Nevertheless, for the most part of the optical propagation path one can consider that aerosol scattering contributes roughly the same amount to extinction as the Rayleigh scattering. For the total extinction coefficient in (10), we adopt therefore the value $\beta_{\mathrm{ext}}^{0}=5\times 10^{-3}$ km*-1*. We also note that a similar expression to (12) has been derived in Ref. Duntley1948 .

IV Atmospheric turbulence

In the previous section we considered the influence of regular refraction on the elongation of the optical path length for light propagating in the atmosphere. The regular refraction is caused by the variation of the air refraction index with the altitude and is predictable provided the altitude variation of temperature and pressure is known. Other types of refractive and diffraction phenomena arise due to the irregular variation in space and time of the refractive index. Such random variations are connected with temperature fluctuations and wind shear and have statistical properties of turbulent scalar fields. The strength of turbulent refractive index fluctuations also varies with the altitude. Hence, the irregular disturbances of the optical signal in satellite-mediated communication depend on the zenith angle. In this section we discuss the dependence of aperture-averaged scintillations, beam broadening, and beam wandering on the zenith angle.

IV.1 Statistical description of optical turbulence

Turbulent air motion consists of a set of vortices or eddies of various diameters, ranging from extremely large with characteristic so-called outer scale $L_{o}$ to extremely small with a scale $l_{o}$. Under the influence of inertial forces, larger eddies break up into smaller ones. This cascade process continues until the minimal scale $l_{o}$, referred to as the inner scale, is reached and dissipation of turbulent flow energy takes place. This evolution of turbulent air vortices leads to a random variability of the refractive index,

[TABLE]

where $\boldsymbol{\rho}=(x\quad y\quad z)^{T}$ and $n_{i}$ is the regular part of the refractive index within $i$th atmospheric layer [cf. Eq. (96) of Appendix B]. In the following we will omit the index $i$ for convenience and adopt the notation $n(\boldsymbol{\rho},t)$ for the refractive index (13). In Eq. (13) the $z$ axis is chosen along the optical ray while $\boldsymbol{r}{=}(x\quad y)^{T}$ lies in the transverse plane to the ray direction. At optical wavelengths, the fluctuating part $\delta n$ is of the order of $10^{-6}$ or less. When the Taylor’s ”frozen turbulence” hypothesis Taylor1938 holds, implying that the random field $n(\boldsymbol{\rho},t)$ is transported with constant velocity $\boldsymbol{v}$, therefore remaining stationary in the moving coordinate system, $n(\boldsymbol{\rho},t)=n({\boldsymbol{\rho}}-{\boldsymbol{v}}t)$, the time dependence of refractive index is incorporated in the spatial variable.

In the statistical theory of optical turbulence, the refractive index fluctuations $\delta n$ are described by the power spectrum

[TABLE]

where the $B_{n}$ is the correlation function

[TABLE]

For a locally isotropic random field the correlation function simplifies further and the turbulent spectrum can be written as

[TABLE]

where $\boldsymbol{\kappa}{=}(k_{x}\quad k_{y})^{T}$ and $\boldsymbol{r}{=}\boldsymbol{r}_{1}{-}\boldsymbol{r}_{2}$. For the inertial range of spatial wave-numbers values, $|\boldsymbol{\kappa}|\in[2\pi/L_{o},2\pi/l_{o}]$, the spectrum for Kolmogorov turbulence is

[TABLE]

Here, $C_{n}^{2}$ is the refractive index structure parameter characterizing the strength of refractive index fluctuations.

In the following considerations we will require the specific form of correlation function, called structure function, for optical phase fluctuations:

[TABLE]

where $\boldsymbol{r}_{i}$, $\boldsymbol{r}_{i}^{\prime}$, $i=1,2$, are the transverse components of the corresponding spatial vectors $\boldsymbol{\rho}_{i}=(\boldsymbol{r}_{i}\quad z)$ and $\boldsymbol{\rho}_{i}^{\prime}=(\boldsymbol{r}_{i}^{\prime}\quad z)$. Here, the fluctuating phase is found from the first approximation of geometric optics

[TABLE]

for uplink and

[TABLE]

for downlink communication scenarios. In this article we focus our attention on the downlink communication configuration and omit the corresponding superscripts for simplicity in the notations. The formulas for uplink can be then obtained by replacing $z^{\prime}\rightarrow L-z^{\prime}$ in integrals similar to the one in Eq. (21).

For a locally isotropic and homogeneous random field $\delta n$, the Markov approximation,

[TABLE]

is well justified Tatarskii2016 ; Baskov . The phase structure function (19) is evaluated then as

[TABLE]

For the Kolmogorov spectrum (18), in the case of downlink the structure function (IV.1) reduces to Gurvich

[TABLE]

where

[TABLE]

is the radius of spatial coherence of a plane wave in the atmosphere, $L_{\mathrm{turb}}$ is the propagation length within the optically active turbulent atmospheric layer, $L_{r}$ is the total propagation length, $\xi$ is the dimensionless integration variable, and

[TABLE]

is the refractive index structure function taken at some reference height $h_{0}$ above the ground (see also Appendix D). The structure function (24) for uplink is obtained by the change of the variable $\xi\rightarrow 1-\xi$.

IV.2 Aperture-averaged scintillations

The dependence of optical intensity fluctuations, i.e., scintillations, on the zenith angle has attracted attention in the connection with astronomical photometry. Early measurements of the scintillation index $\sigma_{\eta}^{2}$ Butler1952 ; Ellison1952 ; Briggs1963 have shown that the aperture-averaged scintillations (or power scintillations) grow with the growth of the zenith angle as

[TABLE]

where

[TABLE]

is the transmittance of the light intensity $I(\boldsymbol{r},L_{r})$ through the receiver aperture with opening area $\mathcal{A}$, and $\boldsymbol{r}$ is the spatial variable transversal to propagation direction. The exponent $\gamma$ in Eq. (27) is related to the statistics of turbulent fluctuations of the refractive index and depends on the characteristics of the receiver telescope. Theoretical considerations Tatarskii2016 based on the Rytov approximation yield $\gamma=11/6$ for small receiving apertures and $\gamma=3$ for large receiving apertures, respectively. These results agree reasonably with the experiments Stecklum1985 .

The later investigations Burke1970 ; Stecklum1985 ; Kucherov ; Fuentes1987 have shown that dependence (27) is valid for small and moderate zenith angles or for highly elevated optical ground stations. For large zenith angles measured stellar scintillations exhibit a saturation or decrease in its value. A similar behavior has been reported for optical signals from satellites Moll2014 ; Moll2015 . On the other hand, the Rytov approximation (27) yields the divergent scintillation index for $Z_{a}{\rightarrow}90^{\circ}$ and does not account for the saturation effect. The phenomenon of saturation of scintillations has been theoretically studied in Refs. Young1969 ; Young1970 ; Young1970a ; Beran1988 ; Andrews and is attributed to multiple scattering phenomena on turbulent inhomogeneities.

In this section we derive the scintillation index of the aperture-averaged optical signal from the satellite and discuss its dependence on the zenith angle. As an accompanying result we obtain the expressions for two first moments of the transmittance (28), the mean width of the beam spot at the receiver and the beam wandering variance. These parameters will be used in Sec. V for the derivation of the probability distribution of quantum channel transmittance.

The scintillation index (27) is derived from the first two moments of the transmittance (28), namely

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

$a$ is the receiving aperture radius, and $u(\boldsymbol{r};L)$ is the optical field amplitude. Obviously, the field correlation functions $\Gamma_{2}$ and $\Gamma_{4}$ turn to be important ingredients for the evaluation of $\sigma_{\eta}^{2}$. On the other hand, their moments are related to the beam wandering variance,

[TABLE]

and to the mean short-term beam spot radius of the transmitted beam,

[TABLE]

Here the variable $x$ denotes the $x$ component of the transverse vector $\boldsymbol{r}$. The short-term beam spot radius is associated with the intensity distribution observed during small exposure times while the long-term radius $W_{\mathrm{LT}}$ includes broadening effects due to beam wandering and is associated with long detection times. In the following section we show that the quantities (29), (30), (33), (34) are of primary importance for the description of quantum atmospheric channels.

The correlation functions $\Gamma_{2}$ and $\Gamma_{4}$ we calculate using the phase approximation of the Huygens-Kirchhoff method Banakh1977 ; Banakh1979 . This approximative method neglects the fluctuations of the field amplitude. On the other hand, the phase fluctuations characterized by Eqs. (20) and (21) are caused by random diffraction and refraction in the turbulent medium, which are included up to the terms of order $\langle\delta n^{2}\rangle^{1/2}$. The former condition is justified for communication scenarios through links with saturated turbulence  Holmes1980 , which is the case for satellite-mediated atmospheric channels. More rigorous methods Klyatskin1974 and experiments Tatarskii2016 ; Artemev1971 suggest that the latter condition is satisfied for arbitrary propagation paths and turbulence strengths. Thus, the phase approximation of the Huygens-Kirchhoff method accounts for the diffraction and refraction effects arising due to the optical beam propagation in the turbulent atmosphere in most practically important cases.

For the TE00 mode of the laser beam in the plane of the radiating aperture,

[TABLE]

with $W_{0}$ being the initial beam spot size and $F$ being the initial beam wavefront radius, the correlation functions (31) and (IV.2) read as

[TABLE]

[TABLE]

[TABLE]

where we have used the relative and the center-of-mass coordinates $\boldsymbol{r}=\boldsymbol{r}_{1}-\boldsymbol{r}_{2}$ and $\boldsymbol{R}=(\boldsymbol{r}_{1}+\boldsymbol{r}_{2})/2$, respectively. Here, $\Omega=kW_{0}^{2}/(2L_{r})$ is the Fresnel number of the transmitter aperture, $g^{2}=1+\Omega^{2}[1-L_{r}/F]^{2}$ is the generalized diffraction parameter, and the phase structure function $\mathcal{D}_{S}$ is given by (24). We also note that Eqs. (36) and (IV.2) incorporate multiple-scattering effects.

IV.2.1 Mean transmittance

The average intensity at the receiver aperture is found from Eqs. (31) and (36). We use the quadratic approximation Andrews2005

[TABLE]

for the spatial dependence of the phase structure function (24). This approximation gives a good accuracy for small values of the radius of spatial coherence $\rho_{0}$, i.e., for long propagation distances and strong optical turbulence, which is always the case for satellite-mediated links. Performing the integration in Eq.  (36), we obtain the Gaussian distribution of the intensity distribution at the receiver aperture plane,

[TABLE]

where the long-term beam spot radius at the receiver is [cf. also Eq. (34)]

[TABLE]

Here, $\rho_{0}$ is defined in Eq. (25) and

[TABLE]

is the weighting factor that depends on the slant profile of the refractive index structure function $C_{n}^{2}(L_{r},1-\xi)$. This profile is related with the vertical profile of turbulence as $C_{n}^{2}(L_{r},1{-}\xi)=C_{n}^{2}[h(\xi){=}L_{r}(1-\xi)\cos Z_{a}]$ (see Appendix D for details). The natural diffraction of the laser beam is included in Eq. (41) by considering the initial beam waist $W_{0}$, wavelength $\lambda$ (Fresnel number $\Omega$), and propagation distance $L_{r}$. The resulting natural laser beam divergence is usually the main source of signal loss. The obtained formula (41) coincides with the long-term beam spot radius given in Ref. Fante1975 .

The obtained mean intensity allows one to calculate the mean transmittance straightforwardly. Inserting Eq. (40) in Eq. (28) and performing the integration with respect to the spatial variable, we obtain

[TABLE]

This is the transmittance of Gaussian beam with the beam-spot radius $W_{\mathrm{LT}}$ through a circular aperture of radius $a$. The approximation sign is used due to the used quadratic approximation (39). This formula can serve for estimative calculations only and for the precise evaluations from Eqs. (29) and (36) we derive

[TABLE]

where $J_{n}(x)$ is the Bessel function of $n$th order.

For the sake of completeness, we give the expression of the short-time beam spot size (34). Following Ref. Yura1973 it is calculated as

[TABLE]

Here, the term proportional to $\mathcal{X}$ arises due to the diffraction-induced beam broadening caused by the turbulent atmosphere. In the limit $\mathcal{X}\rightarrow 0$ we obtain from Eq. (45) the diffraction-induced beam broadening in vacuum. We also note that Eq. (45) does not account for the beam broadening due to the random scattering on particles, aerosols, and precipitations. This additional broadening might be relevant under the condition of high moisture or low visibility, as has been shown in Ref. Vasylyev2017 . In this article, however, we deal with ideal weather conditions for optical communication. The short-term beam spot radius (45) is important for the derivation of the probability distribution of the channel transmittance as will be shown in the next section.

IV.2.2 Second moment of transmittance

The second moment of transmittance is obtained by substituting Eq. (IV.2) in (30). Unfortunately, the integration cannot be performed in analytic form. However, for satellite-mediated atmospheric links we can simplify the integral kernel (38) for the further numerical integration (for details see Appendix E):

[TABLE]

This expression can be simplified further by noting that most of the optical propagation path lies in vacuum as well as in the atmosphere with negligibly small turbulence, i.e., $L_{\mathrm{turb}}\ll L_{r}$. As a consequence, for small values of the integration variable $\xi$ in (IV.2.2), when the contributions from the $\boldsymbol{r}(1-\xi)$ term are dominant, the value of the refractive index structure constant is equal to zero. On the other hand, for $\xi\approx 1$ the structure constant is finite but the contribution from $\boldsymbol{r}(1-\xi)$ is negligibly small. Therefore, one can neglect the dependence on $\boldsymbol{r}$ in (IV.2.2) along the whole propagation path

[TABLE]

We stress that this formula is justified for downlink configurations only. Inserting (IV.2) in (30) and using the approximation (47), we obtain

[TABLE]

The further evaluation of integrals in (IV.2.2) should be performed numerically.

Figure 4 shows the aperture-averaged scintillation index $\sigma_{\eta}^{2}=\langle\Delta\eta^{2}\rangle/\langle\eta\rangle^{2}$ as a function of the zenith angle calculated by using the phase approximation of the Huygens-Kirchhoff method. For the calculation of $\langle\eta^{2}\rangle$ we have used the approximate expression (IV.2.2), whereas $\langle\eta\rangle$ is calculated from Eq. (44). Table 1 lists the atmospheric and optical beam parameters used in the calculation of the scintillation index. The curves shown for three inclination angles yield the same functional dependence on zenith angle for medium and large $Z_{a}$ due to the properties discussed in Sec. III.1 and differ on the minimal value of the zenith angle $Z_{\mathrm{min}}^{\Delta\iota}$ [ see Eq. (2)]. For small and moderate zenith angles the scintillation index $\sigma_{\eta}^{2}$ shows the asymptotic behavior given by Eq. (27). The scintillation index calculated within the phase approximation of Huygens-Kirchhoff method shows the saturation and decrease of intensity fluctuations for large zenith angles.

This result can be qualitatively compared with the experimental data that have been taken within five measurement campaigns between 2006 and 2016 and shown in Fig. 5 (shaded area). The direct quantitative comparison of theoretical curves in Fig. 4 with any individual experimental curve contributing to Fig. 5 is hardly possible due to the lack of all needed parameters of turbulence and the hardly measurable profile of $C_{n}^{2}(z)$ inherent to the given meteorological conditions. Measurements of intensity scintillation index were conducted with an optical ground station nearby Munich, Germany, with elevation of about 602 m above sea level. Measurement wavelength was 847 and 1550 nm, depending on the used satellite (see Table 2 for an overview of measurement campaigns).

We also note that the difference in zenith angles with maximal scintillation index in Figs. 4 and 5 arises due to the different elevations of the observers.

The measurement device that delivered the data for the scintillation index analysis is an infrared camera located in the exit pupil of the telescope. Thus, images of the intensity field incident onto the telescope aperture are recorded and analyzed. The experimental curves in Fig. 5 are obtained based on data from a single camera pixel. The effective radius of the pixel, taking magnification of the optical system into account, is $a=3.2$ mm. The small value of the detector aperture as well as short exposure times (of order 0.1-1 ms) reduce the telescope aperture smoothing effect, yielding the intensity scintillation index

[TABLE]

For the sake of better comparison, the 1550-nm data are re-calculated to 847-nm wavelength using weak scattering theory. This is possible since the 1550 nm measurements lie well within the weak scattering regime, i.e. at low zenith angles. Twelve measurements are analyzed to form the mean run of scintillation index as shown in Fig. 5. The grey area indicates the confidence bound defined by the standard deviation. At low and high zenith angles, the standard deviation is not illustrated since only a single measurement track is recorded in these regimes and, thus, determination of standard deviation is not possible. Further description of the individual measurement campaigns and the data analysis method are found in Moll2015 ; Moll2016 .

In Appendix F we consider an estimation for the aperture-averaged scintillation index based on approximate phenomenological expressions for field correlation functions. This can give some insight into saturation effects observed at large zenith angles. The consideration yields

[TABLE]

Here, $a$ is the receiving aperture radius, $k=2\pi/\lambda$ is the optical wave number, and ${}_{2}F_{3}(a,b;c,d,e;x)$ is the hypergeometric function. This simple analytic formula contains three phenomenological parameters: the refractive index structure parameter at the ground $C_{n,0}^{2}$, the characteristic height of the atmospheric turbulent layer $H_{0}$, and the dimensionless proportionality parameter $\mu$. The theoretical curve based on Eq. (49) shows a reasonable agreement with the experimental data in Fig. 5. The model parameters are: $C_{n,0}^{2}=2.5{\times}10^{-17}$ m*-2/3*, $H_{0}=0.5$ km, $a=3.2$ mm, $\lambda=847$ nm, and $\mu=0.92$.

IV.2.3 Beam wandering variance

For the sake of completeness, we give the expression for the beam wandering variance (33). Beam wandering phenomenon depends strongly on the outer scale of turbulence, $L_{o}$. The finite $L_{o}$ defines the upper bound on the size of turbulent inhomogeneities that are able to deflect the beam as a whole. Since the Kolmogorov spectrum (18) has a discontinuity at the turbulent wave numbers, $|\boldsymbol{\kappa}|\approx\kappa_{o}=2\pi/L_{o}$, it is more suitable to use the smoothed spectrum

[TABLE]

For slant paths the outer turbulence scale varies with the height $h$. Very well known is the empirical Coulman-Vernin profile that reads as Coulman1988

[TABLE]

where the outer scale is given in meters.

Using the spectrum (50) and following the derivation steps of Ref. Mironov1977 we obtain

[TABLE]

Here $W_{\mathrm{ST}}(L_{r})$ is given by (34) and the altitude dependence of the outer scale is given by (51) with $h\approx L_{r}\cos Z_{a}$. Appendix D summarizes the model for the $C_{n}^{2}$ profile needed for the calculation of the beam wandering variance (IV.2.3).

Figure 6 shows the divergence angle $W_{\mathrm{ST}}/L_{r}$ and the standard deviation of the angle-of-arrival fluctuations $\sigma_{\mathrm{BW}}/L_{r}$ as a function of the apparent zenith angle and several inclination angles $\Delta\iota$ of the satellite orbit. For the simulation of turbulent atmosphere we use the AFGL+WK (the Air Force Geophysics Laboratory and Walters-Kunkel model) night model for the altitude variation of the refractive index structure parameter (see Appendix D). Other relevant atmospheric parameters are listed in Table 1. The beam spot radius and the beam wandering variance grow with the growth of the zenith angle. The asymmetry of the curves for inclined orbits is dictated by the geographical position of the observer as discussed in Sec. III.1. For zenith angles near $Z_{a}=80^{\circ}$ some saturation of both quantities appears. This is due to the drop of the turbulence strength in the lower troposphere (3-10 km), the region that has its maximal contribution to the optical path at this zenith angle. Near the observer horizon the contributions of boundary layer turbulence at 1 km as well as of the atmospheric refraction have a maximal effect on the optical beam distortions leading to the growth of both beam-spot size and beam wandering.

V Probability distribution of transmittance

In this section we consider the model of atmospheric quantum channels. For the physical consistence of this model the preservation of the canonical commutation relations for the quantized optical field operators is important. This requirement puts certain restrictions on the probability distribution that governs the statistics of fluctuating channel transmittance. We also show the relationship of the statistical characteristics of this distribution and the moments of the field correlation functions derived in the previous section.

Quantum light transmission through a linear medium, such as Earth’s atmosphere, is conveniently characterized via the input-output relations

[TABLE]

where $\hat{a}_{\mathrm{in(out)}}$ is the input (output) field annihilation operator and $\hat{c}$ is an environmental mode operator. The random transmittance $\eta\in[0,1]$ equals to the instantaneously transmitted normalized intensity truncated by the receiver aperture [cf. Eq. (28)]. The transmission of the quantum state through the atmospheric link depends not only on the characteristics of the channel, but also on the parameters of the receiver aperture. In terms of the Glauber-Sudarshan $P$ function Glauber1963 ; Sudarshan1963 , the input-output relation (53) reads Semenov2009 as

[TABLE]

where we have assumed that the environmental modes are in the vacuum state. Here, $P_{\mathrm{in}}(\alpha)$ and $P_{\mathrm{out}}(\alpha)$ are $P$ functions of the input and output quantum fields and $\mathcal{P}(\eta)$ is the probability distribution of the transmittance (PDT). Hence, the description of quantum-light propagation through the atmosphere reduces merely to identifying $\mathcal{P}(\eta)$.

An important requirement for the PDT is that its domain of definition is restricted to the interval $[0,1]$, which is the consequence of the canonical commutation relation for $\hat{a}_{\mathrm{out}}$ and $\hat{a}_{\mathrm{out}}^{\dagger}$. Violation of this requirement may lead to unphysical effects which have critical influence, e.g., on security bounds of communication protocols. The first consistent model of the PDT Vasylyev2012 considered beam wandering as the main source of fluctuating losses on the receiver. In this case, the PDT can be derived in analytical form and is given by the log-negative Weibull distribution. This beam wandering model was further extended to include effects due to random beam broadening and deformation of the beam profile into an elliptic form Vasylyev2016 . In Ref. Vasylyev2017 this so-called elliptic-beam model was extended to include additional beam broadening and extinction due to random scattering on atmospheric aerosols and dust particles. The elliptic-beam model gives reasonable agreement with the PDT measurements for short-distance links. For long-distance channels, a discrepancy arises between the transmittance moments calculated via the elliptic-beam PDT and the corresponding moments (29) and (30) calculated from the first principles. For elimination of this discrepancy, a PDT model has been proposed based on the law of total probability Vasylyev2018 . This model is most suitable for the description of long-length quantum channels and in this section we extend it to the case of slant propagation paths.

Analyzing experimental PDTs for short and long propagation distances, one observes two limiting behaviors. For short distances, the beam wandering appears to be the major source of fluctuating losses and the corresponding PDT has similar form to the log-negative Weibull distribution. For long propagation distances, beam broadening and deformation are the dominating effects yielding the PDT in the form of truncated log-normal distribution. In the general case, using the law of total probability Schervish1995 , the PDT can be written as

[TABLE]

where the random vector $\boldsymbol{r}_{0}$ transverse to the propagation direction describes the position of the deflected beam centroid relative to the aperture center. The corresponding probability distribution in (55),

[TABLE]

describes the beam-wandering contribution to the total PDT. Here, $\sigma_{\mathrm{BW}}^{2}$ is the beam wandering variance, Eq. (33). The effects of beam-spot distortions are incorporated in the conditional distribution $P(\eta|\boldsymbol{r}_{0})$. Physically it can be interpreted as the conditional PDT for the beam with a centroid position tracked to the position $\boldsymbol{r}_{0}$ relative to the aperture center.

For negligible small beam wandering, $\sigma_{\mathrm{BW}}\rightarrow 0$, the distribution (56) reduces to Dirac delta function and the conditional probability according to (55) reduces to the PDT. In this limit, the PDT resembles the log-normal distribution. Hence, for general situations we can assume that the conditional distribution can be approximated by the truncated log-normal distribution

[TABLE]

where $\mathcal{F}(1)$ is the cumulative function of the (untruncated) log-normal distribution at $\eta=1$. The parameters of this distribution are related to the conditional moments $\langle\eta\rangle_{r_{0}}$ and $\langle\eta^{2}\rangle_{r_{0}}$. These conditional moments can be considered as the corresponding moments of the aperture transmittance of the effective beam with beam-spot radius $W_{\mathrm{ST}}$ whose centroid is displaced to the distance $r_{0}$ from the aperture center. In the limit of weak beam wandering, the conditional moments can be written as (see Ref. Vasylyev2018 for details)

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here, $W_{\mathrm{ST}}$ is given by (45) and ${\rm{I}}_{n}(x)$ is the modified Bessel function of $n$th order. The parameters of the conditional distribution (57) are determined from Eqs. (61) and  (62) approximately:

[TABLE]

[TABLE]

The knowledge of the parameter set $\{\langle\eta\rangle,$ $\langle\eta^{2}\rangle,$ $\sigma_{\mathrm{BW}},$ $W_{\mathrm{ST}}\}$ is therefore sufficient for the determination of the channel PDT. For detection with a Cassegrain-type telescope, the PDT derivation procedure can be further generalized as described in Ref. Vasylyev2018a .

Practical optical communication via satellites is impossible without acquiring and tracking the received signal. Both beam wandering due to atmospheric turbulence Andrews2005 or satellite vibrations Arnon1998 cause changes in the direction of the received beam that result in misalignment between the communication parties. Moreover, the velocity aberration point ahead and the atmospheric dispersion effects Oi2017 should be taken into account. The compensation of these disturbances requires an active beam steering that is accomplished by mechanical means. Quantum key exchange with low mean intensities requires especially precise beam tracking and stable pointing Moll2012 ; Bourgoin2015 . Applying coarse and fine tracking and pointing strategies, one can achieve a tracking accuracy of approximately $\theta_{\mathrm{tr}}\sim 1.2\mu$rad for LEO-ground communication links Liao2017a . We incorporate the tracking procedure in the PDT model (55) by replacing

[TABLE]

in Eq. (56). Finally, we note that the parameters (63) and (66) still contain the beam wandering variance, and we calculate $\sigma_{\mathrm{BW}}^{2}$ by means of Eq. (IV.2.3).

VI Application: Decoy state protocol

We apply the developed theory of satellite-mediated quantum atmospheric channels for the calculation of the secret key rate when decoy states are in use Hwang2003 ; Lo2005 . Conventionally, we refer to communication parties Alice and Bob as to the sender and the receiver, respectively. Based on the BB84 protocol Bennett1984 , the decoy-state method estimates channel parameters by sending two types of states. While one type of states (signal states) is used for transmission of quantum keys, the other is called the decoy state, which is used for the estimation of the number of transmitted single-photon states. Ideally, the single-photon states are most suitable to be used as the signal. Practically perfect single-photon sources are hard to attain and one uses weak coherent states instead. In the security analysis of the decoy-state method, both signal and decoy states possess equal properties except for their intensity.

Usually, only a few decoy states are needed for practical implementations. A simple two-decoy-state protocol with vacuum+weak decoy states gives an optimal key generation rate which is the same as having an infinite number of decoy states Ma2005 . On the first stage of the protocol, Alice, who has a phase-randomized source of coherent states, encodes the bits in the $X$ or $Z$ basis as in the standard BB84 scheme, e.g., by utilizing polarization degrees of freedom. Additionally to the signal field, she generates decoy states in vacuum and weak coherent states. The phase randomization makes the source statistically equivalent to a Poissonian distribution of Fock states such that, when the average photon number from the light source is $\mu$, the probability to send an $n$-photon pulse is $e^{-\mu}\mu^{n}/n!$. We denote the mean photon numbers as $\mu_{j}$, $j=s,d,v$ for signal, weak-decoy, and vacuum states. The following conditions are satisfied: $\mu_{d}<\mu_{s}<1$, $\mu_{v}=0$. After transmission through the free-space channel, Bob performs measurement of transmitted bits in a randomly chosen $X$ or $Z$ basis. The conditional probability of a detection event at Bob’s side given that Alice sends an $i$-input state is referred to as the yield $Y_{i}$ of an $i$-photon state. The vacuum state is used for the estimation of the background detection probability $Y_{0}$ while the weak-decoy state allows one to estimate the single-photon yield $Y_{1}$ and the error rate of the single-photon state $e_{1}$. Since both signal and weak-decoy states propagate through the same channel, the single-photon yield would be the same for these states. Due to this property, the security of the decoy-state protocol against the photon-number splitting attacks Scarani2009 ; Pfister2016 or the Trojan-horse attacks Tamaki2016 can be verified.

As the next step, the parties perform the sifting of the raw key, its error correction, and privacy amplification. Finally, if these steps were successful, Alice and Bob share a shorter but more secure key. Defining the averaging over channel fluctuations with the applied beam tracking procedure (69) as

[TABLE]

we find for the lower bound of the average secure key rate

[TABLE]

Here, $q$ depends on the implemented protocol, $H(x)=-x\log_{2}(x)-(1-x)\log_{2}(1-x)$ is the binary Shannon information function, $Q_{\mu_{s}}$ is the gain of the signal states, $Q_{1}$ is the gain of single-photon states, and $f(x)\geq 1$ is the bidirectional error correction efficiency ($f(x)=1$ corresponds to the perfect error correction case). The quantum bit error rate, QBER, is estimated as

[TABLE]

The gain $Q_{\mu_{s}}$ represents the ratio between the number of events where Bob observes a click under the condition that Alice sends a certain number of signal states. This overall gain with respect to an ideal threshold detector Kok2007 can be evaluated as Ma2005

[TABLE]

where the yield of the $i$-photon state is

[TABLE]

Here

[TABLE]

accounts for deterministic losses such as detector efficiency, channel losses due to atmospheric extinction (10), and absorption by optical components, while $\eta$ is the random transmittance of the free-space channel [cf. Eq. (28)]. Similarly, defining the error rate of the $i$-state as

[TABLE]

one derives for the overall error gain

[TABLE]

Here, $e_{\mathrm{det}}$ is the probability that an incorrect bit value occurred that depends on the alignment and the stability of the optical system. The background error rate is $e_{0}=1/2$ for randomly occurring dark and background counts. The dark count contribution to $Y_{0}$ is of order $10^{-6}$ for a commercially available Geiger-mode APD at room temperature and can be further decreased with proper cooling. The contribution from transmitted vacuum decoy states with the accounting of finite-size effects is given in Appendix G. The total value $Y_{0}$ can be further enhanced due to sky-noise photodetection Erlong2005 ; Gruneisen2015 ; Nordholt2002 . In this article we assume that the background error rate is constant with the value taken from Ref. Liao2017 .

For the determination of the single-photon gain, $Q_{1}$, and the single-photon error rate, $e_{1}$, in Eq. (71) the statistical fluctuations must be considered. Indeed, since the communication link with a LEO satellite can be established for only several minutes, only a finite set of data can be transmitted. In the security analysis the accounting for possible deviations from most probable values must be taken into account. Moreover, statistical fluctuations tend to become more important as the distance of the QKD increases, i.e. for large values of the zenith angle. In Ref. Ma2005 the Gaussian model and in Refs. Curty2014 ; Lim2014 the Chernoff-Hoeffding method have been applied for deriving finite-key security bounds. We adopt the statistical fluctuation analysis of Ref. Zhang2017 that uses the Chernoff bound for establishing the lower bound for single-photon gain,

[TABLE]

and the relation between the upper bounds for single-photon error rates,

[TABLE]

in $X$ and $Z$ bases. Appendix G summarizes the method for obtaining $\theta^{U}$ [cf. Eq. (169)], the lower bound for the single-photon yield $Y_{1}^{\gamma L}$ [cf. Eq (161)] and the upper bound for the bit-flip error rate $e_{1}^{xU}$ [cf. Eq. (162)].

We calculate the secure key rate for the downlink communication scenario. Table 3 lists the values of parameters used in the calculation, whereas Table 1 gives parameters of the atmospheric link and communication system. With the source of repetition rate 150 MHz, the total number of generated bits by Alice during 14 minutes of communication session is $N=10^{11}$. We estimate the mean number of sifted key bits,

[TABLE]

where $\eta_{\mathrm{sift}}=0.5$ is the sifting efficiency for conventional sifting protocols. The number of sifted key bits is averaged over the fluctuations of the channel transmittance and depends on the zenith angle (e.g., for zenith orbit $\langle M^{s}\rangle_{\mathrm{tr}}=1.46{\times}10^{5}$, $\langle M^{d}\rangle_{\mathrm{tr}}=5.64{\times}10^{4}$ for $Z_{a}=0^{\circ}$, and $\langle M^{s}\rangle_{\mathrm{tr}}=1.93{\times}10^{4}$, $\langle M^{d}\rangle_{\mathrm{tr}}=7.43{\times}10^{3}$ for $Z_{a}=70^{\circ}$).

Figure 7 shows the QBER and the lower bound of the secure key rate calculated for several inclination angles of satellite orbits relative to the observer’s meridian plane. For small zenith angles scintillations, atmospheric refraction, and absorption have minimal impact on the performance of decoy-state protocols. The QBER has almost constant value, the phenomenon that was theoretically predicted for satellite links in Ref. Shapiro2011 and experimentally observed in Refs. Casado2017 ; Liao2017 ; Liao2017a . Figure 7 a) shows a similar behavior for the zenith angles where the aperture-averaged scintillation index has the asymptotic behavior (27) and satisfies the condition $\sigma^{2}_{\eta}<1$. In the region of saturated scintillations the QBER grows making quantum communication impossible for large zenith angles. In this region the wavefront distortions are maximal. These distortions together with the extinction losses lower the signal-to-noise ratio and correspondingly the QBER is inflated by the larger relative contributions from the background. The flat region of the QBER diminishes with the growing inclination angle of the satellite orbit allowing smaller time windows for the secure key exchange as shown in Fig. 7 b).

It is worth to note that by placing the observer at high altitudes above sea level this time window can be increased. Indeed, the thickness of the dense ground layer, which is responsible for scintillation, is different at sea level and at high altitudes. This causes a shift of the region with saturated scintillations to higher zenith angles. This phenomenon can be observed if one compares the curves of Fig. 4 calculated for the observer at sea level with the experimental curve of Fig. 5 with observer’s elevation of 602 m above sea level. The saturation region may even vanish at particularly high altitudes, the phenomenon which is known in optical astronomy Kucherov ; Fuentes1987 ; Butler1952 . This observation makes observatories in mountains especially attractive as a OGS node for quantum free-space communication with satellites.

VII Conclusions

In this article we have presented the theoretical analysis of satellite-mediated quantum links. We have discussed the influence of regular refraction, extinction, and turbulence on the transmission properties of optical signals through the Earth’s atmosphere. Since in satellite-mediated communication scenario the position of the satellite changes rapidly for the observer located on the Earth’s surface, the thoughtful analysis is presented of how atmospheric disturbances depend on the observer’s geographical position, observer’s zenith angle, and the orbit inclination angle.

We focus our analysis on low orbit satellites with perfect polar orbit. In particular, we considered the case when the satellite orbit is inclined to the observer meridian and derived the corresponding slant range. The orbit inclination restricts the definition range of the zenith angle by introducing the lower bound of the angle. As a consequence, the most favorable satellite trajectory corresponds to the zenith orbit that passes through the observer meridian. In this case, the smallest zenith angle is zero and the slant range is the shortest one at the observer zenith.

The effect of regular atmospheric refraction increases the optical slant range and changes the value of the true zenith angle to the apparent one. This effect is especially pronounced near the horizon where it can increase the propagation path of the optical signal up to $30\%$. Based on the standard atmosphere model we derived the corresponding path elongation factor as a function of the apparent zenith angle and gave the corresponding analytical fit formula. The resulting slant range has been used for the calculation of the atmospheric extinction factor due to absorption and scattering on atmospheric gases and aerosols.

Another important factor that deteriorates the optical performance of satellite-ground links is the atmospheric turbulence. It appears that the second- and the fourth-order optical field correlation functions play a central role in the description of light propagation through the turbulent media. These functions allow one to derive the aperture-averaged scintillation index, as well as the mean beam-spot radius and the beam wandering variance of the transmitted beam. Based on the properties of intensity covariance function, we have derived an analytical expression for the aperture-averaged scintillation index. For large zenith angles these scintillations saturate and decrease; this effect is well observed in experiment. It arises due to multiple scattering in most turbulent air layers near the ground, the process that degrades the performance of the receiver telescope and leads to the additional aperture-averaging of scintillations. Based on this simple model of aperture averaged scintillations, we have developed the rigorous approach for calculating field correlation functions and their moments.

The developed description of atmospheric channels has been adopted for the description of quantum light propagation through the Earth’s atmosphere. For this sake we use the input-output relations for optical field operators, rewritten in terms of Glauber-Sudarshan $P$ function. Fluctuations of the channel transmittance due to atmospheric turbulence are accounted with the help of the probability distribution of transmittance. We have obtained the latter for satellite-mediated quantum links using the law of total probability.

Finally, the security of quantum decoy-state communication protocols is analyzed with realistic channel parameters and communication conditions. In this connection, the inclusion of the finite-key effects plays an important role. Scintillation phenomena at large zenith angles influence greatly the performance of quantum channels, leading to a growth of the quantum bit error rate. In the region of saturated scintillations that appears close to the horizon, no secure quantum key can be obtained. For small zenith angles the quantum bit error rate has minor variation that allows one to obtain secure key bits within a certain time window. This window increases with the decrease of the satellite inclination angle relative to observer’s meridian plane. We also expect that the increase of the observer’s altitude relative to the sea level will lead to a decrease of the region with saturated scintillations and hence to a better communication performance.

We have omitted several aspects affecting the performance of the satellite-mediated quantum communication. Background radiation from the Sun, Moon, stars, or light reflected from the satellite introduce additional noise. This noise is detected by photodetectors, and therefore, it increases the quantum bit error rate. The use of light buffers, time gate, frequency, and spatial filters can partially mitigate the problem of background noise. In this study we have neglected the parallax-connected errors due to the Earth’s rotation during one communication session. For observers at small geographical latitudes, and for the satellites with small altitude such parallax effects should be included in the rigorous analysis. On the other hand, for the analysis of realistic communication scenarios not only the instantaneous position of a satellite relative to the observer is important, but also the relative position of the Sun or of the Moon. We refer to extensive literature that study the influence of background noise and dark counts on the satellite-mediated quantum key distribution.

Much less studied is the influence of random scattering on atmospheric aerosols and dust particles on the performance of satellite-based communication. Especially in overcast conditions the dominance of Mie scattering makes the optical quantum communication with satellites impossible. We have included in the model atmospheric extinction effects due to absorption and scattering on molecules and aerosols which correspond to average communication conditions. Under the conditions of low visibility and overcast skies, the performed analysis is not applicable.

Finally, here we have considered fundamental spatial modes for the light beams that have Gaussian profiles of the intensity distributions. The efficient generation of diffraction-free beams will allow one to reduce losses associated with regular and random diffraction in free-space and may extend the range of satellite-mediated quantum communication.

Acknowledgements.

The work was supported by the Deutsche Forschungsgemeinschaft through projects VO 501/21-2 and VO 501/22-2. The authors are grateful to A. A. Semenov for useful and enlightening discussions.

Appendix A Slant range of optical beam

In this appendix we derive the length of geometrical distance between the satellite and the optical ground station.

A.1 Geometry of satellite communication links

The geometric path length for a specific communication scenario is determined from the value of the zenith angle $Z$, which in turn is a function of the satellite orbit inclination, the instant position of the satellite, and geometric latitude of the observer. In the Earth surface coordinate system of the observer the position of the satellite is given by the zenith angle $Z$ and the azimuth angle $A$. The latter we measure from the north point eastwards. In the Earth center (geocentric) coordinate system the satellite position is determined by its declination angle and its orbit inclination angle. For circular orbits the declination angle is given by

[TABLE]

where $\omega_{\mathrm{sat}}=2\pi/T_{\mathrm{sat}}$ is the angular speed of the satellite with orbiting period $T_{\mathrm{sat}}$ and $t_{\mathrm{Pole}}$ is the reference time associated with the trajectory crossing either the North or South Pole. During a single communication session the declination angle changes from $\delta_{0}$ to $\delta_{0}+\delta_{\mathrm{com}}$ with $\delta_{0}$ being the initial declination of the satellite when it appears on the observer horizon and

[TABLE]

Here $t_{\mathrm{com}}$ is the maximal duration of the communication session, $H$ is the satellite altitude above the ground, and $R_{\oplus}$ is the Earth radius. Finally, the inclination angle $\Delta\iota$ we define as the angle between the observer meridian plane and the satellite orbit plane and the geographical latitude of the observer we denote as $\Psi$. The set $\{\Psi,\delta,\Delta\iota\}$ determines the instantaneous position of the satellite relative to the observer in the geocentric coordinate system and is the alternative parameter set to the set $\{A,Z\}$ in the Earth surface coordinate system.

The following relation between angles $Z$, $\delta$, $\Delta\iota$, $\Psi$ shown in Fig. 8 can be established using spherical trigonometry

[TABLE]

[TABLE]

Here $\delta$ changes during the communication session from $\delta_{0}=-\arctan(\cos\Delta\iota\cot\Psi)$ to $\delta_{0}+\delta_{\mathrm{com}}$. We note that the finite value of the inclination angle $\Delta\iota$ in (84) restricts the definition domain of $Z$ to $[Z_{\mathrm{min}}^{\Delta\iota},\pi/2]$ for the observer located at the geographical latitude $\Psi$. The minimal value of $Z$ can be found from Eq. (84) as

[TABLE]

which corresponds to the satellite declination angle (counted from the Equator)

[TABLE]

The knowledge of the satellite declination angle, the orbit inclination angle, and geographical position of the observer allows one to determine the instantaneous values of zenith angle from Eqs. (84) as well as to determine the slant range.

A.2 Inclined and zenith orbits

In most practical cases the satellite trajectory is inclined relative to the observer zenith direction. For example, the initial trajectory with zero inclination to the observer’s zenith becomes inclined when the satellite makes one or more revolutions. Let us consider the satellite whose orbit initially has zero inclination relative to the meridian plane of the observer $\widetilde{O}$ (see Fig. 9). Let us also assume that initially the satellite is positioned at observer’s zenith, i.e., at point $S^{\prime}$. For simplicity we consider the ideal polar orbit that passes through both Earth’s celestial poles. After time $T_{\mathrm{sat}}$, i.e., the satellite orbiting period, the satellite is positioned again in the point $S^{\prime}$. Meanwhile the observer position $\widetilde{O}$ is moved to the point $O$ due to Earth’s rotation. The resulting inclination angle reads as

[TABLE]

where $v_{\oplus}$ is the speed of Earth’s rotation at the Equator. Here and in the following we assume that the change of the inclination angle during the satellite transition over the observer horizon is relatively small, so that we can approximate it being constant (87) during the whole communication session. The negative values of $\Delta\iota$ refer to situations when the visible satellite trajectory is positioned eastwards to the observer. We calculate the slant range between $O$ and $S$ along the satellite orbit provided the trajectory $S^{\prime}S$ of the moving satellite lies above the observer horizon $\mathrm{H}$.

The instantaneous position of the satellite relative to the observer is determined by zenith $Z$ and azimuth $A$ angles as shown in Fig. 8. For satellite tracking purposes both these angles are of importance. In the context of optical satellite-based communication we are interested primarily in the instantaneous value of the slant range as a function of zenith angle. Indeed, by assuming that tracking systems point correctly the telescope towards the flying satellite, we can can ignore the dependence of the slant range on the azimuth angle. Then, the slant range $L$ for the inclined orbit can be determined from the triangle $OSC$ in Fig. 2 using the law of cosines

[TABLE]

Solving this equation with respect to $L$ we obtain

[TABLE]

The root of the quadratic equation (88) is chosen to yield the non-negative slant range. At the observer’s horizon the slant range reaches its maximal value L(90^{\circ})$$=\sqrt{(R_{\oplus}+H)^{2}-R_{\oplus}^{2}}\approx\sqrt{2HR_{\oplus}}. The shortest slant range is determined as $L(Z_{\mathrm{min}}^{\Delta\iota})$, where the minimal zenith angle is given by Eq. (85).

The special case of the slant range (89) corresponds to the zenith satellite orbit with zero inclination angle $\Delta\iota$. Some authors refer to this orbit as to the ”best pass” Bourgoin2013 since the slant range, that we denote as $L_{0}$, gives the shortest distance between the observer and the satellite at zenith, namely, $L_{0}(0^{\circ})=H$. This fact makes the zenith orbit as well as the orbits with small inclination angle $\Delta\iota$ most attractive for the optical communication. Denoting the zenith angle within the zenith orbit plane (observer’s meridian plane) as $Z_{0}$, it follows from Eq. (84) that $Z_{0}=\Psi-\delta$. Consequently, we have L_{0}(Z_{0})=L(Z)\bigl{|}_{\Delta\iota=0}.

Appendix B Standard Atmosphere

The standard atmosphere is an idealized model of the Earth’s atmosphere for heights ranging from the surface to 1000 km StandardAtmosphere . The model yields the air density, viscosity, etc., for various altitudes. For our purposes, the most useful values calculated within the standard atmosphere model are the temperature $T$ and the pressure $P$. These variables are important for the calculation of the refractive index variation with altitude.

The altitude dependence of the temperature is approximated with linear segments

[TABLE]

Each segment lies within an atmospheric layer bounded by the surfaces with the altitudes $H_{b-1}$ and $H_{b}$. Table  4 summarizes the reference values $T_{b}$ and the vertical gradient values of temperature (lapse rates)

[TABLE]

which enter Eq. (90). The altitude dependence of pressure can be found from the gas law and hydrostatic equation and reads as

[TABLE]

for constant lapse rate and

[TABLE]

for isothermal layers ($\lambda_{b}=0$). Here $g=9.8~{}\text{m/s}^{2}$ is the gravitational acceleration and $R=287.053~{}\text{J/kg}\cdot\text{K}$ is the gas constant for air. The reference values of pressure $P_{b}$ are given in Table 4.

Following Birch and Downs Birch1994 we adopt the revised form of the Edlén equation for the atmospheric refractive index $n$,

[TABLE]

where $(n-1)_{s}$ is given by the dispersion equation

[TABLE]

Here $\lambda$ is the optical wavelength given in $\mu$m. Since the values of the refractive index are distinct from the vacuum value $n=1$ up to the mesosphere, we restrict our attention only to altitudes ranging from 0 to 85 km.

Figure 10 shows the variation of temperature, pressure and refractive index with altitude. The temperature profile shows several layers where temperature dependence on altitude can be approximated by linear relations. We use these specific altitudes for the calculation of the refractive index profile using linear segmentation. Such a procedure simplifies the calculation of atmospheric refraction but introduces errors that grow for large zenith angles $Z\sim 90^{\circ}$. The altitude dependence of the refractive index within the $i$th segment reads as

[TABLE]

Table 5 summarizes the corresponding values of the refractive index values and gradients for the corresponding segments.

For the calculation of the deterministic atmospheric extinction, the number density of absorbing and scattering particles is required. The standard atmosphere model gives the following dependence of the relative number density on altitude $h$ (given in meters) Duntley1948 ,

[TABLE]

Here $N_{0}$ is the number density at the observer level ($N_{0}=2.55\times 10^{25}$ m*-3* at sea level) and $\mathcal{H}_{0}=6\,600$ m.

Appendix C Path elongation due to atmospheric refraction

In the presence of Earth’s atmosphere, the slant path $L$ is elongated due to refraction on interfaces of atmospheric layers with different refractive indices (see Fig. 11). In the real atmosphere the gradient of the refractive index is a continuous function of height. We use the standard atmosphere model and consider 10 atmospheric layers (c.f. Appendix B, Table 5). The refractive index is linearly segmented, such that it is a linear function of height within one layer [cf. Eq. (96)].

As a consequence of atmospheric refraction the apparent zenith angle $Z_{a}$ starts to deviate from the true zenith angle $Z$. These angles are related with each other as

[TABLE]

where $n_{0}=1.00027$ is the refractive index of the lower layer of atmosphere Birch1994 . Near sea level the two zenith angles differ by approximately a minute of arc at $Z=45^{\circ}$ and half a degree near the horizon.

Refraction of optical rays happens on the $i$th interface of two adjoined $i$th and $(i+1)$th layers of heights $H_{i}$ and $H_{i+1}$, correspondingly, and is characterized by the angle of incidence $\pi/2-\beta_{i}$ and the angle of refraction $\psi_{i}$ (see Fig. 11). If the ray enters the observer telescope at the zenith angle $Z_{a}$, Snell’s law yields the geometric invariant

[TABLE]

where

[TABLE]

is the elevation angle at the observer $O$ and $R_{\oplus}$ is the Earth radius. Using the notation

[TABLE]

with $H$ being the altitude of the satellite above the ground we derive from Eq. (99)

[TABLE]

where the index $N$ corresponds to the last atmospheric layer at altitude 85 km which is accounted in our calculations (in our case $N=10$). We also use the index $i=1,...,N$ to denote the atmospheric layers above the ground level.

Using simple geometric considerations and Snell’s law, we derive the following relations for the angles relevant for the calculation of the ray path length,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The remaining angles to be determined are $\delta_{i}$ and $r_{i}$. The former angle is associated with the local elevation angle error and can be found from the law of sines

[TABLE]

Solving these equations with respect to $\delta_{i}$ we derive

[TABLE]

[TABLE]

The remaining angle $r_{i}$ is the bending angle within an $i$th layer which is calculated from the refraction integral Mahan1962 ,

[TABLE]

The total angle $r=\sum_{i}r_{i}$ is known in optical astronomy as atmospheric refraction. It can be shown that the refraction integral (113) can be simplified to the following approximate analytical form,

[TABLE]

provided that $H_{i}-H_{i-1}\ll H_{i-1}$ and if the refractive index gradient ${\rm{d}}n/{\rm{d}}h$ is constant within the atmospheric layer.

The length of the optical ray path inside the $i$th layer can be found from the triangle $CP_{i-1}P_{i}$ applying the law of cosines

[TABLE]

[TABLE]

For the optical path length in vacuum we find from the triangle $CP_{N}S$

[TABLE]

Herein we need the angle

[TABLE]

The ratio of the total length of the ray trajectory to the geometric path length [c.f. Eq. (1)] reads as,

[TABLE]

which is the elongation factor due to atmospheric refraction.

The elongation factor (119) is calculated without taking into account the finite radius of curvature of optical rays. Actually, the beam curvature near the ground is the smallest one, i.e. the ray bending is the greatest. In this case the empirical function for ray curvature $\mathcal{K}$ is given by Gaifillia2016 ; Hirt2010

[TABLE]

where $T$ is the temperature in K, $P$ is the pressure in mb, $\lambda(h)={\rm{d}}T/{\rm{d}}h$ is the lapse rate in K/km. For the standard atmosphere at $Z_{a}=90^{\circ}$ the curvature reaches its minimum value of $4.4R_{\oplus}$. Since this value is still larger than the Earth’s radius, we can assume that the optical rays within each atmospheric layer can be considered to be straight lines.

Appendix D Atmospheric model of optical turbulence

Dewan et al. Dewan1993 proposed a simplified version of the multiparameter VanZandt model VanZandt1978 for the refractive index structure parameter variation with altitude. This so-called Air Force Geophysics Laboratory (AFGL) model utilizes the meteorological data derived from radiosondes and yields the following analytic formula for the refractive-index structure parameter profile:

[TABLE]

Here the factors $(0.1)^{4/3}10^{Y(h)}$ determine the outer scale of turbulence in a statistical manner. The function $Y(h)$ is empirically related with the altitude-dependent wind shear $S(h)$ and the lapse rate $\lambda(h)$ as Jackson2004

[TABLE]

The parameter $M$ is connected with the gradient of the index of refraction,

[TABLE]

where $P$ is the pressure in mb, the temperature $T$ is given in K, the buoyancy frequency (Brunt-Väisälä frequency) in s*-1*. The buoyancy frequency reads as

[TABLE]

where $g$ is the acceleration of gravity, and $\gamma=9.8$ K/km is the dry air adiabatic lapse rate. The altitude dependencies of pressure and temperature are governed by Eqs. (90), (92), and (93).

The altitude variation of the refractive index structure parameter can be determined from Eq. (121), provided the altitude variation of the wind shear $S(h)$ is known. In this study we assume horizontal homogeneity of the atmosphere, which means that the mean wind properties do not depend on the horizontal position of the observer. Thus, we assume a flat terrain and neglect any spatial inhomogeneity of Earth’s surface. This assumption is equivalent to the independence of the vertical wind component on altitude Wyngaard . For the determination of the vertical wind shear components due to meridional and zonal winds we use the HWM93 thermospheric wind model Hedin1996 . The examples of winds and corresponding shear components are shown in Fig. 12 for summer night and the observer located near Munich Program .

For the atmospheric boundary layer close to the Earth’s surface, the refractive index structure parameter given by Eq. (121) is not applicable. We use the Walters and Kunkel model (WK) Walters1981 for the $C_{n}^{2}$ profile within the boundary layer,

[TABLE]

for nighttime and

[TABLE]

for daytime. Here $h_{i}$ is the height of the inversion layer above the ground ($h_{i}\sim 0.5$ km at nighttime and $\sim 1$ km at daytime), $h_{0}$ is the reference height referred to as the Monin-Obukhov scale ($h_{0}\sim 10$ m at nighttime and $\sim 5$ m at daytime).

Figure 13 shows the profile of the refractive-index structure parameter calculated using the AFGL and the WK models [cf. Eqs. (121), (127), and (131)]. For comparison, the widely used Hufnagel model Hufnagel1964 as well as the estimation of $C_{n}^{2}$ based on thermosonde data Frehlich2010 are shown. The AFGL+WK model agrees better with the experimental profile than the Hufnagel model. The difference between nighttime and daytime profiles is negligibly small for the altitudes above the inversion layer height $h_{i}$. Figure 12 shows the different behavior of the wind shear for meridional and zonal winds. However this anisotropy has almost no influence on the $C_{n}^{2}$ profile. As a consequence, we could neglect the dependence of this turbulence characteristics on the direction of view of the observer.

Finally, we relate the vertical profile $C_{n}^{2}(h)$ to the corresponding profile along the slant range. For the uplink configuration, the relation between the height $h$ of the certain point along the slant path and the distance between this point and the light source, $L_{r}\xi$, $\xi\in[0,1]$, can be found from the law of cosines [cf. Ref. Gurvich and see Fig. 2 (a)]

[TABLE]

Similarly, for the downlink we derive by replacing $\xi\rightarrow 1-\xi$

[TABLE]

Under the condition $L_{r}/R_{\oplus}\ll 1$, these equations reduce to

[TABLE]

[TABLE]

The refractive index structure parameter given by Eqs. (121), (127), and (131) maps for the uplink slant range as $C_{n}^{2}(h)\rightarrow C_{n}^{2}(h^{\mathrm{UL}})=C_{n}^{2}(L_{r},\xi)$ and for the downlink as $C_{n}^{2}(h)\rightarrow C_{n}^{2}(L_{r},1-\xi)$.

Appendix E Approximation of Eq. (38) for satellite-mediated links

In this appendix we give the approximation for the integral kernel (38). We consider the downlink configuration and just note that the formulas for uplink are obtained in the similar footing with the variable replacement $\xi\rightarrow 1-\xi$. For the case of strong turbulence or long propagation distances in turbulence the coherence radius $\rho_{0}$ [cf. Eq. (25)] is small. This condition is well satisfied for satellite-mediated atmospheric links. In this case, the exponential in Eq. (38), $\mathcal{J}(\mathbf{r},\mathbf{r}_{1}^{\prime},\mathbf{r}_{2}^{\prime},\mathbf{r}_{3}^{\prime})$, differs significantly from zero in the following regions:

[TABLE]

[TABLE]

[TABLE]

The function (38) can be approximated then as linear combination of three terms Banakh1979 ; Charnotskii2010

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The first term (140) accounts for the contributions from regions (136) and (138). The term (141) accounts for the regions (137) and (138), while the term (142) eliminates the double counting of the region (138) in the integral kernel $\mathcal{J}$. In the first approximation ($n=0$) we have

[TABLE]

which already gives a good approximation to the kernel $\mathcal{J}$ (cf. Ref. Banakh1979 ).

Appendix F Phenomenological model of aperture-averaged scintillations

The rigorous analysis of satellite-mediated quantum communication links should account for turbulent disturbances for a wide range of the zenith angle. Here we calculate the aperture averaged scintillation index with accounting for saturation effects. The aperture-averaged scintillation index can be obtained by substitution of Eq. (28) in Eq. (27), which leads to Tatarskii2016 333Originally, the expression for the scintillation index for the log transmittance, $\sigma_{\log\eta}^{2}$, is derived in the chapter 13 of Ref. Tatarskii2016 . Equation (144) is then obtained from the relation $\sigma_{\eta}^{2}=\exp[\sigma_{\log\eta}^{2}]-1$ (cf. Ref. Andrews2005 ).

[TABLE]

where

[TABLE]

is the aperture filter function, and $\langle I\rangle$ is the mean intensity. The Fourier transformed correlation function of the intensity fluctuations is given by

[TABLE]

For a circular aperture of radius $a$ the filter function (145) is easily evaluated to be

[TABLE]

where $\kappa=|\boldsymbol{\kappa}|$ and $J_{n}(x)$ is a Bessel function of the first kind.

In the limiting case of the vanishingly small Fresnel number $\Omega=kW_{0}^{2}/(2L_{r})\rightarrow 0$ the collimated beam with the initial beam spot size $W_{0}$ can be considered as a plane wave. Indeed, for the satellite-mediated link under consideration the slant range $L_{r}$ is large enough even at the zenith and $\Omega\rightarrow 0$. In this limiting case the intensity spectral density $F_{I}(\boldsymbol{\kappa},L_{r})$ can be obtained by solving the equations of geometrical optics Tatarskii2016 444Reference Tatarskii2016 derives the formula for the logarithmic amplitude of the light wave $\chi=\log(A/A_{0})$. Under weak fluctuation conditions the relation $F_{I}(\boldsymbol{\kappa},L_{r})=4\langle I\rangle^{2}F_{\chi}(\boldsymbol{\kappa},L_{r})$ is valid Andrews2005 .

[TABLE]

where the turbulence spectrum $\Phi_{n}$ is given in Eq. (16). For the Kolmogorov turbulence spectrum (18) the intensity spectral density $F_{I}(\boldsymbol{\kappa},L_{r})$ is proportional to $\kappa^{1/3}$ in the inertial range. Thus, for large values of the transverse wave number the intensity spectral function diverges. In order to remedy this nonphysical effect we introduce the cutoff of higher spatial frequencies $\Delta\kappa$ such that $\kappa\in[0,\Delta\kappa]$. The specifical choice of $\Delta\kappa$ will be discussed later in this appendix. Taking into account this cutoff effect the aperture averaged scintillation index is calculated straightforwardly,

[TABLE]

The obtained result shows that the intensity fluctuations are directly related to spectrum of turbulent fluctuations of the refractive index.

For the Kolmogorov turbulence spectrum (18) the integrations can be performed, provided the altitude dependence of structure constant is known. Here we adopt a simple dependence on the $z$ variable assuming that the strongest turbulence is near the Earth surface and decreases exponentially with the height Young1969 ,

[TABLE]

Here $H_{0}$ is the characteristic length at zenith and the factor $\sec Z_{a}$ accounts for the increase of this length for the slant paths. The integral with respect to $z$ in Eq. (149) is evaluated as

[TABLE]

In the limit of infinitely distant light source, $L_{r}\rightarrow\infty$, we get $\mathcal{I}(Z_{a})=2(H_{0}\sec Z_{a})^{3}$. Consequently, the scintillation index (149) for $\Delta\kappa\rightarrow\infty$ has the same dependence on the zenith angle as those obtained within the Rytov approximation, cf. Eq. (27).

The cutoff spatial frequency $\Delta\kappa$ in Eq. (149) can be estimated from the following considerations. The rigorous asymptotic analysis shows Fante1983 that for long propagation paths the covariance function $B_{I}(\boldsymbol{r},L_{r})$ in Eq. (148) differs from zero if $|\boldsymbol{r}|\leq\Delta r$, where

[TABLE]

Here, $\lambda$ is the optical wavelength, $L_{\mathrm{turb}}$ is the propagation path length in turbulence, and $\rho_{0}$ is the spatial coherence length given by Eq. (25). According to Eq. (150) the turbulence is sufficient for heights $h\leq H_{0}$. Hence, without loss of generality we set $L_{\mathrm{turb}}=H_{0}\sec Z_{a}$. The finiteness of the definition domain for the covariance function $B_{I}(\boldsymbol{r},L_{r})$ introduces the cutoff frequency $\Delta\kappa$ for the spectral density $F_{I}(\boldsymbol{\kappa},L_{r})$. Using the relation of Fourier analysis $\Delta\kappa\Delta r\geq\mu$, we may estimate $\Delta\kappa=\mu/\Delta r$, where $\mu\sim 1$ is a phenomenological parameter. It is easy to see that $\Delta\kappa\propto(\sec Z_{a})^{-8/5}$ and decreases while the zenith angle approaches the horizon. For further insight into physical meaning of the parameters $\Delta r$ and $\Delta\kappa$ we point the reader’s attention on Ref. Andrews1999 . Finally, performing the integration of (149) over $\kappa$ and over $z$ variable [cf. Eq. (F)] we obtain Eq. (49), where only the first term of $\mathcal{I}(Z_{a})$ is accounted. At the cost of analytic simplicity, the obtained result can be further generalized for the case of Cassegrain-type apertures with the inner circular obscuration Vasylyev2018a .

Appendix G Decoy state: statistical fluctuation analysis

In this appendix, we summarize the method for estimation the lower (upper) bounds of single-photon gain and error rate taking into account the finite-key effects Zhang2017 . On quantum state preparation stage, Alice generates each bit in her raw key by randomly choosing the encoding basis ($X$ or $Z$) and the intensity (corresponding to the vacuum state, weak decoy state, and signal state). The total number of bits sent by Alice to Bob is given by

[TABLE]

where the superscripts correspond to signal ($s$), weak-decoy ($d$), and vacuum ($v$) states. We denote $q^{a}=N^{a}/N$ as the rate with which Alice encodes a state with intensity $\mu_{a}$, $a=s,d,v$. The conditional probability that an $i$-photon state corresponds to a coherent pulse with the intensity $\mu_{a}$ is

[TABLE]

where the approximation sign is due to statistical fluctuations. The finiteness of numbers of generated bits yields

[TABLE]

for the $q$ parameter in Eq. (71).

Bob measures the received states in $X$ or $Z$ basis chosen randomly. After basis reconciliation and key sifting, Bob possesses the total number of sifted bits

[TABLE]

For the atmospheric quantum channel, the bit numbers $N$ and $M$ are related as

[TABLE]

where the factor $1/2$ is due to the sifting procedure, $\eta_{\mathrm{d}}$ is the deterministic channel loss including the detector efficiency given by Eq. (75), and $\eta$ is the fluctuating channel transmittance. The bit number corresponding to the $i$-photon state is determined as

[TABLE]

where $p_{1}^{a}$ is the same as the probability (154) chosen by Alice.

After error correction and error verification the secure key rate depends on a lower bound, on the gain of single-photon components of the signal state $Q_{1}^{L}$ and an upper bound, on the corresponding error rate $e_{1}^{U}$. The estimation of $Q_{1}^{L}$ and $e_{1}^{U}$ should be performed in each basis. The corresponding components we denote by superscripts $x$ and $z$ referring to the $X$ and $Z$ bases, correspondingly. For the $Z$ basis the bounds on the single-photon gain and the error rate are given by

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

Here $e_{0}=1/2$ and the lower and upper bounds

[TABLE]

[TABLE]

[TABLE]

[TABLE]

are estimated by using the Chernoff bound method and are related to the overall gain (73), with the overall quantum bit error rate (77) components, and with the counting rate for vacuum decoy states and dark count contributions

[TABLE]

The deviation function

[TABLE]

can be determined for the specified failure probability $\varepsilon$.

Finally, the upper bound $\theta^{U}$ in (160) is obtained by numerically solving

[TABLE]

with respect to $\theta^{U}$  Fung2010 ; Zhang2017 . Here, the averaging $\langle...\rangle_{\mathrm{tr}}$ is performed according to Eq. (70) and the lower bounds of sifted key numbers are

[TABLE]

[TABLE]

with $q^{a}=N^{a}/N$. In Eq. (169) the exponential function $\xi(\theta)$ read as Fung2010

[TABLE]

where $q^{x}{=}\langle M_{1}^{xL}\rangle_{\mathrm{tr}}/(\langle M_{1}^{xL}\rangle_{\mathrm{tr}}{+}\langle M_{1}^{szL}\rangle_{\mathrm{tr}})$ is the bias ratio.