Satellite-based links for Quantum Key Distribution: beam effects and weather dependence

TL;DR

This paper develops a detailed model for atmospheric losses affecting satellite-based quantum communication links and analyzes the performance of various QKD implementations under different weather conditions.

Contribution

It provides a comprehensive atmospheric loss model and evaluates QKD performance considering finite-key effects for practical satellite communication scenarios.

Findings

Weather significantly impacts quantum signal transmission.

Finite-key effects influence QKD protocol performance.

The model enables more accurate feasibility assessments for satellite QKD.

Abstract

The establishment of a world-wide quantum communication network relies on the synergistic integration of satellite-based links and fiber-based networks. The first are helpful for long-distance communication, as the photon losses introduced by the optical fibers are too detrimental for lengths greater than about 200 km. This work aims at giving, on the one hand, a comprehensive and fundamental model for the losses suffered by the quantum signals during the propagation along an atmospheric free-space link. On the other hand, a performance analysis of different Quantum Key Distribution (QKD) implementations is performed, including finite-key effects, focusing on different interesting practical scenarios. The specific approach that we chose allows to precisely model the contribution due to different weather conditions, paving the way towards more accurate feasibility studies of satellite-based QKD missions.

Tables

Table 1. Table 1: Parameters related to the optical and technical properties of the link.

Parameter	Value	Brief description
$W_0$	15 cm, 50 cm	down-links, up-links
$W_0$	5 cm, 50 cm	CubeSat down-links, up-links
$a$	50 cm, 15 cm	down-links, up-links
$a$	50 cm, 5 cm	CubeSat down-links, up-links
$\lambda$	785 nm	Wavelength of the signal light
$\beta$	0.7	Parameter in ${\chi }_{\mathrm{ext}}\left(\theta \right)$
$\alpha$	${\mathrm{1.2\; 10}}^{-6}\mathrm{rad}$	Pointing error
${\eta }_{\det }$	0.5	Detector efficiency
$T_{\mathrm{opt}}$	0.8	Transmittance of the optical system

Table 2. Table 2: Parameters related to the atmospheric weather conditions.

Parameter	Value	Brief description
$h$	20 km	Atmosphere thickness
$L$	500 km	Minimum altitude (zenith)
$C_n^2$	${\mathrm{1.12\; 10}}^{-16}{m}^{-2/3}$	Night-time, condition 1
$C_n^2$	${\mathrm{1.64\; 10}}^{-16}{m}^{-2/3}$	Day-time, condition 1
$C_n^2$	${\mathrm{5.50\; 10}}^{-16}{m}^{-2/3}$	Night-time, condition 2
$C_n^2$	${\mathrm{8.00\; 10}}^{-16}{m}^{-2/3}$	Day-time, condition 2
$C_n^2$	${\mathrm{1.10\; 10}}^{-15}{m}^{-2/3}$	Night-time, condition 3
$C_n^2$	${\mathrm{1.60\; 10}}^{-15}{m}^{-2/3}$	Day-time, condition 3
$n_0$	0.61 $m^{-3}$	Night-time, condition 1
$n_0$	0.01 $m^{-3}$	Day-time, condition 1
$n_0$	3.00 $m^{-3}$	Night-time, condition 2
$n_0$	0.05 $m^{-3}$	Day-time, condition 2
$n_0$	6.10 $m^{-3}$	Night-time, condition 3
$n_0$	0.10 $m^{-3}$	Day-time, condition 3

Table 3. Table 3: Parameters related to stray photons and environmental light.

