Satellite-based links for quantum key distribution: beam effects and weather dependence

Quantum key distribution (QKD) and quantum communication in general have the potential to revolutionise the way we communicate confidential information over the internet. The natural carriers for quantum information are photons, that are already widely used in classical networks of optical fibers to achieve high communication rates. Unfortunately, even though enormous improvements have been obtained in the last years [1, 2], scaling quantum communication protocols over long distances is very challenging, due to the losses experienced during the propagation inside the optical fibers. Several schemes for the realization of quantum repeaters have been proposed in recent years, that could allow to bridge long distances and naturally be implemented inside a quantum communication network [3–7]. Considering the important technological hurdles that quantum repeaters should overcome before becoming useful, satellite-based free-space links look like the most practical way to achieve long-distance QKD in the short term [8]. They can take advantage of the satellite technology and the optical communication methods developed in the last decades in the classical case. Various feasibility studies had addressed this topic in the last twenty years [8–11] and several experiments have definitely proved that the technology involved is ready for deployment [12–16].

Optical satellite-based links have the important drawback of being strongly dependent on the weather conditions [17–20]. The presence of turbulent eddies and scattering particles like haze or fog generates random fluctuations of the relative permittivity of the air, on different length- and time-scales. This phenomenon affects the light propagation in a complicated way, inducing random deviations and deformations of any optical beam sent through the atmosphere. It results in reduced transmittance, because of geometrical losses due to the finite collection aperture, and random modifications of the phase front. A comprehensive model of these effects is then necessary, in order to precisely evaluate the performance of the link when used for quantum communication protocols.

In this work we generalize the approach proposed in [21, 22] to satellite-based links and we evaluate their losses in several practical cases, under different weather conditions. This information is then used to assess the performance of the link in terms of the achievable key rates using different implementations of QKD. The case of Low Earth Orbit (LEO) satellites is addressed, assuming different payloads and sizes of the optical elements.

This work is organized as follows. In section 2 we introduce the problem of free-space optical links and an analytical method to study them. The discussion continues in appendix A. In section 3 a detailed description of the model used to simulate the satellite-based link is presented. Then, the main results are shown and discussed, together with pros and cons of our approach. In section 4 we use the analysis of the transmittance of the channel conducted in the previous section to study the key rate achievable by different QKD implementations, in some interesting real-life scenarios. The analysis concerning the use of smaller and more affordable satellites is performed in section 5. Finally, the results are summarized and discussed in section 6. The appendix starts with a recap of the results of [21, 22] (appendix A) and their application to the problem at hand (appendix B). Then two models for the estimation of the stray light satellite links [11, 23] are presented in appendix C. Appendix D is devoted to the definition of the QKD protocols we use in section 4 and the expression of the correspondent key rates. In appendix E we report the parameters chosen for the simulations and we discuss their pertinence.

The problem that we address in the first part of the work is the following. A Gaussian beam is sent, either from an orbiting transmitter or from a ground station, through a non-uniform link partially inside the atmosphere and partially in vacuum. We are interested in the transmittance of the received beam through a circular aperture of radius a (the receiving telescope)

$$\begin{eqnarray}&&\eta ={\int }_{| {\boldsymbol{\rho }}{| }^{2}={a}^{2}}{{\rm{d}}}^{2}{\boldsymbol{\rho }}\ | u({\boldsymbol{\rho }},L){| }^{2},\end{eqnarray} \tag{ 1 }$$

which is a random variable, because of the intrinsic randomness of the fluctuations in the medium. Here $u({\boldsymbol{\rho }},L)$ is the beam envelope at the receiver plane (at distance L from the transmitter, with ρ the position in the transverse plane).

The so-called Elliptic Beam Approximation [21] greatly simplifies the analysis: the atmosphere is assumed to generate only:

• deflection of the beam as a whole (Beam Wandering);
• elliptic deformations of the beam profile;
• extinction losses due to back-scattering and absorption.

In this case the state of the beam at the receiver plane is completely described by the vector of parameters (refer to figure 1)

$$\begin{eqnarray}&&{\bf{v}}=({x}_{0},{y}_{0},{W}_{1},{W}_{2},{\varphi }_{0}),\end{eqnarray} \tag{ 2 }$$

representing the beam-centroid coordinates, the principal semi-axes of the elliptic profile and the angle of orientation of the ellipse. The transmittance is then a function of these beam parameters and the radius of the receiving aperture.

Zoom In Zoom Out Reset image size

The fluctuations of the relative permittivity of the atmospheric air can be statistically modeled [24–32]. The probability distribution of the parameters in equation (2) can then be analytically estimated, as shown in [21, 22]. A brief recap of the derivation and the main results is presented in appendix A. This allows, through random sampling, to obtain the probability distribution of the transmittance (PDT), an important figure of merit for fluctuating links. This approach gives no information about the phase of the wavefront, but this is not a problem when phase-insensitive measurements are considered (e.g. the BB-84 QKD protocol that we analyze in section 4).

The atmosphere can in general be divided into several layers, depending on the properties of different physical parameters, like density of the air, pressure, temperature, density of ionized particles, and so on. This structure is site-dependent, especially regarding the thickness of the different layers. For this reason, in this work we assume a simplified version of a satellite-based optical link: a uniform atmosphere up to a certain altitude $\bar{h}$, then vacuum all the way up to the satellite (at altitude $\bar{L}$), as pictured in figure 2. Instead of a continuum of values describing the physical quantities as a function of the altitude, we now have only two parameters, namely the value of the quantity inside the uniform atmosphere and the effective thickness $\bar{h}$. This is likely to be a good approximation, because the atmospheric effects are prominent only in the first 10–20 km from the ground, while usual orbit height for LEO satellites are above 400 km. For the remainder of the paper we choose a minimum altitude of the satellite $\bar{L}=500$ km, achieved exactly above the ground station. In this case, the extension of the orbit of the satellite which can be usable for key distribution corresponds roughly to the interval $L\in [500,2000]$ km, corresponding to angles from the zenith in the interval [0°, 80°]. The effective thickness of the atmosphere $\bar{h}$ is fixed here to 20 km, for the considerations above.

Zoom In Zoom Out Reset image size

As introduced in section 2, we want to generalize the model proposed in [21, 22] to the just described case of a non-uniform link between the ground and a satellite. The computation follows the same steps and is described in appendices A and B. First of all we need to evaluate equations (A.13), (A.14) and (A.15) in order to compute the moments of the distributions of the elliptic beam parameters (equation (2)). To do so, an integration along the propagation path must be performed (equations (A.19) and (A.20)). Here we introduce the considerations of the previous paragraph, imposing that the parameters measuring the strength of the atmospheric effects are constant (greater than 0) inside the atmosphere and 0 outside. In particular we assume

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{Down}-\mathrm{links}\qquad {C}_{n}^{2}(z) & = & {C}_{n}^{2}\ {\rm{\Theta }}(z-(L-h))\\ {n}_{0}(z) & = & {n}_{0}\ {\rm{\Theta }}(z-(L-h))\\ \mathrm{Up}-\mathrm{links}\qquad \quad {C}_{n}^{2}(z) & = & {C}_{n}^{2}\ {\rm{\Theta }}(h-z)\\ {n}_{0}(z) & = & {n}_{0}\ {\rm{\Theta }}(h-z),\end{array}\end{eqnarray} \tag{ 3 }$$

where C_n² is the value of the refractive index structure constant and n₀ is the density of scattering particles. Θ(z) is the so-called Heaviside step-function, z is the longitudinal coordinate, L is the total length of the link and h is the length traveled inside the atmosphere, as shown in figure 2. A down-link corresponds to the situation of satellite-to-ground communication, so the atmospheric effects kick-in only for z > (L − h) (final section of the propagation), while for up-links it is limited to z < h. We remark that some models for the altitude-dependence of the optical quantities, like C_n², are available in the literature [33–37], but they are correct only in the geographical site and in the atmospheric conditions in which they had been experimentally extracted (more details in appendix E). Additional extinction losses due to back-scattering and absorption in the atmosphere are modeled by a parameter χ_ext, as described in appendix A. Its value is adjusted from the analysis performed in [10] based on the MODTRAN5 software [38]. In this model, the values of C_n² and n₀ completely describe the atmospheric conditions together with the thickness $\bar{h}$ and the extinction factor χ_ext.

Following the analysis of appendix A (in particular equations (A.13), (A.14), (A.15)), we compute the first and second moments of the beam parameters in equation (2) for the link described in equation (3). The distribution of the angle of orientation of the elliptical profile φ₀ is assumed uniform in [0, π/2] as in [21, 22]. The mean value and variance of the beam centroid position for up-links are the same for x and y directions and equal to (equation (A.13))

$$\begin{eqnarray}&&\langle {x}_{0}\rangle =\langle {y}_{0}\rangle =0,\qquad \langle {x}_{0}^{2}\rangle =\langle {y}_{0}^{2}\rangle =0.419\ {\sigma }_{R}^{2}\ {W}_{0}^{2}\ {{\rm{\Omega }}}^{-\tfrac{7}{6}}\displaystyle \frac{h}{L},\end{eqnarray} \tag{ 4 }$$

where the quantity ${\sigma }_{R}^{2}=1.23\ {C}_{n}^{2}\ {k}^{\tfrac{7}{6}}{L}^{\tfrac{11}{6}}$ is the so-called Rytov parameter and ${\rm{\Omega }}=\tfrac{{{kW}}_{0}^{2}}{2L}$ is the Fresnel number. The condition $\langle {x}_{0}\rangle =0$ is achieved by proper pointing. The first two moments of the semi-axes of the ellipse squared, W²_i with i = 1, 2, are instead estimated from equations (A.14) and (A.15)

$$\begin{eqnarray}&&\langle {W}_{i}^{2}\rangle =\displaystyle \frac{{W}_{0}^{2}}{{{\rm{\Omega }}}^{2}}\left(1+\displaystyle \frac{\pi }{8}\ L\ {n}_{0}\ {W}_{0}^{2}\displaystyle \frac{h}{L}+2.6\ {\sigma }_{R}^{2}\ {{\rm{\Omega }}}^{\tfrac{5}{6}}\displaystyle \frac{h}{L}\right),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\langle {\rm{\Delta }}{W}_{i}^{2}{\rm{\Delta }}{W}_{j}^{2}\rangle =(2{\delta }_{{ij}}-0.8)\ \displaystyle \frac{{W}_{0}^{4}}{{{\rm{\Omega }}}^{\tfrac{19}{6}}}\left(1+\displaystyle \frac{\pi }{8}\ L\ {n}_{0}\ {W}_{0}^{2}\displaystyle \frac{h}{L}\right){\sigma }_{R}^{2}\displaystyle \frac{h}{L}.\end{eqnarray} \tag{ 6 }$$

Similar expressions hold for down-links, for the beam centroid position

$$\begin{eqnarray}&&\langle {x}_{0}\rangle =\langle {y}_{0}\rangle =0\qquad \langle {x}_{0}^{2}\rangle =\langle {y}_{0}^{2}\rangle =\alpha \ L\,\end{eqnarray} \tag{ 7 }$$

and for the semi-axes of the elliptical profile

$$\begin{eqnarray}&&\langle {W}_{i}^{2}\rangle =\displaystyle \frac{{W}_{0}^{2}}{{{\rm{\Omega }}}^{2}}\left(1+\displaystyle \frac{\pi }{24}\ L\ {n}_{0}\ {W}_{0}^{2}{\left(\displaystyle \frac{h}{L}\right)}^{3}+1.6\ {\sigma }_{R}^{2}\ {{\rm{\Omega }}}^{\tfrac{5}{6}}{\left(\displaystyle \frac{h}{L}\right)}^{\displaystyle \frac{8}{3}}\right),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\langle {\rm{\Delta }}{W}_{i}^{2}{\rm{\Delta }}{W}_{j}^{2}\rangle =(2{\delta }_{{ij}}-0.8)\displaystyle \frac{3}{8}\displaystyle \frac{{W}_{0}^{4}}{{{\rm{\Omega }}}^{\tfrac{19}{6}}}\left(1+\displaystyle \frac{\pi }{24}\ L\ {n}_{0}\ {W}_{0}^{2}{\left(\displaystyle \frac{h}{L}\right)}^{3}\right){\sigma }_{R}^{2}{\left(\displaystyle \frac{h}{L}\right)}^{\displaystyle \frac{8}{3}},\end{eqnarray} \tag{ 9 }$$

where α ∼ 2 μrad is the angular pointing error.

There are two main differences between the expressions related to the up-link and down-link configurations. First, they depend on a different power of the ratio $\tfrac{h}{L}$. As $\tfrac{h}{L}\ll 1$, we deduce, as expected, that the atmospheric effects are much stronger for up-links than for down-links. The phenomena involved here (beam deflection and broadening) are angular effects, whose contribution on the final size of the beam (and thus, on the losses of the channel) are proportional to the distance traveled after the 'kick in' of the effect. For up-links, these effects happen very close to the transmitter, and then the beam broadens for hundreds of km before being detected. In the down-link scenario, instead, the beam travels in vacuum for the largest portion of the distance, and the atmospheric effects take place only at the end of the propagation, in the last tens of km before the receiver. The second difference resides in the origin of the fluctuations of the beam centroid position x₀. For up-links, in fact, the deflections induced by the atmospheric effects are usually much stronger than the pointing error, which we neglect. For down-links, instead, at the top of the atmosphere the beam dimensions are already much larger than any turbulent inhomogeneity. In this case the induced beam wandering can be neglected and the pointing error becomes the main contribution.

The knowledge of the probability distribution of the elliptic beam parameters is then used to compute the PDT, through equation (A.22) and random sampling. Two examples are shown in figures 3 and 4 for a down-link and an up-link, respectively. The considerations of the previous paragraph can naturally be used to explain the difference in the shape of these two distributions. For down-links, especially at high elevation angles, like the case shown in figure 3, the value of the beam width at the receiver is comparable to the wandering induced by pointing errors. This means that it can happen that the beam wanders completely off the receiving aperture, giving values of transmittance close to 0. In the up-link case, instead, the beam broadening gets the upper hand: the beam at the receiver is so large that the wandering induced by the atmosphere cannot change the total transmittance very much. It results in a rather narrow distribution, peaked at much lower values of transmittance with respect to the down-link case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Now we want to study the expected loss introduced by the link as a function of the total link length. We show in figures 5 and 6 the mean value of the PDT as a function of the angle from the zenith and the total link length, for down-links and up-links, under different weather conditions. Every point in the graph has been obtained from 1000 samples of the parameters in equation (2) and using equation (A.22). The asymmetric nature of the PDT for some configurations of the link can make the use of the mean value partially misleading, however, the full PDT will be used in the next section to compute the secret key rates.

Zoom In Zoom Out Reset image size

The critical parameters here are, apart from the ones related to the atmospheric effects, the diameter of the sending and receiving telescopes and the signal wavelength. We chose D_sat = 30 cm for the orbiting one, D_grnd = 1 m for the ground station telescope and λ = 800 nm. These are demanding values, consistent with the Chinese mission Micius (see [12–15] for details). Further analysis are reported in section 5. We notice here that the assumption of perfect Gaussian beams sent by the transmitter is not very realistic. Standard telescopes generate beams with intensity distributions rather close to a circular Gaussian profile but with some imperfections, introduced for example by the truncation at the border of the optical elements. The main downside is that such beams will exhibit larger intrinsic beam broadening due to diffraction. In our model this effect can be taken into account by adjusting the value of the initial beam waist W₀, in order to match the far-field divergence expected from the imperfect quasi-Gaussian beam.

Our analysis confirms that, at least for the parameters chosen for the simulation, down-links are much preferable over up-links for quantum communication due to the smaller losses. However, up-links can still achieve losses below the threshold for the accomplishment of quantum communication tasks, QKD included. Particularly interesting is the comparison between night- and day-time operation. During the day, the higher temperatures bring stronger wind and more active mixing between the different layers of the atmosphere, leading to more pronounced turbulence effects. However, on average, during clear days the moisture content of the lower atmosphere is smaller than at night, resulting in weaker beam spreading due to scattering particles. At night, instead, the lower temperature results, on one hand, in a less turbulent atmosphere and, on the other, in the formation of haze and mist. In this situation, the contribution of scattering over such particulate can be stronger than the turbulence-induced effects.

We point out that in this analysis the path elongation due to refraction in the atmosphere has not been taken into account, as it gives substantial effects only at low elevation angles (≤10°). This and other zenith-angle-dependent effects have been thoroughly studied in [39].

Many different models for atmospheric channels and satellite-based links had already been proposed in the literature, due to the increasing interest in free-space optical communication. A comparison with them can highlight the strengths of the approach we reported in this section. Many feasibility studies [40] rely on models that average the intensity over sufficiently long times, so that the only atmospheric effect is overall a broadening of the beam. This approach gives no information on the PDT of the channel, that can be useful in many instances (for example, to apply post-selection techniques). A different approach has been chosen by [10], based on convolution between the beam envelope and the time-averaged pointing errors and beam broadening, leading again to no information about the PDT. A popular technique, that involves heavy numerical computations, is based on simulating the effect of the atmosphere by random phase screens regularly distributed along the propagation path in vacuum [41–43]. Many theoretical works have been devoted to find the analytical probability distribution that better fits the experimentally measured transmittance of free-space optical links. Mainly used are the log-normal [44, 45], Gamma–Gamma [46] and Double Weibull [47] distributions. Each of them appears to be more suitable depending on the strength of the turbulence, the length of the link and the configuration of the transmitting and receiving telescopes. On the contrary, the approach used here is a constructive method that allows to determine the PDT starting from the characteristics of the beam and the atmospheric conditions.