Parameter	Value	Brief description
Sky brightness $H_b$	${\mathrm{1.5\; 10}}^{-6}W{m}^{-2}{\mathrm{sr}}^{-1}{\mathrm{nm}}^{-1}$	Night, clear sky, full Moon [23]
Sky brightness $H_b$	${\mathrm{1.5\; 10}}^{-3}W{m}^{-2}{\mathrm{sr}}^{-1}{\mathrm{nm}}^{-1}$	Day, clear sky [23]
Field of view ${\Omega }_{\mathrm{fov}}$	${\left({\mathrm{100\; 10}}^{-6}\right)}^2$ sr	Night-time down-link
Field of view ${\Omega }_{\mathrm{fov}}$	${\left({\mathrm{10\; 10}}^{-6}\right)}^2$ sr	Day-time down-link
Field of view ${\Omega }_{\mathrm{fov}}$	${\left({\mathrm{10\; 10}}^{-6}\right)}^2$ sr	Night-time up-link
Time-window $\Delta t$	1 ns	Night- and day-time
Spectral filter width $B_f$	1 nm	Night-time down-link
Spectral filter width $B_f$	0.2 nm	Day-time down-link
Spectral filter width $B_f$	1 nm	Night-time up-link
$H_{\mathrm{sun}}$	${\mathrm{4.610\; 10}}^{18}\mathrm{phot}{s}^{-1}{\mathrm{nm}}^{-1}{m}^{-2}$	Solar spectral irradiance
$A_e$	0.300	Earth’s albedo
$A_m$	0.136	Moon’s albedo
$R_M$	${\mathrm{1.737\; 10}}^6$ m	Moon’s radius
$d_{EM}$	${\mathrm{3.600\; 10}}^8$ m	Earth-Moon distance

Equations

$$\eta ={\int }_{|\mathbf{\rho }{|}^2=a^2}{d}^2\mathbf{\rho }|u(\mathbf{\rho },L){|}^2,$$

$$\mathbf{v}=(x_0,y_0,W_1,W_2,{\phi }_0),$$

$$\mathrm{Down}-\mathrm{links}$$

$$\mathrm{Up}-\mathrm{links}$$

$$⟨x_0⟩=⟨y_0⟩=0,⟨x_0^2⟩=⟨y_0^2⟩=0.419{\sigma }_R^2{W}_0^2{\Omega }^{-\frac{7}{6}}\frac{h}{L},$$

$$\langle W_{i}^{2}\rangle=\frac{W_{0}^{2}}{\Omega^{2}}\bigg{(}1+\frac{\pi}{8}\ L\ n_{0}\ W_{0}^{2}\frac{h}{L}+2.6\ \sigma_{R}^{2}\ \Omega^{\frac{5}{6}}\frac{h}{L}\bigg{)}$$

$$\langle\Delta W_{i}^{2}\Delta W_{j}^{2}\rangle=(2\delta_{ij}-0.8)\ \frac{W_{0}^{4}}{\Omega^{\frac{19}{6}}}\bigg{(}1+\frac{\pi}{8}\ L\ n_{0}\ W_{0}^{2}\frac{h}{L}\bigg{)}\sigma_{R}^{2}\frac{h}{L}\ .$$

$$⟨x_0⟩=⟨y_0⟩=0⟨x_0^2⟩=⟨y_0^2⟩=\alpha L$$

$$\langle W_{i}^{2}\rangle=\frac{W_{0}^{2}}{\Omega^{2}}\bigg{(}1+\frac{\pi}{24}\ L\ n_{0}\ W_{0}^{2}\bigg{(}\frac{h}{L}\bigg{)}^{3}+1.6\ \sigma_{R}^{2}\ \Omega^{\frac{5}{6}}\bigg{(}\frac{h}{L}\bigg{)}^{\frac{8}{3}}\bigg{)}$$

$$\langle\Delta W_{i}^{2}\Delta W_{j}^{2}\rangle=(2\delta_{ij}-0.8)\frac{3}{8}\frac{W_{0}^{4}}{\Omega^{\frac{19}{6}}}\bigg{(}1+\frac{\pi}{24}\ L\ n_{0}\ W_{0}^{2}\bigg{(}\frac{h}{L}\bigg{)}^{3}\bigg{)}\sigma_{R}^{2}\bigg{(}\frac{h}{L}\bigg{)}^{\frac{8}{3}}\ ,$$

$$\bar{R}={\int }_0^{1}R(\eta )\mathcal{P}(\eta )d\eta =\sum _{i=1}^{N_{\mathrm{bins}}}R({\eta }_i)\mathcal{P}({\eta }_i).$$