It has been shown that a post-selection of the time-intervals with greater transmittance can help to increase the secret key rates [48–50]: in this context, the ability of our approach to simulate not only the expected value of the transmittance, but its probability distribution too, may prove to be of great interest. In [51] the authors showed that the detrimental effects of asymmetric and fluctuating losses in measurement-device-independent QKD with decoy states can be counteracted by means of additional losses introduced by the central node. In this context, when the quantum links used are free-space, the information about the PDT allows to optimize such compensation losses to maximize the key rate.

Finally, we effectively take into account the contribution due to scattering particles, like fog or haze, making possible to model the effect of different weather conditions, a problem usually not addressed in previous works. It is particularly important during night-time operation, where a substantial amount of beam deformations can be imputed to scattering on moisture particles.

The transmittance shown in figures 5 and 6 can now be used to compute the expected secret key rates of a QKD protocol. In the following we analyze the performance of the BB-84 protocol [52] with polarization encoding, implemented using either a true single photon (SP) source or weak coherent pulses (WCPs). We use modern techniques to compute the secret key rates for SPs [53] and WCPs with decoy states [54–57], taking into account finite-key effects. The key rates are averaged over the PDT computed for different link lengths and configurations

$$\begin{eqnarray}&&\bar{R}={\int }_{0}^{1}R(\eta )\ { \mathcal P }(\eta )\ {\rm{d}}\eta =\displaystyle \sum _{i=1}^{{N}_{\mathrm{bins}}}R({\eta }_{i})\ { \mathcal P }({\eta }_{i}).\end{eqnarray} \tag{ 10 }$$

Zoom In Zoom Out Reset image size

Here $\bar{R}$ is the averaged key rate, R(η) the key rate at the specific value of the transmittance, ${ \mathcal P }(\eta )$ is the PDT. The integral average is approximated dividing the range [0, 1] in N_bins bins, centered in η_i for i = 1, N_bins, and taking the weighted sum of the rates. ${ \mathcal P }({\eta }_{i})$ is estimated through random sampling, as pointed out in section 3. The expressions for the key rates R(η) for the different implementations are given in appendix D, see equation (D.1) for SPs and equation (D.2) for WCPs.

The biggest source of noise in free-space optical links is represented by environmental light entering in the receiver telescope together with the signal photons. Simple models to estimate the amount of stray light [11, 23] are given in appendix C for down-links and up-links. In the following analysis we consider the number of stray photons to be independent of the position of the satellite. Particular situations concerning light pollution, like the presence of a city close to the ground station, may require a more specific model for low elevation angles.

The secret key rate resulting from a down-link and an up-link are reported in figures 7 and 8, for both night-time and day-time operation, under good weather conditions, corresponding to situation 1 in figure 5. We also report in the graphs the correspondent quantum bit error rate (QBER), defined in equation (C.4) in appendix C, averaged over the PDT. In the following we set the block sizes at 10⁶ for SPs and at 10⁸ for WCPs in down-link. This difference is justified by the higher repetition rates obtainable by modern WCP sources with respect to (still under development) true SP sources. Consider that the total link duration is around 300 s, corresponding to the complete passage of a LEO satellite over the ground station. In this time span, assuming a repetition rate of 10 MHz for SP sources and 1 GHz for WCP sources, several blocks of the size specified above can be exchanged in the down-link configuration. Due to the higher losses encountered in an up-link, the block size is lowered to 10⁵ for SPs and at 10⁷ for WCPs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At night it is possible to establish a non-zero key rate in down-link during the whole passage of the satellite in the SP implementation. Using WCPs, instead, the key rate drops to 0 when the satellite is around 20° over the horizon. In the daytime, instead, due to the stronger background light, the key rate vanishes at higher elevation angles, even considering improved spatial, spectral and temporal filtering. According to [23], the typical brightness of the sky background (see also section appendix C) in a clear day is about 3 orders of magnitude higher than during a full-moon night. We found that, in order to achieve a non-zero key rate for a reasonable portion of the transit, the filtering of the stray light must be tighter than during the night of a factor ≳100, obtained acting on the field of view (FOV) of the telescope and the width of the spectral filters (refer to table E3 in appendix E for details).

Up-links have poorer performance due to higher losses, but we are still able to distill a secret key with non-zero rates during the night, with slightly improved filtering (table E3 in appendix E). The SP implementation reaches almost the same range (in elevation angle) as the down-link configuration, while the difference with WCPs is greater because of the smaller block size. For day-time operation the stronger background light makes the QBER too high and the key rate vanishes, therefore we omit the corresponding graph. We stress that here (we refer to appendix C for details) we did not consider artificial light pollution. So these results reliably simulate only ground stations which are isolated and far from big cities.

Note that the finite key effects can be very detrimental when the number of exchanged signals becomes too small. Particular attention must be payed when up-links are considered. In order to reproduce the results reported in figure 8, the block length used in the security analysis is of the same order of magnitude of the number of signals exchanged during the whole passage of the satellite. This means that all the signals exchanged in a QKD session are processed in a single block in this case.

The simulations reported in sections 3 and 4 assume a quite demanding value of the optical aperture of the orbiting telescope. It is compatible with the Micius satellite [12–15], operated by the Chinese Academy of Science, as part of the Quantum Experiments at Space Scale research project. The complexity and high cost of the mission make the use of such big satellites unfeasible for the establishment of a world-wide quantum communication network. Many recent proposals foresee the use of nano-satellites (e.g. CubeSats [58–60]) for QKD implementation [9, 61–64]. The possibility to deploy many of such satellites in a single mission, to share the vector with other payloads and the modular nature lowers considerably the launch and building cost of these devices. They are usually loaded with smaller optics, of diameter ≤10 cm, even if larger apertures can be achieved with the use of deployable optics or bigger CubeSats (like the 12-Units satellite proposed in [63]). When used as transmitter, in the down-link configuration, the smaller aperture creates beams with much higher intrinsic divergence than the case studied in section 3. In the up-link configuration, instead, smaller transmittance is due to the smaller collecting area. We show in figure 9 the results of the link simulation for down-links and up-links, in good weather conditions.

Zoom In Zoom Out Reset image size

We see that the effect of the smaller optics diameter amounts to a difference in transmittance of about 5 dB for down-link and to 10 dB for up-links. Even though this result favors the down-link configuration even more, we have to take into account that a smaller aperture will collect not only less signal light, but also less stray light. The resulting QBER for up-links, then, will be almost independent of the diameter of the receiving telescope (compare figures 8 and 11), if other parameters, such as the FOV of the telescope, are kept fixed (see equation (C.1) in appendix C for details).

The key rates achievable for nano-satellites in the down-link and up-link configurations are reported in figures 10 and 11. As expected, the range of angles over which a non-zero key rate can be exchanged shrinks with respect to the case of section 4. We point out that we kept the block length fixed at the values reported in 4 even if, especially in the up-link configuration, that number of signals cannot be exchanged in a single transit of the Cube-sat.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We provide a general and fundamental model to simulate the losses introduced by a satellite-based optical link, useful for feasibility and performance analysis of future free-space QKD experiments. The ability to precisely evaluate the contribution due to different weather conditions will be crucial in many situations. The geographical sites with better conditions can be more precisely mapped, in order to optimize the structure of future global quantum networks [65–69]. Through the use of this model, the data from meteorological predictions can directly be linked to the key rate achievable by the QKD link, allowing more accurate statistical studies of the number of operative days per year. The characterization of the transmittance of the channel has then be used to evaluate the performance of the link in terms of achievable secret key rates. We focused on two implementations of the BB-84 cryptographic protocol, using single photons and weak coherent pulses. The noise expected in interesting real-life scenarios, during night-time and day-time, has been modeled and taken into account. We also pointed out the importance of finite-key effects, which can be very detrimental due to the short duration of the link between ground station and satellite. The simulations confirm that long-distance quantum communication can be achieved not only using medium-sized satellites, like the Chinese Micius, but also nano-satellites, allowing to considerably cut the cost of a space-based global quantum network. Ultimately, such links are expected to be integrated with a repeater-based quantum network on the ground, to complement it and enhance the key rate when long distances need to be bridged. The analysis of such a configuration and the optimization of its topology and structure are still under study and represent a crucial milestone towards the realization of the dreamt quantum internet.

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 675662.

In this section we summarize the analysis of atmospheric optical channels proposed in [21, 22]. We will discuss the background, show the main steps of the derivation and recap some results, that will be used as starting point for the simulations described in section 3.

We start from the introduction to the problem given in section 1. The solution of the paraxial wave equation, with phase-approximation using the Huygens–Kirchhoff method [70], can be written in the following way

$$\begin{eqnarray}&&u({\boldsymbol{\rho }},L)={\displaystyle \int }_{{\mathbb{R}}}{{\rm{d}}}^{2}{\boldsymbol{\rho }}^{\prime} {u}_{0}({\boldsymbol{\rho }}^{\prime} ){G}_{0}({\boldsymbol{\rho }},{\boldsymbol{\rho }}^{\prime} ;L,0)\exp [{\rm{i}}{S}({\boldsymbol{\rho }},{\boldsymbol{\rho }}^{\prime} ;z,z^{\prime} )].\end{eqnarray} \tag{ A.1 }$$

Since the losses due to back-scattering and absorption cannot be included in the paraxial approximation of the Helmholtz equation, we treat them phenomenologically multiplying the beam envelope $u({\boldsymbol{\rho }},L)$ by $\sqrt{{\chi }_{\mathrm{ext}}}$. The extinction factor ${\chi }_{\mathrm{ext}}\in [0,1]$ accounts for absorption and back-scattering losses and can be considered as a non-fluctuating quantity (see [22]). In equation (A.1) ${u}_{0}({\boldsymbol{\rho }}^{\prime} )$ is the Gaussian envelope at the transmitter plane (z = 0, z is the longitudinal coordinate)

$$\begin{eqnarray}&&{u}_{0}({\boldsymbol{\rho }})=\sqrt{\displaystyle \frac{2}{\pi {W}_{0}^{2}}}\ \exp \left[-\displaystyle \frac{1}{{W}_{0}^{2}}| {\boldsymbol{\rho }}{| }^{2}-\displaystyle \frac{{\rm{i}}{k}}{2F}| {\boldsymbol{\rho }}{| }^{2}\right],\end{eqnarray} \tag{ A.2 }$$

with W₀ the beam spot radius at the transmitter, k the optical wavenumber and F the focal length of the beam. G₀ is a Gaussian integral kernel

$$\begin{eqnarray}&&{G}_{0}({\boldsymbol{\rho }},{\boldsymbol{\rho }}^{\prime} :z,z^{\prime} )=\displaystyle \frac{k}{2\pi {\rm{i}}(z-z^{\prime} )}\ \exp \left[\displaystyle \frac{{\rm{i}}{k}| {\boldsymbol{\rho }}-{\boldsymbol{\rho }}^{\prime} {| }^{2}}{2(z-z^{\prime} )}\right]\end{eqnarray} \tag{ A.3 }$$

while S contains all the atmospheric effects

$$\begin{eqnarray}&&S({\boldsymbol{\rho }},{\boldsymbol{\rho }}^{\prime} ;z,z^{\prime} )=\displaystyle \frac{k}{2}{\int }_{z^{\prime} }^{z}{\rm{d}}\zeta \ \delta \varepsilon \left({\boldsymbol{\rho }}\displaystyle \frac{\zeta -z^{\prime} }{z-z^{\prime} }+{\boldsymbol{\rho }}^{\prime} \displaystyle \frac{z-\zeta }{z-z^{\prime} },\zeta \right).\end{eqnarray} \tag{ A.4 }$$

Here $S({\boldsymbol{\rho }},{\boldsymbol{\rho }}^{\prime} ;z,z^{\prime} )$ gives the phase contribution due to inhomogeneities of the relative permittivity of the air $\delta \epsilon ({\boldsymbol{\rho }}^{\prime\prime} ,\xi )$ from $z^{\prime}$ to z. Note that δ ε can be separated in two contributions, related to turbulence and scattering

$$\begin{eqnarray}&&\delta \varepsilon =\delta {\varepsilon }_{\mathrm{turb}}+\delta {\varepsilon }_{\mathrm{scat}}.\end{eqnarray} \tag{ A.5 }$$

Assuming the two contributions to be statistically independent, the same factorization holds for the permittivity fluctuation spectrum

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\varepsilon }({\bf{K}})={{\rm{\Phi }}}_{\varepsilon }^{\mathrm{turb}}({\bf{K}})+{{\rm{\Phi }}}_{\varepsilon }^{\mathrm{scat}}({\bf{K}}),\end{eqnarray} \tag{ A.6 }$$

defined as the Fourier transform of the correlation function of $\delta \varepsilon ({\bf{r}})$

$$\begin{eqnarray}&&\langle \delta \varepsilon ({{\bf{r}}}_{1})\delta \varepsilon ({{\bf{r}}}_{2})\rangle =\int {{\rm{d}}}^{3}{\bf{K}}\ {{\rm{\Phi }}}_{\varepsilon }({\bf{K}})\ \exp [{\rm{i}}{\bf{K}}\cdot ({{\bf{r}}}_{1}-{{\bf{r}}}_{2})].\end{eqnarray} \tag{ A.7 }$$

In the previous equations K denotes the momentum and is a 3-dimensional vector. The Markov approximation (that in our case corresponds to assume delta-correlation in the z direction) simplify this expression [25, 27]

$$\begin{eqnarray}&&\langle \delta \varepsilon ({{\bf{r}}}_{1})\delta \varepsilon ({{\bf{r}}}_{2})\rangle =2\pi \delta ({z}_{1}-{z}_{2})\int {{\rm{d}}}^{2}{\bf{k}}\ {{\rm{\Phi }}}_{\varepsilon }({\bf{k}})\ \exp [{\rm{i}}{\bf{k}}\cdot ({{\boldsymbol{\rho }}}_{1}-{{\boldsymbol{\rho }}}_{2})],\end{eqnarray} \tag{ A.8 }$$

where k represents the momentum in the plane transverse to the propagation direction. ρ₁ and ρ₂ are the components of the vectors r₁ and r₂ in the transversal plane, while δ(z) is the Dirac-delta. The Kolmogorov model allows us to write the turbulence-related part of the relative permittivity fluctuation spectrum as [24–27]

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\varepsilon }^{\mathrm{turb}}({\bf{k}})=0.132\ {C}_{n}^{2}\ | {\bf{k}}{| }^{-\tfrac{11}{3}}.\end{eqnarray} \tag{ A.9 }$$

The refractive index structure constant C_n² characterizes the strength of turbulence in the optical domain and is an important parameter of the model. The scattering term in equation (A.6) can be approximated as a Gaussian function [28–32]

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\varepsilon }^{\mathrm{scat}}({\bf{k}})=\displaystyle \frac{{n}_{0}{\zeta }_{0}^{4}}{8\pi {k}^{2}}\ \exp \left[-{\zeta }_{0}^{2}| {\bf{k}}{| }^{2}\right],\end{eqnarray} \tag{ A.10 }$$

with ζ₀ correlation length of the fluctuations due to scattering particles. Here n₀ is the mean number of scatterers per unit volume and represents the main parameter in describing the strength of the scattering contribution.

Now we want to use these ingredients to calculate the probability distribution of the parameters in the elliptic beam approximation, introduced in section 1

$$\begin{eqnarray}&&{\bf{v}}=({x}_{0},{y}_{0},{W}_{1},{W}_{2},{\varphi }_{0}).\end{eqnarray} \tag{ A.11 }$$

First of all we define normalized variables from the ellipse semi-axes

$$\begin{eqnarray}&&{{\rm{\Theta }}}_{i}=\mathrm{ln}\left(\displaystyle \frac{{W}_{i}^{2}}{{W}_{0}^{2}}\right)\quad i=1,2,\end{eqnarray} \tag{ A.12 }$$

where W₀ is the beam spot radius at the transmitter. Now we assume that, in the case of uniform turbulence and scatterers density, the probability distribution of ${x}_{0},{y}_{0},{{\rm{\Theta }}}_{1},{{\rm{\Theta }}}_{2}$ is Gaussian, while the angle of orientation φ₀ is uniformly distributed in [0, π/2] (see appendix of [22] for details). The mean value and the variance of these distributions can be analytically computed. We recall the main steps of the derivation in the following paragraphs.

Starting from the beam centroid position (x₀, y₀), we can choose the reference frame such that $\langle {x}_{0}\rangle =\langle {y}_{0}\rangle =0$ and [25, 71]

$$\begin{eqnarray}&&\langle {x}_{0}^{2}\rangle =\langle {y}_{0}^{2}\rangle ={\int }_{{{\mathbb{R}}}^{4}}{{\rm{d}}}^{2}{{\boldsymbol{\rho }}}_{1}{{\rm{d}}}^{2}{{\boldsymbol{\rho }}}_{2}{x}_{1}{x}_{2}{{\rm{\Gamma }}}_{4}({{\boldsymbol{\rho }}}_{1},{{\boldsymbol{\rho }}}_{2};L).\end{eqnarray} \tag{ A.13 }$$