$$u(\mathbf{\rho },L)$$

$$u_{0}(\boldsymbol{\rho})=\sqrt{\frac{2}{\pi W_{0}^{2}}}\ {\rm exp}\bigg{[}-\frac{1}{W_{0}^{2}}|\boldsymbol{\rho}|^{2}-\frac{ik}{2F}|\boldsymbol{\rho}|^{2}\bigg{]}\ ,$$

$$\displaystyle G_{0}(\boldsymbol{\rho},\boldsymbol{\rho}^{\prime}:z,z^{\prime})=\frac{k}{2\pi i(z-z^{\prime})}\ {\rm exp}\bigg{[}\frac{ik|\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime}|^{2}}{2(z-z^{\prime})}\bigg{]}$$

$$S(\boldsymbol{\rho},\boldsymbol{\rho}^{\prime};z,z^{\prime})=\frac{k}{2}\int_{z^{\prime}}^{z}d\zeta\ \delta\varepsilon\bigg{(}\boldsymbol{\rho}\frac{\zeta-z^{\prime}}{z-z^{\prime}}+\boldsymbol{\rho}^{\prime}\frac{z-\zeta}{z-z^{\prime}},\zeta\bigg{)}\ .$$

$$\delta \epsilon =\delta {\epsilon }_{\mathrm{turb}}+\delta {\epsilon }_{\mathrm{scat}}.$$

$${\Phi }_{\epsilon }(\mathbf{K})={\Phi }_{\epsilon }^{\mathrm{turb}}(\mathbf{K})+{\Phi }_{\epsilon }^{\mathrm{scat}}(\mathbf{K}),$$

$$⟨\delta \epsilon ({\mathbf{r}}_1)\delta \epsilon ({\mathbf{r}}_2)⟩=\int d^3\mathbf{K}{\Phi }_{\epsilon }(\mathbf{K})\exp [i\mathbf{K}\cdot ({\mathbf{r}}_1-{\mathbf{r}}_2)].$$

$$⟨\delta \epsilon ({\mathbf{r}}_1)\delta \epsilon ({\mathbf{r}}_2)⟩=2\pi \delta (z_1-z_2)\int d^2\mathbf{k}{\Phi }_{\epsilon }(\mathbf{k})\exp [i\mathbf{k}\cdot ({\mathbf{\rho }}_1-{\mathbf{\rho }}_2)]$$

$${\Phi }_{\epsilon }^{\mathrm{turb}}(\mathbf{k})=0.132C_n^2|\mathbf{k}{|}^{-\frac{11}{3}}.$$

$$\Phi_{\varepsilon}^{\rm scat}({\bf k})=\frac{n_{0}\zeta_{0}^{4}}{8\pi k^{2}}\ {\rm exp}\big{[}-\zeta_{0}^{2}|{\bf k}|^{2}\big{]}\ ,$$

$$\mathbf{v}=(x_0,y_0,W_1,W_2,{\phi }_0).$$

$$\Theta_{i}=\ln\bigg{(}\frac{W_{i}^{2}}{W_{0}^{2}}\bigg{)}\quad i=1,2\ ,$$

$$⟨x_0^2⟩=⟨y_0^2⟩={\int }_{{\mathbb{R}}^4}{d}^2{\mathbf{\rho }}_1{d}^2{\mathbf{\rho }}_2{x}_1{x}_2{\Gamma }_4({\mathbf{\rho }}_1,{\mathbf{\rho }}_2;L).$$

$$\displaystyle\langle W_{1/2}^{2}\rangle=4\bigg{[}\int_{\mathbb{R}^{2}}d^{2}\boldsymbol{\rho}\ x^{2}\ \Gamma_{2}(\boldsymbol{\rho};L)-\langle x_{0}^{2}\rangle\bigg{]}\ ,$$

$$\displaystyle\langle\Delta W_{i}^{2}\Delta W_{j}^{2}\rangle=-8\biggl{\{}2\bigg{(}\int_{\mathbb{R}^{2}}d^{2}\boldsymbol{\rho}\ x^{2}\ \Gamma_{2}(\boldsymbol{\rho};L)\bigg{)}^{2}$$

$$-{\int }_{{\mathbb{R}}^4}{d}^2{\mathbf{\rho }}_1{d}^2{\mathbf{\rho }}_2[x_1^2{x}_2^2(4{\delta }_{ij}-1)-x_1^2{y}_2^2(4{\delta }_{ij}-3)]$$

$$\displaystyle\times\Gamma_{4}(\boldsymbol{\rho}_{1},\boldsymbol{\rho}_{2};L)\bigg{\}}-16[4\delta_{ij}-1]\langle x_{0}^{2}\rangle^{2}\ ,$$

$${\Gamma }_2(\mathbf{\rho }$$

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

$${\mathcal{D}}_S={\mathcal{D}}_S^{\mathrm{turb}}+{\mathcal{D}}_S^{\mathrm{scat}}.$$

$${\mathcal{D}}_S^{\mathrm{turb}}(\mathbf{\rho },{\mathbf{\rho }}^{\prime })=2.95k^{2}L{\int }_0^{1}d\xi C_n^2(\xi )|\mathbf{\rho }\xi +{\mathbf{\rho }}^{\prime }(1-\xi ){|}^{\frac{5}{3}}$$

Full text

Keywords: Satellite links, Quantum Key Distribution, Atmospheric effects

1 Introduction

Quantum Key Distribution (QKD) and Quantum Communications 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, 4, 5, 6, 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 [9, 10, 8, 11] and several experiments have definitely proved that the technology involved is ready for deployment [12, 13, 14, 15, 16].

Optical satellite-based links have the important drawback of being strongly dependent on the weather conditions [17, 18, 19, 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 performances 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 performances 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 Sec. 2 we introduce the problem of free-space optical links and an analytical method to study them. The discussion continues in A. In Sec. 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 Sec. 4 we use the analysis of the transmittance of the channel conducted in the previous section to study the performances of different QKD implementations, in some interesting real-life scenarios. The analysis concerning the use of smaller and more affordable satellites is performed in Sec. 5. Finally, the results are summarized and discussed in Sec. 6. The appendix starts with a recap of the results of [21, 22] (A) and their application to the problem at hand (B). Then two models for the estimation of the stray light satellite links [11, 23] are presented in C. D is devoted to the definition of the QKD protocols we use in Sec. 4 and the expression of the correspondent key rates. In E we report the parameters chosen for the simulations and we discuss their pertinence.

2 Free-space optical links and the Elliptic Beam Approximation

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)

[TABLE]

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 $\boldsymbol{\rho}$ 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 Fig 1)

[TABLE]

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.

The fluctuations of the relative permittivity of the atmospheric air can be statistically modeled [24, 25, 26, 27, 28, 29, 30, 31, 32]. The probability distribution of the parameters in Eq. (2) can then be analytically estimated, as shown in [21, 22]. A brief recap of the derivation and the main results is presented in 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 Sec. 4).

3 Satellite-based links: model and results

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 Fig.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 to 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^{\circ}]$. The effective thickness of the atmosphere $\bar{h}$ is fixed here to $20$ km, for the considerations above.

As introduced in Sec.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 A and B. First of all we need to evaluate Eqs. (23), (24) and (25) in order to compute the moments of the distributions of the elliptic beam parameters (Eq. (2)). To do so, an integration along the propagation path must be performed (Eqs. (30) and (30)). 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

[TABLE]

where $C_{n}^{2}$ is the value of the refractive index structure constant and $n_{0}$ is the density of scattering particles. $\Theta(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 Fig. 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}^{2}$, are available in the literature [33, 34, 35, 36, 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 E). Additional extinction losses due to back-scattering and absorption in the atmosphere are modeled by a parameter $\chi_{\rm ext}$, as described in 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}^{2}$ and $n_{0}$ completely describe the atmospheric conditions together with the thickness $\bar{h}$ and the extinction factor $\chi_{ext}$.