Here ${{\rm{\Gamma }}}_{4}({{\boldsymbol{\rho }}}_{1},{{\boldsymbol{\rho }}}_{2};z)=\langle {u}^{* }({{\boldsymbol{\rho }}}_{1},z)u({{\boldsymbol{\rho }}}_{1},z){u}^{* }({{\boldsymbol{\rho }}}_{2},z)u({{\boldsymbol{\rho }}}_{1},z)\rangle$ is the fourth-order field correlation function.

The means and covariances of the squared ellipse semi-axes W_i² have the following form (see appendix of [21] for details)

$$\begin{eqnarray}&&\langle {W}_{1/2}^{2}\rangle =4\left[{\displaystyle \int }_{{{\mathbb{R}}}^{2}}{{\rm{d}}}^{2}{\boldsymbol{\rho }}\ {x}^{2}\ {{\rm{\Gamma }}}_{2}({\boldsymbol{\rho }};L)-\langle {x}_{0}^{2}\rangle \right],\end{eqnarray} \tag{ A.14 }$$

$$\begin{eqnarray}\begin{array}{rcl}\langle {\rm{\Delta }}{W}_{i}^{2}{\rm{\Delta }}{W}_{j}^{2}\rangle & = & -8\left\{2{\left({\displaystyle \int }_{{{\mathbb{R}}}^{2}}{{\rm{d}}}^{2}{\boldsymbol{\rho }}{x}^{2}{{\rm{\Gamma }}}_{2}({\boldsymbol{\rho }};L)\right)}^{2}\right.\\ & & -{\displaystyle \int }_{{{\mathbb{R}}}^{4}}{{\rm{d}}}^{2}{{\boldsymbol{\rho }}}_{1}\ {{\rm{d}}}^{2}{{\boldsymbol{\rho }}}_{2}\ [{x}_{1}^{2}{x}_{2}^{2}(4{\delta }_{{ij}}-1)-{x}_{1}^{2}{y}_{2}^{2}(4{\delta }_{{ij}}-3)]\\ & & \left.\times {{\rm{\Gamma }}}_{4}({{\boldsymbol{\rho }}}_{1},{{\boldsymbol{\rho }}}_{2};L)\right\}-16[4{\delta }_{{ij}}-1]\langle {x}_{0}^{2}{\rangle }^{2},\end{array}\end{eqnarray} \tag{ A.15 }$$

where the second-order field correlation function ${{\rm{\Gamma }}}_{2}({\boldsymbol{\rho }};z)=\langle {u}^{* }({\boldsymbol{\rho }},z)u({\boldsymbol{\rho }},z)\rangle$ has been used.

The next step is the calculation of the field correlation function, for which we use the expression of the beam envelope given in equation (A.1). We report the calculations only for ${{\rm{\Gamma }}}_{2}({\boldsymbol{\rho }};L)$, the equivalent but more cumbersome expressions for ${{\rm{\Gamma }}}_{4}({{\boldsymbol{\rho }}}_{1},{{\boldsymbol{\rho }}}_{2};L)$ can be found in [22], appendix B. Substituting equation (A.1) in the definition of ${{\rm{\Gamma }}}_{2}({\boldsymbol{\rho }};z)$ yields

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Gamma }}}_{2}({\boldsymbol{\rho }};L) & = & {\displaystyle \int }_{{{\mathbb{R}}}^{4}}{{\rm{d}}}^{2}\ {{\boldsymbol{\rho }}}_{1}^{{\prime} }{{\rm{d}}}^{2}{{\boldsymbol{\rho }}}_{2}^{{\prime} }{u}_{0}({{\boldsymbol{\rho }}}_{1}^{{\prime} }){u}_{0}^{* }({{\boldsymbol{\rho }}}_{2}^{{\prime} }){G}_{0}({\boldsymbol{\rho }},{{\boldsymbol{\rho }}}_{1}^{{\prime} };L,0)\\ & & \times {G}_{0}^{* }({\boldsymbol{\rho }},{{\boldsymbol{\rho }}}_{2}^{{\prime} };L,0)\ \exp \left[-\displaystyle \frac{1}{2}{{ \mathcal D }}_{S}(0,{{\boldsymbol{\rho }}}_{1}^{{\prime} }-{{\boldsymbol{\rho }}}_{2}^{{\prime} })\right],\end{array}\end{eqnarray} \tag{ A.16 }$$

with the last term embodying the phase fluctuations due to the atmosphere (remember the definition of $S({\boldsymbol{\rho }},{\boldsymbol{\rho }}^{\prime} ;z,z^{\prime} )$ in equation (A.4))

$$\begin{eqnarray}&&{{ \mathcal D }}_{S}({{\boldsymbol{\rho }}}_{k}-{{\boldsymbol{\rho }}}_{l},{{\boldsymbol{\rho }}}_{k}^{{\prime} }-{{\boldsymbol{\rho }}}_{l}^{{\prime} })=\left\langle {\left[S({{\boldsymbol{\rho }}}_{k},{{\boldsymbol{\rho }}}_{k}^{{\prime} };z,z^{\prime} )-S({{\boldsymbol{\rho }}}_{l},{{\boldsymbol{\rho }}}_{l}^{{\prime} };z,z^{\prime} )\right]}^{2}\right\rangle .\end{eqnarray} \tag{ A.17 }$$

Substituting equation (A.4) and exploiting again the Markov approximation, the factorization in equations (A.5) and (A.6) can be carried over

$$\begin{eqnarray}&&{{ \mathcal D }}_{S}={{ \mathcal D }}_{S}^{\mathrm{turb}}+{{ \mathcal D }}_{S}^{\mathrm{scat}}.\end{eqnarray} \tag{ A.18 }$$

We can now introduce the models for the permittivity fluctuations spectrum related to turbulence (equation (A.9)) and scatterers (equation (A.10)), obtaining

$$\begin{eqnarray}&&{{ \mathcal D }}_{S}^{\mathrm{turb}}({\boldsymbol{\rho }},{\boldsymbol{\rho }}^{\prime} )=2.95\ {k}^{2}\ L{\int }_{0}^{1}{\rm{d}}\xi \ {C}_{n}^{2}(\xi )\ | {\boldsymbol{\rho }}\xi +{\boldsymbol{\rho }}^{\prime} (1-\xi ){| }^{\tfrac{5}{3}},\end{eqnarray} \tag{ A.19 }$$

$$\begin{eqnarray}&&{{ \mathcal D }}_{S}^{\mathrm{scat}}({\boldsymbol{\rho }},{\boldsymbol{\rho }}^{\prime} )=\displaystyle \frac{\pi }{8}L{\int }_{0}^{1}{\rm{d}}\xi \ {n}_{0}(\xi )\ | {\boldsymbol{\rho }}\xi +{\boldsymbol{\rho }}^{\prime} (1-\xi ){| }^{2},\end{eqnarray} \tag{ A.20 }$$

where we introduced the rescaled longitudinal coordinate $\xi \in [0,1]$, where ξ = 1 corresponds to z = L. We allowed for a dependence on longitudinal coordinate in ${C}_{n}^{2}(\xi )$ and n₀(ξ) for later use. We recall the definition of the so-called Rytov parameter ${\sigma }_{R}^{2}=1.23\ {C}_{n}^{2}\ {k}^{\tfrac{7}{6}}{L}^{\tfrac{11}{6}}$. Substituting in equation (A.16) the definition of the Gaussian envelope ${u}_{0}({\boldsymbol{\rho }})$ (equation (A.2)) and the integral kernel ${G}_{0}({\boldsymbol{\rho }},{\boldsymbol{\rho }}^{\prime} :L,0)$ (equation (A.3)), the second-order field correlation function reads

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Gamma }}}_{2}({\boldsymbol{\rho }};L) & = & \displaystyle \frac{{{\rm{\Omega }}}^{2}}{{\pi }^{2}{W}_{0}^{4}}{\displaystyle \int }_{{{\mathbb{R}}}^{2}}{{\rm{d}}}^{2}{\boldsymbol{\rho }}^{\prime} \ {{\rm{e}}}^{-\tfrac{\alpha }{2{W}_{0}^{2}}| {\boldsymbol{\rho }}^{\prime} {| }^{2}-2i\tfrac{{\rm{\Omega }}}{{W}_{0}^{2}}{\boldsymbol{\rho }}\cdot {\boldsymbol{\rho }}^{\prime} }\\ & & \times \exp \left[-\displaystyle \frac{1}{2}{{ \mathcal D }}_{S}^{{\rm{turb}}}(0,{\boldsymbol{\rho }}^{\prime} )\right]\exp \left[-\displaystyle \frac{1}{2}{{ \mathcal D }}_{S}^{{\rm{scat}}}(0,{\boldsymbol{\rho }}^{\prime} )\right].\end{array}\end{eqnarray} \tag{ A.21 }$$