Following the analysis of A (in particular equations Eq. (23), (24), (25)), we compute the first and second moments of the beam parameters in Eq. (2) for the link described in Eq. (3). The distribution of the angle of orientation of the elliptical profile $\varphi_{0}$ is assumed uniform in $[0,\pi/2]$ as in [21, 22]. The mean value and variance of the beam centroid position are the same for $x$ and $y$ directions and equal to (Eq. (23))

[TABLE]

where the quantity $\sigma_{R}^{2}=1.23\ C_{n}^{2}\ k^{\frac{7}{6}}L^{\frac{11}{6}}$ is the so-called Rytov parameter and $\Omega=\frac{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^{2}_{i}$ with $i=1,2$, are instead estimated from Eq. (24) and (25)

[TABLE]

[TABLE]

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

[TABLE]

and for the semi-axes of the elliptical profile

[TABLE]

[TABLE]

where $\alpha\sim 2\ \mu$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 $\frac{h}{L}$. As $\frac{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_{0}$. 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 Eq. (A) and random sampling. Two examples are shown in Fig. 3 and Fig. 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 Fig. 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.

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, just like in Fig. 3 and Fig. 4, from 1000 samples of the parameters in Eq. (2) and using Eq. (A). 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.

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_{\rm sat}=30$ cm for the orbiting one, $D_{\rm grnd}=1$ m for the ground station telescope and $\lambda=800$nm. These are demanding values, consistent with the Chinese mission Micius (see [12, 13, 14, 15] for details). Further analysis are reported in Sec.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_{0}$, 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.

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 [39] 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 [40, 41, 42]. 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 [43, 44], Gamma-Gamma [45] and Double Weibull [46] 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 [47, 48, 49]: 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. 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.

4 Performances of QKD implementations

The transmittance shown in Fig. 5 and 6 can now be used to compute the expected secret key rates of a QKD protocol. In the following we analyze the performances of the BB-84 protocol [50] 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 [51] and WCPs with decoy states [52, 53, 54, 55], taking into account finite-key effects. The key rates are averaged over the PDT computed for different link lengths and configurations

[TABLE]

Here $\bar{R}$ is the averaged key rate, $R(\eta)$ 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_{\rm bins}$ bins, centered in $\eta_{i}$ for $i=1,\ N_{\rm bins}$, and taking the weighted sum of the rates. $\mathcal{P}(\eta_{i})$ is estimated through random sampling, as pointed out in Sec. 3. The expressions for the key rates $R(\eta)$ for the different implementations are given in D, see Eq. (D) for SPs and Eq. (D) 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 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 Fig. 7 and Fig. 8, for both night-time and day-time operation, under good weather conditions, corresponding to situation 1 in Fig. 5. In the following we set the block sizes at $10^{6}$ for SPs and at $10^{8}$ 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^{5}$ for SPs and at $10^{7}$ for WCPs.

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^{\circ}$ 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 (refer to Tab. 3 in E).

Up-links have poorer performances 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 (Tab. 3 in 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 quantum bit error rate too high and the key rate vanishes, therefore we omit the corresponding graph. We stress that here (we refer to 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 Fig. 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.

5 Cube-sat performance analysis

The simulations reported in Sec. 3 and Sec. 4 assume a quite demanding value of the optical aperture of the orbiting telescope. It is compatible with the Micius satellite [12, 13, 14, 15], operated by the Chinese Academy of Science, as part of the Quantum Experiments at Space Scale (QUESS) 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., Cube-sats [56, 57, 58]) for QKD implementation [9, 59, 60, 61, 62]. The possibility to deploy many of such satellites in a single mission, or to share the vector with other payloads, lowers considerably the launch cost of these devices. They must be loaded with much smaller optics, of diameter $\leq 10$ cm, even if bigger aperture can be achieved with the use of deployable optics. When used as transmitter, in the down-link configuration, the smaller aperture creates beams with much higher intrinsic divergence than the case studied in Sec. 3. In the up-link configuration, instead, smaller transmittance is due to the smaller collecting area. We show in Fig. 9 the results of the link simulation for down-links and up-links, in good weather conditions.

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 Quantum Bit Error Rate for up-links, then, will be almost independent of the diameter of the receiving telescope.

The key rates achievable for nano-satellites in the down-link and up-link configurations are reported in Fig. 10 and Fig. 11. As expected, the range of angles over which a non-zero key rate can be exchanged shrinks with respect to the case of Sec. 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 can’t be exchanged in a single transit of the Cube-sat.

6 Conclusion

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 performances can be more precisely mapped, in order to optimize the structure of future global quantum networks [63, 64, 65, 66, 67]. Through the use of this model, the data from meteorological predictions can directly be linked to the performances 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 performances 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 communications 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 its performances 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.

Appendix A Free-space links with turbulence and scatterers

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 Sec. 3.

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

[TABLE]

Since the losses due to back-scattering and absorption can’t 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_{\rm ext}}$. The extinction factor $\chi_{\rm ext}\in[0,1]$ accounts for absorption and back-scattering losses and can be considered as a non-fluctuating quantity (see [22]). In Eq. (11) $u_{0}(\boldsymbol{\rho}^{\prime})$ is the Gaussian envelope at the transmitter plane ($z=0$, $z$ is the longitudinal coordinate)

[TABLE]

with $W_{0}$ the beam spot radius at the transmitter, $k$ the optical wavenumber and $F$ the focal length of the beam. $G_{0}$ is a Gaussian integral kernel

[TABLE]

while $S$ contains all the atmospheric effects

[TABLE]

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 $\delta\varepsilon$ can be separated in two contributions, related to turbulence and scattering

[TABLE]

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

[TABLE]

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

[TABLE]

In the previous equations ${\bf 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]

[TABLE]

where ${\bf k}$ represents the momentum in the plane transverse to the propagation direction. $\boldsymbol{\rho}_{1}$ and $\boldsymbol{\rho}_{2}$ are the components of the vectors ${\bf r}_{1}$ and ${\bf r}_{2}$ in the transversal plane, while $\delta(z)$ is the Dirac-delta. The Kolmogorov model allows us to write the turbulence-related part of the relative permittivity fluctuation spectrum as [24, 25, 26, 27]

[TABLE]

The refractive index structure constant $C_{n}^{2}$ characterizes the strength of turbulence in the optical domain and is an important parameter of the model. The scattering term in Eq. (16) can be approximated as a Gaussian function [28, 29, 30, 31, 32]

[TABLE]

with $\zeta_{0}$ correlation length of the fluctuations due to scattering particles. Here $n_{0}$ 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 Sec. 1

[TABLE]

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

[TABLE]

where $W_{0}$ 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},\Theta_{1},\Theta_{2}$ is Gaussian, while the angle of orientation $\varphi_{0}$ is uniformly distributed in $[0,\pi/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_{0},y_{0})$, we can choose the reference frame such that $\langle x_{0}\rangle=\langle y_{0}\rangle=0$ and [25, 69]

[TABLE]

Here $\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}^{2}$ have the following form (see Appendix of [21] for details)

[TABLE]

where the second-order field correlation function $\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 Eq. (11). We report the calculations only for $\Gamma_{2}(\boldsymbol{\rho};L)$, the equivalent but more cumbersome expressions for $\Gamma_{4}(\boldsymbol{\rho}_{1},\boldsymbol{\rho}_{2};L)$ can be found in [22], Appendix B. Substituting Eq. (11) in the definition of $\Gamma_{2}(\boldsymbol{\rho};z)$ yields

[TABLE]

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 Eq. (14))

[TABLE]

Substituting Eq. (14) and exploiting again the Markov approximation, the factorization in Eq. (15) and (16) can be carried over

[TABLE]

We can now introduce the models for the permittivity fluctuations spectrum related to turbulence (Eq. (19)) and scatterers (Eq. (20)), obtaining

[TABLE]

where we introduced the rescaled longitudinal coordinate $\xi\in[0,1]$, where $\xi=1$ corresponds to $z=L$. We allowed for a dependence on longitudinal coordinate in $C_{n}^{2}(\xi)$ and $n_{0}(\xi)$ for later use. We recall the definition of the so-called Rytov parameter $\sigma_{R}^{2}=1.23\ C_{n}^{2}\ k^{\frac{7}{6}}L^{\frac{11}{6}}$. Substituting in Eq. (26) the definition of the Gaussian envelope $u_{0}(\boldsymbol{\rho})$ (Eq. (12)) and the integral kernel $G_{0}(\boldsymbol{\rho},\boldsymbol{\rho}^{\prime}:L,0)$ (Eq. (13)), the second-order field correlation function reads