Here $\alpha =1+{{\rm{\Omega }}}^{2}{\left(1-\tfrac{L}{F}\right)}^{2}$ with the Fresnel number defined as ${\rm{\Omega }}=\tfrac{{{kW}}_{0}^{2}}{2L}$.

Integrating equation (A.21) (and the equivalent one for Γ₄) and then equations (A.14), we obtain the first and second moments of the probability distribution of W_i². Then, the moments for the variables Θ_i are easily obtained from equation (A.12).

We consider now the transmittance, defined in equation (1), of an elliptic beam impinging on a circular aperture of radius a. It can be written as

$$\begin{eqnarray}\begin{array}{rcl}\eta ({x}_{0},{y}_{0},{W}_{1},{W}_{2},{\varphi }_{0}) & = & \displaystyle \frac{2\ {\chi }_{\mathrm{ext}}}{\pi {W}_{1}{W}_{2}}{\displaystyle \int }_{0}^{a}{\rm{d}}\rho {\displaystyle \int }_{0}^{2\pi }{\rm{d}}\theta \ {{\rm{e}}}^{-2{A}_{1}{\left(\rho \cos \theta -{\rho }_{0}\right)}^{2}}\\ & & \times \,{{\rm{e}}}^{-2{A}_{2}{\rho }^{2}{\sin }^{2}\theta }{{\rm{e}}}^{-2{A}_{3}(\rho \cos \theta -{\rho }_{0})r\sin \theta },\end{array}\end{eqnarray} \tag{ A.22 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}{A}_{1} & = & \left(\displaystyle \frac{{\cos }^{2}({\varphi }_{0}-{\theta }_{0})}{{W}_{1}^{2}}+\displaystyle \frac{{\sin }^{2}({\varphi }_{0}-{\theta }_{0})}{{W}_{2}^{2}}\right)\\ {A}_{2} & = & \left(\displaystyle \frac{{\sin }^{2}({\varphi }_{0}-{\theta }_{0})}{{W}_{1}^{2}}+\displaystyle \frac{{\cos }^{2}({\varphi }_{0}-{\theta }_{0})}{{W}_{2}^{2}}\right)\\ {A}_{3} & = & \left(\displaystyle \frac{1}{{W}_{1}^{2}}-\displaystyle \frac{1}{{W}_{2}^{2}}\right)\ \sin \ 2({\varphi }_{0}-{\theta }_{0}).\end{array}\end{eqnarray} \tag{ A.23 }$$

In the previous equations (ρ, θ) are the integration variables in the area of the circular aperture, while $({x}_{0},{y}_{0})=({\rho }_{0}\ \cos {\theta }_{0},{\rho }_{0}\ \sin {\theta }_{0})$ is the beam-centroid position.

The PDT is then easily reconstructed. Extract at random M 5-tuples of values for (x₀, y₀, Θ₁, Θ₂, φ₀), according to the correct probability distribution. Compute first the values of the ellipse semi-axes W_i from Θ_i and then the value of the transmittance for every tuple. Collect the statistics in an histogram and compute statistical estimators (e.g. the median). Two examples of the simulated PDT are shown in section 3 of the main text (figures 3 and 4).

In this section we are going to apply the model described in the previous section to a satellite-based link, as described in section 3 of the main text. We will discuss some details about the calculations involved and show an example of how to proceed with the integration of the expressions in appendix A. In particular, we will focus on the first term of the quantity $\langle {W}_{1/2}^{2}\rangle$ defined in equation equation (A.14), which only contains the second order correlation function ${{\rm{\Gamma }}}_{2}({\boldsymbol{\rho }};L)$. The computations involving the integration of the fourth order correlation function ${{\rm{\Gamma }}}_{4}({{\boldsymbol{\rho }}}_{1},{{\boldsymbol{\rho }}}_{2};L)$ are much more cumbersome and will not be reported here.

Inserting equation (A.21) in the first term of equation (A.14) we obtain the following integration, where all the quantities are defined in the previous section appendix A

$$\begin{eqnarray}\begin{array}{rcl}{\displaystyle \int }_{{{\mathbb{R}}}^{2}}{{\rm{d}}}^{2}{\boldsymbol{\rho }}\ {x}^{2}\ {{\rm{\Gamma }}}_{2}({\boldsymbol{\rho }};L) & = & \displaystyle \frac{{{\rm{\Omega }}}^{2}}{{\pi }^{2}{W}_{0}^{4}}{\displaystyle \int }_{{{\mathbb{R}}}^{4}}{{\rm{d}}}^{2}{\boldsymbol{\rho }}\ {{\rm{d}}}^{2}{\boldsymbol{\rho }}^{\prime} \ {x}^{2}\ {{\rm{e}}}^{-\tfrac{\alpha }{2{W}_{0}^{2}}| {\boldsymbol{\rho }}^{\prime} {| }^{2}-2{\rm{i}}\tfrac{{\rm{\Omega }}}{{W}_{0}^{2}}{\boldsymbol{\rho }}\cdot {\boldsymbol{\rho }}^{\prime} }\\ & & \times \exp \left[-\displaystyle \frac{1}{2}{{ \mathcal D }}_{S}^{\mathrm{turb}}(0,{\boldsymbol{\rho }}^{\prime} )\right]\exp \left[-\displaystyle \frac{1}{2}{{ \mathcal D }}_{S}^{\mathrm{scat}}(0,{\boldsymbol{\rho }}^{\prime} )\right].\end{array}\end{eqnarray} \tag{ B.1 }$$

We assume that the beam is focused (F = L) so that α = 1. First of all we can compute the terms ${{ \mathcal D }}_{S}^{\mathrm{turb}}(0,{\boldsymbol{\rho }}^{\prime} )$ and ${{ \mathcal D }}_{S}^{\mathrm{scat}}(0,{\boldsymbol{\rho }}^{\prime} )$ defined in equations (A.19) and (A.20)

$$\begin{eqnarray}&&{{ \mathcal D }}_{S}^{\mathrm{turb}}(0,{\boldsymbol{\rho }}^{\prime} )=2.95\ {k}^{2}\ L{\int }_{0}^{1}{\rm{d}}\xi \ {C}_{n}^{2}(\xi )\ | {\boldsymbol{\rho }}^{\prime} (1-\xi ){| }^{\tfrac{5}{3}},\end{eqnarray} \tag{ B.2 }$$

$$\begin{eqnarray}&&{{ \mathcal D }}_{S}^{\mathrm{scat}}(0,{\boldsymbol{\rho }}^{\prime} )=\displaystyle \frac{\pi }{8}L{\int }_{0}^{1}{\rm{d}}\xi \ {n}_{0}(\xi )\ | {\boldsymbol{\rho }}^{\prime} (1-\xi ){| }^{2}.\end{eqnarray} \tag{ B.3 }$$

In the rescaled longitudinal coordinate ξ, the conditions in equation (3) in the main text become

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{Down}-\mathrm{links}\qquad {C}_{n}^{2}(\xi ) & = & {C}_{n}^{2}\ {\rm{\Theta }}(\xi -(1-h/L))\\ {n}_{0}(\xi ) & = & {n}_{0}\ {\rm{\Theta }}(\xi -(1-h/L))\\ \mathrm{Up}-\mathrm{links}\qquad \quad {C}_{n}^{2}(\xi ) & = & {C}_{n}^{2}\ {\rm{\Theta }}(h/L-\xi )\\ {n}_{0}(\xi ) & = & {n}_{0}\ {\rm{\Theta }}(h/L-\xi ).\end{array}\end{eqnarray} \tag{ B.4 }$$

Inserting equation (B.4) (we consider down-link in this example) into equations (B.2) and (B.3) we can solve the integration and obtain

$$\begin{eqnarray}&&{{ \mathcal D }}_{S}^{\mathrm{turb}}(0,{\boldsymbol{\rho }}^{\prime} )=2.95\ {k}^{2}\ L\ | {\boldsymbol{\rho }}^{\prime} {| }^{\tfrac{5}{3}}{C}_{n}^{2}{\displaystyle \int }_{1-h/L}^{1}{\rm{d}}\xi \ \ {\left(1-\xi \right)}^{\tfrac{5}{3}}=2.4\ {\sigma }_{R}^{2}\ {k}^{5/6}\ {L}^{-5/6}\ | {\boldsymbol{\rho }}^{\prime} {| }^{\tfrac{5}{3}}\left[\displaystyle \frac{3}{8}{\left(\displaystyle \frac{h}{L}\right)}^{8/3}\right],\end{eqnarray} \tag{ B.5 }$$