[TABLE]

Here \alpha=1+\Omega^{2}\big{(}1-\frac{L}{F}\big{)}^{2} with the Fresnel number defined as $\Omega=\frac{kW_{0}^{2}}{2L}$.

Integrating Eq. (A) (and the equivalent one for $\Gamma_{4}$) and then Eqs. (24), we obtain the first and second moments of the probability distribution of $W_{i}^{2}$. Then, the moments for the variables $\Theta_{i}$ are easily obtained from Eq. (22).

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

[TABLE]

with

[TABLE]

In the previous equations $(\rho,\theta)$ 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 Probability Distribution of the Transmittance (PDT) is then easily reconstructed. Extract at random $M$ 5-tuples of values for $(x_{0},y_{0},\Theta_{1},\Theta_{2},\varphi_{0})$, according to the correct probability distribution. Compute first the values of the ellipse semi-axes $W_{i}$ from $\Theta_{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 Sec. 3 of the main text (Fig. 3 and Fig. 4).

Appendix B Application of the model to a satellite-based link

In this section we are going to apply the model described in the previous section to a satellite-based link, as described in Sec. 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 A. In particular, we will focus on the first term of the quantity $\langle W_{1/2}^{2}\rangle$ defined in equation Eq. (24), which only contains the second order correlation function $\Gamma_{2}(\boldsymbol{\rho};L)$. The computations involving the integration of the fourth order correlation function $\Gamma_{4}(\boldsymbol{\rho}_{1},\boldsymbol{\rho}_{2};L)$ are much more cumbersome and will not be reported here.

Inserting Eq. (A) in the first term of Eq. (24) we obtain the following integration, where all the quantities are defined in the previous section A

[TABLE]

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

[TABLE]

In the rescaled longitudinal coordinate $\xi$, the conditions in Eq. (3) in the main text become

[TABLE]

Inserting Eq. (B) (we consider down-link in this example) into Eq. (35) and (36) we can solve the integration and obtain

[TABLE]

where $\sigma^{2}_{R}$ has been defined in Sec. 3. From this passage we clearly see where the dependency on $\frac{h}{L}$ in Eq. (4) to Eq. (9) originates from. When we introduce Eq. (38) and Eq. (39) in Eq. (34), we recognize that it only contains Gaussian integrals of the form

[TABLE]

with $c=\{0,2\}$ and that can be readily solved. The only exception is the turbulence term, which contains $|\boldsymbol{\rho}^{\prime}|^{\frac{5}{3}}$. We can simplify the computation introducing the approximation [22, 69] $|\boldsymbol{\rho}^{\prime}/W_{0}|^{\frac{5}{3}}\simeq|\boldsymbol{\rho}^{\prime}/W_{0}|^{2}$. Then, one just has to solve the multiple Gaussian integrals and insert it in Eq. (24) to obtain the value of $\langle W_{1/2}^{2}\rangle$ for down-links, as in Eq. (8). Similar techniques can be used to compute all the other moments of the beam variables.

Eq. (4), (5), (7) and (8) have been computed specifically for the problem at hand, the non-uniform link described at the beginning of Sec. 3. Eq. (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 Eq. (4), (5), (7) and (8) the corrections due to the non-uniformity of the link (of the form $A(h/L)^{\beta}$, where $A$ is a constant and $\beta=\{1,8/3,3\}$) act like multiplicative factors on the parameters $\sigma^{2}_{R}$ and $n_{0}$. So, we started from the calculation of the quantity $\langle\Delta W_{i}^{2}\Delta W_{j}^{2}\rangle$ in [22] and attached the multiplicative corrections found above, in order to obtain Eq. (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 Eq. (6) and (9), as the interplay between beam wandering (Eq. (4) and (7)) and beam spreading (Eq. (5) and (8)) is much more significant in this context. Finally, we point out that the computation of the quantities in Eq. (4), (5), (7) and (8) have been carried out without the introduction of the weak turbulence approximation used in [21, 22]. For Eq. (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.

Appendix C Error model and environmental photons

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 expected 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]

[TABLE]

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_{\rm sun}$ is the solar spectral irradiance in ${\rm photons}\ {\rm s}^{-1}{\rm nm}^{-1}{\rm m}^{-2}$ at the wavelength of interest. $\Omega_{\rm fov}$ and $a$ are angular field of view and radius of the receiving telescope. $B_{f}$ is the width of the spectral filtering and $\Delta 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]

[TABLE]

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 Eq. (42) we derive the number of photons per time window

[TABLE]

where $h$ is the Planck constant and $\nu$ 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}\ {\rm sr}\ \mu{\rm m}$ during a full-Moon night and $H_{b}=1\ W\ {\rm m}^{-2}\ {\rm 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 Quantum bit error rate (QBER) is computed assuming the noisy photons to be completely unpolarized

[TABLE]

Here $Q_{0}$ 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_{0}=2\%$. $N_{\rm noise}$ and $N_{\rm sig}$ are, respectively, the expected 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 $\Omega_{\rm fov}$, $B_{f}$ and $\Delta 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.

Appendix D Rates for BB-84 with single photons and Weak Coherent Pulses

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=\{\ket{0},\ket{1}\}$ or $Z=\{\ket{+},\ket{-}\}$ at random, with $\ket{\pm}=(\ket{0}\pm\ket{1})/\sqrt{2}$. B measures the received qubits in the bases $X$ or $Z$, at random. The results of [51] state that a secret key of length $l$ can be shared, if

[TABLE]

out of $n$ successfully exchanged single photon signals, where the function $h_{2}$ 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_{\rm 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 ${\rm 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 $\varepsilon_{\rm sec}$ and $\varepsilon_{\rm cor}$.

An alternative protocol based on decoy states [52, 54] is used when the source emits Weak Coherent Pulses instead of real single photons. We follow the analysis of [53, 55], 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

[TABLE]

$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. $\phi_{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 [55]). 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$.

Appendix E Choice of parameters for the satellite-based link

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

The parameters $C_{n}^{2}$, $n_{0}$ 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}^{2}$ are available in the literature [34, 35, 36, 37]. The widely used parametric fit due to Hufnagel and Valley [34, 35] reliably replicates the behaviour of $C_{n}^{2}$ in mid-latitude climate

[TABLE]

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\ 10^{-14}\ {\rm m}^{-2/3}$ and $v=21\ {\rm m}/{\rm s}$, although $v=57{\rm m}/{\rm s}$ is sometimes used for stronger wind conditions. The value of $C_{n}^{2}$ inside the atmosphere in our model is estimated by the integral average of this function in $[0,\infty]$, rescaled by the fixed thickness $h$

[TABLE]

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

The scattering particles described by the density $n_{0}$ mainly consist of water droplets, so, in order to estimate the value of $n_{0}$, we start from the profile of the water vapour content in the atmosphere. The absolute humidity vertical profile $\tau(z)$ in the range $[0,10\ {\rm km}]$ can be written as a double exponential [70, 71]

[TABLE]

with the two scale heights $H_{1}$ and $H_{2}$. The contribution of the region with $z>10\ {\rm km}$ is rather low and we neglect it here. The parameters $H_{1}$ and $H_{1}$ 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 U. S. Standard Atmosphere (1962) [72], $H_{1}=2.243$ and $H_{2}=1.414$. We obtain a rescaling factor $\omega$ in the same way as we did in the previous case

[TABLE]

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

The extinction factor $\chi_{\rm ext}(\theta)$ varies as a function of the elevation angle in the following way

[TABLE]

The value of the parameter $\beta$ reported in Tab. 1 has been chosen to match the amount of extinction used in [10], based on the MODTRAN5 software [38].