$$\begin{eqnarray}&&{{ \mathcal D }}_{S}^{\mathrm{scat}}(0,{\boldsymbol{\rho }}^{\prime} )=\displaystyle \frac{\pi }{8}{{Ln}}_{0}| {\boldsymbol{\rho }}^{\prime} {| }^{2}{\displaystyle \int }_{1-h/L}^{1}{\rm{d}}\xi \ \ | {\boldsymbol{\rho }}^{\prime} (1-\xi ){| }^{2}=\displaystyle \frac{\pi }{8}{{Ln}}_{0}| {\boldsymbol{\rho }}^{\prime} {| }^{2}\left[\displaystyle \frac{1}{3}{\left(\displaystyle \frac{h}{L}\right)}^{3}\right],\end{eqnarray} \tag{ B.6 }$$

where ${\sigma }_{R}^{2}$ has been defined in section 3. From this passage we clearly see where the dependency on $\tfrac{h}{L}$ in equations (4)–(9) originates from. When we introduce equations (B.5) and (B.6) in (B.1), we recognize that it only contains Gaussian integrals of the form

$$\begin{eqnarray}&&{\int }_{-\infty }^{\infty }{\rm{d}}{x}\ {x}^{c}\ \exp [a\ {x}^{2}+i\ b\ x],\end{eqnarray} \tag{ B.7 }$$

with c = {0, 2} and that can be readily solved. The only exception is the turbulence term, which contains $| {\boldsymbol{\rho }}^{\prime} {| }^{\tfrac{5}{3}}$. We can simplify the computation introducing the approximation [22, 71] $| {\boldsymbol{\rho }}^{\prime} /{W}_{0}{| }^{\tfrac{5}{3}}\simeq | {\boldsymbol{\rho }}^{\prime} /{W}_{0}{| }^{2}$. Then, one just has to solve the multiple Gaussian integrals and insert it in equation (A.14) to obtain the value of $\langle {W}_{1/2}^{2}\rangle$ for down-links, as in equation (8). Similar techniques can be used to compute all the other moments of the beam variables. Equations (4), (5), (7) and (8) have been computed specifically for the problem at hand, the non-uniform link described at the beginning of section 3. Equations (6) and (9), on the other hand, have been deduced from the equivalent results obtained in [22] for a uniform link. We see that in equations (4), (5), (7) and (8) the corrections due to the non-uniformity of the link (of the form $A{\left(h/L\right)}^{\beta }$, where A is a constant and β = {1, 8/3, 3}) act like multiplicative factors on the parameters ${\sigma }_{R}^{2}$ and n₀. So, we started from the calculation of the quantity $\langle {\rm{\Delta }}{W}_{i}^{2}{\rm{\Delta }}{W}_{j}^{2}\rangle$ in [22] and attached the multiplicative corrections found above, in order to obtain equations (6) and (9). This inconsistency should not be considered too detrimental regarding the reliability of the model. We checked through the simulation that the mean value and the shape of the PDT are not very sensitive to variations of the value of the quantities in equations (6) and (9), as the interplay between beam wandering (equations (4) and (7)) and beam spreading (equations (5) and (8)) is much more significant in this context. Finally, we point out that the computation of the quantities in equations (4), (5), (7) and (8) have been carried out without the introduction of the weak turbulence approximation used in [21, 22]. For equations (6) and (9), instead, we used the results obtained in [22] in the weak turbulence regime, which we verified to be still valid in the case of satellite-based links.

In a free-space link, environmental photons are usually the most important source of noise. In this section we summarize the analysis of [11, 23] regarding the amount of environmental photons that hit the detector for down-links and up-links, that we use to calculate the QBER. We suppose that an accurate time synchronization had been operated between sender and receiver, in order to tag the photons and perform a time filtering on the incoming signal. On top of that, wavelength filtering is applied to further reduce the amount of detected noisy photons.

For up-links, we only consider the case of night-time operation. If the ground station site has a low level of light pollution, the biggest fraction of environmental photons comes from the sunlight reflected first by the Moon and then by the Earth [11]

$$\begin{eqnarray}&&{N}_{\mathrm{night}}^{\mathrm{up}}={A}_{E}{A}_{M}{R}_{M}^{2}{a}^{2}\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{fov}}}{{d}_{{EM}}^{2}}{B}_{f}\ {\rm{\Delta }}t\ {H}_{\mathrm{sun}}.\end{eqnarray} \tag{ C.1 }$$

Here A_M and R_M are the albedo and the radius of the Moon, while A_E is the albedo of the Earth and d_EM is the Earth–Moon distance. H_sun is the solar spectral irradiance in $\mathrm{photons}\,{{\rm{s}}}^{-1}\,{\mathrm{nm}}^{-1}\,{{\rm{m}}}^{-2}$ at the wavelength of interest. Ω_fov and a are angular field of view and radius of the receiving telescope. B_f is the width of the spectral filtering and Δt is the detection time-window. We assumed Lambertian diffusion on the Moon and the Earth.

For down-links, the evaluation of the background photons is strongly site-dependent. The power received by the telescope can be expressed as follows [23]

$$\begin{eqnarray}&&{P}_{b}={H}_{b}{{\rm{\Omega }}}_{\mathrm{fov}}\pi {a}^{2}{B}_{f}.\end{eqnarray} \tag{ C.2 }$$

The parameter H_b is the total brightness of the sky background and it depends on the hour of the day and the weather conditions. From equation (C.2) we derive the number of photons per time window

$$\begin{eqnarray}&&{N}^{\mathrm{down}}=\displaystyle \frac{{H}_{b}}{h\nu }{{\rm{\Omega }}}_{\mathrm{fov}}\pi {a}^{2}{B}_{f}\ {\rm{\Delta }}t,\end{eqnarray} \tag{ C.3 }$$

where h is the Planck constant and ν is the frequency of the background photons (after filtering). Typical values of the brightness of the sky are ${H}_{b}={10}^{-3}\ W\ {{\rm{m}}}^{-2}\ \mathrm{sr}\ \mu {\rm{m}}$ during a full-Moon night and ${H}_{b}=1\ W\ {{\rm{m}}}^{-2}\ \mathrm{sr}\ \mu {\rm{m}}$ for a clear sky in day-time. This analysis assumes that neither the Moon during the night nor the Sun during the day are included in the field of view of the collecting aperture.

The QBER is computed assuming the noisy photons to be completely unpolarized

$$\begin{eqnarray}&&{\rm{QBER}}={Q}_{0}+\displaystyle \frac{1}{2}\displaystyle \frac{{N}_{\mathrm{noise}}}{{N}_{\mathrm{noise}}+{N}_{\mathrm{sig}}}.\end{eqnarray} \tag{ C.4 }$$

Here Q₀ corresponds to the error rate associated with depolarization in the encoding degree of freedom or imperfection of the preparation or detection stage leading to incorrect state discrimination. We chose a conservative value of Q₀ = 2%. N_noise and N_sig are, respectively, the number of photons per time window associated to noise and signal. As expected, the number of collected environmental photons are proportional to the area of the receiving aperture, but so is the intensity of the signal. To reduce the noise and at the same time raise the signal to noise ratio, we can act on Ω_fov, B_f and Δt. Reducing the field of view involves a better pointing and tracking system, while a very good time synchronization allows the use of short time windows.

We report here the expression of the secret key rates we used in the performance study of section 4. The set-up is the usual one for QKD: two parties, A and B, are connected through a completely insecure quantum channel and an authenticated classical channel. After many uses of the links, their goal is to share an identical key, which is secret regardless of the attack strategy that an hypothetical eavesdropper could implement. For the single-photon implementation of the BB-84 protocol (using, e.g. polarization encoding), party A sends qubits in the basis $X=\{\left|0\right\rangle ,\left|1\right\rangle \}$ or $Z=\{\left|+\right\rangle ,\left|-\right\rangle \}$ at random, with $\left|\pm \right\rangle =(\left|0\right\rangle \pm \left|1\right\rangle )/\sqrt{2}$. B measures the received qubits in the bases X or Z, at random. The results of [53] state that a secret key of length l can be shared, if

$$\begin{eqnarray}\begin{array}{rcl}l & \leqslant & n(q-{h}_{2}({Q}_{\mathrm{tol}}+\mu ))-{\mathrm{leak}}_{{EC}}-\alpha ({\varepsilon }_{\sec },{\varepsilon }_{\mathrm{cor}})\\ \alpha ({\varepsilon }_{\sec },{\varepsilon }_{\mathrm{cor}}) & = & {\mathrm{log}}_{2}\displaystyle \frac{2}{{\varepsilon }_{\sec }^{2}{\varepsilon }_{\mathrm{cor}}}\ \ \mu =\sqrt{\displaystyle \frac{n+k}{{nk}}\displaystyle \frac{k+1}{k}\mathrm{ln}\displaystyle \frac{2}{{\varepsilon }_{\sec }}},\end{array}\end{eqnarray} \tag{ D.1 }$$

out of n successfully exchanged single photon signals, where the function h₂ denotes the binary entropy. Here q is a parameter describing the preparation quality of the initial states of the signal sent by A. In the qubit case it is connected to the maximum fidelity allowed between states prepared in the X and Z bases. In a perfect implementation of the BB-84 protocol, like the one considered here, the two bases are mutually unbiased, for which the maximum q = 1 is achieved. Q_tol is the channel error tolerance and k is the number of bits of the raw key used for parameter estimation. The achievable key rate is obtained by maximizing over these two parameters. The term leak_EC gives the amount of information in bits that the parties had to exchange during the error correction phase. The desired security and correctness thresholds are specified by the parameters ε_sec and ε_cor.

An alternative protocol based on decoy states [54, 56] is used when the source emits weak coherent pulses instead of real single photons. We follow the analysis of [55, 57], where two decoy states are used. The bases used for the encoding are X and Z as in the single photon implementation. A secret key of length l can be extracted, with

$$\begin{eqnarray}&&l\leqslant {s}_{X,0}+{s}_{X,1}(1-{h}_{2}({\phi }_{X}))-{\mathrm{leak}}_{{EC}}-6\ {\mathrm{log}}_{2}\displaystyle \frac{21}{{\varepsilon }_{\sec }}-{\mathrm{log}}_{2}\displaystyle \frac{2}{{\varepsilon }_{\mathrm{cor}}}.\end{eqnarray} \tag{ D.2 }$$

s_X,0 and s_X,1 represent the number of bits in the raw key generated by vacuum events and single photon events, respectively. ϕ_X instead is the phase error rate measured in the channel during parameter estimation. The subscript X means that these estimations are valid for the events in which both A and B chose the basis X and they include the corrections due to finite key effects (for the actual expressions we refer to [57]). In this case the maximization is over the portion of signals used for parameter estimation, the intensity of the signal and decoy states and the probability of sending different intensities.

In both cases, the key rates are obtained taking the ratio between the length in bit of the final secret key l and total number of signals sent n.

In this section we show the values of the parameters utilized throughout the paper and we discuss about their pertinence. They are reported in tables E1, E2 and E3, together with a brief explanatory description, where necessary. More detailed explanations about particular parameters are in the remainder of this section.

Table E3.  Parameters related to stray photons and environmental light.

Parameter	Value	Brief description
Sky brightness Hb	$1.5\,\times \,{10}^{-6}\ {\rm{W}}\ {{\rm{m}}}^{-2}\ {\mathrm{sr}}^{-1}\,{\mathrm{nm}}^{-1}$	Night, clear sky, full Moon [23]
Sky brightness Hb	$1.5\,\times \,{10}^{-3}\ {\rm{W}}\ {{\rm{m}}}^{-2}\ {\mathrm{sr}}^{-1}\,{\mathrm{nm}}^{-1}$	Day, clear sky [23]
Field of view Ωfov	${\left(100\times {10}^{-6}\right)}^{2}\,\mathrm{sr}$	Night-time down-link
Field of view Ωfov	${\left(10\times {10}^{-6}\right)}^{2}\,\mathrm{sr}$	Day-time down-link
Field of view Ωfov	${\left(30\times {10}^{-6}\right)}^{2}\,\mathrm{sr}$	Night-time up-link
Time-window Δt	1 ns	Night- and day-time
Spectral filter width Bf	1 nm	Night-time down-link
Spectral filter width Bf	0.2 nm	Day-time down-link
Spectral filter width Bf	1 nm	Night-time up-link
Hsun	$4.610\,\times \,{10}^{18}\ \mathrm{phot}\ {{\rm{s}}}^{-1}\,{\mathrm{nm}}^{-1}\,{{\rm{m}}}^{-2}$	Solar spectral irradiance
Ae	0.300	Earth'salbedo
Am	0.136	Moon's albedo
RM	1.737 × 106 m	Moon's radius
dEM	3.600 × 108 m	Earth–Moon distance

The parameters C_n², n₀ and h should in general be fixed by fitting the experimental data. However, in order to have a predictive model, we want to estimate these parameters in a reasonable way. First of all, in order to estimate the effective thickness of the atmosphere, we start from the variation of density of the air as a function of the altitude. We chose h = 20 km, as a layer around the Earth with this thickness contains on average 95% of the total mass of the atmosphere. As already stated in the main text, some models for the altitude-dependence of the refractive index structure constant C_n² are available in the literature [34–37]. The widely used parametric fit due to Hufnagel and Valley [34, 35] reliably replicates the behaviour of C_n² in mid-latitude climate

$$\begin{eqnarray}&&{C}_{n}^{2}(z)=5.94\,\times \,{10}^{-53}{\left(\displaystyle \frac{v}{27}\right)}^{2}{z}^{10}\ \exp [-z/1000]+2.7\,\times \,{10}^{-16}\ \exp [z/1500]+A\ \exp [z/100].\end{eqnarray} \tag{ E.1 }$$

Here z is the altitude coordinate, v is a parameter related to high-altitude wind speed and A describes the relative strength of the turbulence near the ground level. Typical values are $A=1.7\,\times \,{10}^{-14}\ {{\rm{m}}}^{-2/3}$ and v = 21 m s⁻¹, although v = 57 m s⁻¹ is sometimes used for stronger wind conditions. The value of C_n² inside the atmosphere in our model is estimated by the integral average of this function in $[0,\infty ]$, rescaled by the fixed thickness h

$$\begin{eqnarray}&&{C}_{n}^{2}=\displaystyle \frac{1}{h}{\int }_{0}^{\infty }{C}_{n}^{2}(z)\ {\rm{d}}{z}.\end{eqnarray} \tag{ E.2 }$$

The parameter v is kept fixed to the recommended value of 21 m s⁻¹. A is chosen to match the values of ${C}_{n}^{2}(0)$ measured in [22], ${A}_{n}=1.10\,\times \,{10}^{-14}\ {{\rm{m}}}^{-2/3}$ at night and ${A}_{d}=2.75\,\times \,{10}^{-14}\ {{\rm{m}}}^{-2/3}$ during the day. Through equation (E.1), the first corresponds to ${C}_{n}^{2}=1.12\,\times \,{10}^{-16}\ {{\rm{m}}}^{-2/3}$ and the latter to ${C}_{n}^{2}=1.64\,\times \,{10}^{-16}\ {{\rm{m}}}^{-2/3}$.

The scattering particles described by the density n₀ mainly consist of water droplets, so, in order to estimate the value of n₀, we start from the profile of the water vapour content in the atmosphere. The absolute humidity vertical profile τ(z) in the range [0, 10 km] can be written as a double exponential [72, 73]

$$\begin{eqnarray}\begin{array}{rcl}\tau (z) & = & \tau (0)\ \exp [-z/{H}_{1}]\qquad \qquad \qquad \ \ \mathrm{for}\ 0\leqslant z\leqslant 5\ \mathrm{km}\\ & = & \tau ({H}_{1})\ \exp [-(z-5\ \mathrm{km})/{H}_{2}]\qquad \mathrm{for}\ 5\ \mathrm{km}\leqslant z\leqslant 10\ \mathrm{km}\end{array}\end{eqnarray} \tag{ E.3 }$$

with the two scale heights H₁ and H₂. The contribution of the region with z > 10 km is rather low and we neglect it here. The parameters H₁ and H₁ can on average vary in the range [1.53, 2.8] and [1.19, 1.82], respectively, depending on the geographical position and the season. We choose in the following the values stated in the US Standard Atmosphere (1962) [74], H₁ = 2.243 and H₂ = 1.414. We obtain a rescaling factor ω in the same way as we did in the previous case

$$\begin{eqnarray}&&\omega =\displaystyle \frac{1}{h\ \tau (0)}{\int }_{0}^{10\ \mathrm{km}}\tau (z)\ {\rm{d}}{z}.\end{eqnarray} \tag{ E.4 }$$

Then, the value of the parameter ${n}_{0}^{* }$ in our case is obtained multiplying by the factor ω the value found in [22] for night- and day-time, ${n}_{0}^{* }=\omega \ {n}_{0}$. For the given values of the scale heights ω ≃ 0.107.

The extinction factor χ_ext(θ) varies as a function of the elevation angle in the following way

$$\begin{eqnarray}&&{\chi }_{\mathrm{ext}}(\theta )=\exp \left[-\beta \ \sec (\theta )\right].\end{eqnarray} \tag{ E.5 }$$

The value of the parameter β reported in table E1 has been chosen to match the amount of extinction used in [10], based on the MODTRAN5 software [38].