Optimal design and performance evaluation of free-space quantum key distribution systems

The first term entering the channel efficiency in equation (1) takes into account the attenuation due to atmospheric absorption. The channel absorption efficiency η_A for a link distance z depends on the absorption coefficient A(λ) for a specific wavelength λ and can be modeled as

$$\begin{equation}{\eta }_{\text{A}}=1{0}^{-A(\lambda )\cdot z},\end{equation} \tag{ 4 }$$

where the absorption coefficient is assumed constant for a horizontal link. An established tool for calculating the spectral properties of the atmosphere is the LOWTRAN software package [33], which can be used to predict atmospheric absorption and scattering over a wide wavelength range, on horizontal or slanted paths, taking into account both geographical and seasonal atmospheric variations. In figure 2, we show the atmospheric absorption coefficient computed by LOWTRAN as a function of wavelength for a horizontal link, considering a sub-arctic winter atmosphere. Clearly, signal wavelengths for which the atmospheric absorption coefficient is high should be avoided in free-space links.

Zoom In Zoom Out Reset image size

The second term entering the overall channel efficiency in equation (1) is the collection efficiency ${\eta }_{{D}_{\text{Rx}}}$ and takes into account the finite size of the receiver aperture D_Rx that is used to collect the incoming beam. We now detail the impact of turbulence-induced beam broadening, and beam wander and scintillation on the collection efficiency ${\eta }_{{D}_{\text{Rx}}}$.

2.2.1. Turbulence-induced beam broadening

In vacuum, the beam size W(z) of a collimated Gaussian beam of waist W₀ and wavelength λ propagating for a distance z is given by the formula for diffraction-limited propagation:

$$\begin{equation}W(z)={W}_{0}\sqrt{1+{\left(\frac{\lambda z}{\pi {W}_{0}^{2}}\right)}^{2}}.\end{equation} \tag{ 5 }$$

According to the Kolmogorov's theory of turbulence [34], when the propagation happens through the turbulent air, the index of refraction $\text{n}$ is treated as a fluctuating random field around a mean value which induces a perturbation of the wavefront, resulting in an overall loss of coherence of the optical wave. The atmospheric perturbation is captured by the so-called power spectral density Φ_n(κ), which is the Fourier transform of the refractive-index covariance function in terms of the spatial frequency κ. In our model, we use the power spectral density Φ_n(κ) for refractive-index fluctuations given by the well-known Kolmogorov spectrum of atmospheric turbulence

$$\begin{equation}{{\Phi}}_{\text{n}}(\kappa )=0.033{C}_{\text{n}}^{2}{\kappa }^{-11/3},\end{equation} \tag{ 6 }$$

which is widely used in theoretical calculations. The strength of the turbulence is parametrized by the refractive-index structure constant ${C}_{\text{n}}^{2}$, which may be considered constant along a horizontal link.

The effect of turbulence-induced coherence loss on Gaussian beam propagation has been studied by [35, 36], who found that the following formula holds:

$$\begin{equation}W(z)={W}_{0}\sqrt{1+\left(1+\frac{2{W}_{0}^{2}}{{\rho }_{0}^{2}(z)}\right){\left(\frac{\lambda z}{\pi {W}_{0}^{2}}\right)}^{2}},\end{equation} \tag{ 7 }$$

where

$$\begin{equation}{\rho }_{0}(z)={(0.55{C}_{\text{n}}^{2}{k}^{2}z)}^{-3/5}\end{equation} \tag{ 8 }$$

is the spherical-wave atmospheric spatial coherence radius (with k = 2π/λ the wave-number). A parameter typically used to quantify the strength of the turbulence on a channel is the atmospheric coherence width r₀, also known as the Fried parameter [34]. It describes the maximum size of a telescope aperture whose resolution is comparable to the diffraction limited one under that specific turbulence, and is related to ρ₀ by:

$$\begin{equation}{r}_{0}=2.1{\rho }_{0}.\end{equation} \tag{ 9 }$$

While for the diffraction-limited case a larger beam waist W₀ implies a smaller intrinsic divergence θ₀ = λ/πW₀, in the turbulence-affected case the larger the ratio W₀/ρ₀, the stronger is the loss of coherence. Indeed, for z ≫ 1 we have that the two competing effects, diffraction-broadening and turbulence broadening, compensate each other and the beam size tends to a value that is independent of the initial waist W₀:

$$\begin{equation}W(z)\stackrel{z\gg 1}{\sim }\frac{\lambda \sqrt{2}}{\pi }{\left(0.55{C}_{\text{n}}^{2}{k}^{2}\right)}^{3/5}{z}^{8/5},\end{equation} \tag{ 10 }$$

as shown in figure 3 for different values of W₀ and ${C}_{\text{n}}^{2}$.

Zoom In Zoom Out Reset image size

Assuming for simplicity that the pointing error is negligible, the average contribution to the collection efficiency ${\eta }_{{D}_{\text{Rx}}}$ caused by diffraction and beam broadening for a receiver of finite aperture diameter D_Rx is given by

$$\begin{equation}\left\langle {\eta }_{{D}_{\text{Rx}}}\right\rangle =1-\mathrm{exp}\left[-\frac{{D}_{\text{Rx}}^{2}}{2W{(z)}^{2}}\right],\end{equation} \tag{ 11 }$$

which is the integral of a Gaussian distribution of standard deviation W(z)/2 over a concentric circular area of diameter D_Rx. Although the easiest way to optimize this term is to increase D_Rx, this strategy decreases the fiber coupling efficiency, as we will study in section 2.3.4, and it might not be a viable solution since it requires larger and more expensive optics.

2.2.2. Beam wander and scintillation

A consequence of the asymptotic trend of W(z) shown in figure 3 is that, for relatively long links, the average collection efficiency becomes independent of the choice of the waist at the transmitter. However, this does not necessarily imply that a smaller transmitter waist should be preferred when, faced with equivalent performance at the receiver, a smaller telescope is less complex and less costly.

This is because the beam-size W(z) appearing in equation (11) is the so-called long-term spot size, which represents the size of the beam averaged over a timescale much longer than the turbulence dynamics. The instantaneous short-term (ST) beam size W_ST(z) at a distance z is given by:

$$\begin{equation}{W}_{\text{ST}}(z)=\sqrt{{W}^{2}(z)-\left\langle {r}_{\text{c}}^{2}\right\rangle },\end{equation} \tag{ 12 }$$

where $\left\langle {r}_{\text{c}}^{2}\right\rangle$ is the variance of beam-wander fluctuations at the receiver aperture plane. Beam wandering is caused by the larger-sized turbulence eddies, and results in a shift of the ST spot on the receiver aperture. For a collimated beam transmitted along a horizontal channel one has [34]:

$$\begin{equation}\left\langle {r}_{\text{c}}^{2}\right\rangle =2.42{C}_{\text{n}}^{2}{z}^{3}{W}_{0}^{-1/3}.\end{equation} \tag{ 13 }$$

We can then see that, since the beam wandering depends on ${W}_{0}^{-1/3}$, a smaller transmitter waist always leads to a higher beam wandering variance at the receiver. The effect of beam wandering may be removed by automatically tracking and correcting the tip/tilt fluctuation at both the transmitter and receiver using of two counter-propagating beacons [37, 38]. This link configuration assumes a symmetry between the front- and back-propagating channels, and its treatment for QKD links is covered in [29]. In this work, we will restrict ourselves to the case of a single pointing feedback at the receiver, which does not enable the removal of beam wandering.

In addition to beam wander, one must also consider the effect of atmospheric scintillation, which introduces random fluctuations in the beam irradiance profile. The pdf ${p}_{{D}_{\text{Rx}}}$ of ${\eta }_{{D}_{\text{Rx}}}$ was derived analytically in [28, 29] exploiting the law of total probability and separating the contributions from turbulence-induced beam wandering and atmospheric scintillation. The distribution ${p}_{{D}_{\text{Rx}}}$ varies depending on the strength of turbulence, which is parametrized by the Rytov variance ${\sigma }_{\text{R}}^{2}$, which is defined as:

$$\begin{equation}{\sigma }_{\text{R}}^{2}=1.23{C}_{\text{n}}^{2}{k}^{7/6}{z}^{11/6}.\end{equation} \tag{ 14 }$$

For weak turbulence $({\sigma }_{\text{R}}^{2}< 1)\ {p}_{{D}_{\text{Rx}}}$ resembles a log-negative Weibull distribution, while for stronger turbulence $({\sigma }_{\text{R}}^{2} > 1)$ one finds a truncated log-normal distribution. The exact form and derivation of ${p}_{{D}_{\text{Rx}}}$ can be found in [29] and requires the knowledge of the following quantities: the average collection efficiency $\left\langle {\eta }_{{D}_{\text{Rx}}}\right\rangle$ of (11), the collection efficiency calculated using the ST waist W_ST of (12), the beam-wander variance $\left\langle {r}_{\text{c}}^{2}\right\rangle$, and the mean-squared efficiency $\left\langle {\eta }_{{D}_{\text{Rx}}}^{2}\right\rangle$, that is:

$$\begin{equation}\left\langle {\eta }_{{D}_{\text{Rx}}}^{2}\right\rangle =\left\langle {\eta }_{{D}_{\text{Rx}}}\right\rangle (1+{\sigma }_{\text{I}}^{2}({D}_{\text{Rx}})),\end{equation} \tag{ 15 }$$

where ${\sigma }_{\text{I}}^{2}({D}_{\text{Rx}})$ is the aperture-averaged scintillation index (flux variance) [34]:

$$\begin{equation}{\sigma }_{\text{I}}^{2}({D}_{\text{Rx}})=\mathrm{exp}\left[\frac{0.49{\tilde{{\beta }_{0}}}^{2}}{{\left(1+0.18{d}^{2}+0.56{\tilde{{\beta }_{0}}}^{12/5}\right)}^{7/6}}+\frac{0.51{\tilde{{\beta }_{0}}}^{2}{\left(1+0.69{\tilde{{\beta }_{0}}}^{12/5}\right)}^{-5/6}}{1+0.90{d}^{2}+0.62{d}^{2}{\tilde{{\beta }_{0}}}^{12/5}}\right]-1,\end{equation} \tag{ 16 }$$

with $d=\sqrt{\frac{k{D}_{\text{Rx}}^{2}}{4z}}$ and ${\tilde{{\beta }_{0}}}^{2}=0.4065{\sigma }_{\text{R}}^{2}$. Some examples of probability distributions ${p}_{{D}_{\text{Rx}}}$ calculated in this way are provided in figure 4.

Zoom In Zoom Out Reset image size

The third term entering the channel efficiency in equation (1) is due to the coupling of the free-space beam into a SMF at the receiver. The fiber coupling efficiency is given by the normalized overlap integral between the fiber mode, and the incident optical field $U(\vec{r},t)$,

$$\begin{equation}U(\vec{r},t)={U}_{0}(\vec{r},t)\mathrm{exp}\left[\chi (\vec{r},t)+i{\Psi}(\vec{r},t)\right],\end{equation} \tag{ 17 }$$

where χ is the log-amplitude perturbation term, and Ψ is the wavefront phase term. The different origins of χ and Ψ perturbations imply the statistical independence of scintillation and phase effects [32, 39], and the coupling efficiency η_SMF can be factorized into three terms:

$$\begin{equation}{\eta }_{\text{SMF}}={\eta }_{0}\ {\eta }_{\text{AO}}\ {\eta }_{\text{S}}\ ,\end{equation} \tag{ 18 }$$

where η₀ is the optical efficiency of the receiver telescope, η_AO is the coupling efficiency due to wavefront perturbations that may be partially corrected by AO, and η_S is the coupling efficiency due to the spatial structure of atmospheric scintillation. We now discuss the three terms separately.

2.3.1. Optical coupling efficiency

The optical coupling efficiency η₀ of an optical system measures the matching between an unperturbed received beam and the MFD of the SMF. η₀ is determined by the design optics of the receiving telescope, and particularly by the ratio α = D_Obs/D_Rx between the diameters of the central obscuration and of the telescope aperture [40]. The ideal coupling efficiency can be parametrized by

$$\begin{equation}{\eta }_{0}(\alpha ,\beta )=2{\left[\frac{\mathrm{exp}(-{\beta }^{2})-\mathrm{exp}(-{\beta }^{2}{\alpha }^{2})}{\beta \sqrt{1-{\alpha }^{2}}}\right]}^{2},\end{equation} \tag{ 19 }$$

where β is given by

$$\begin{equation}\beta =\frac{\pi {D}_{\text{Rx}}}{4\lambda }\frac{\mathrm{M}\mathrm{F}\mathrm{D}}{f},\end{equation} \tag{ 20 }$$

with D_Rx the receiver diameter, f the effective focal length of the optical system, and MFD the MFD of the SMF. Given a particular α, the value of β can be optimized to achieve the optimal coupling efficiency ${\eta }_{0}^{(\mathrm{o}\mathrm{p}\mathrm{t})}={\eta }_{0}({\beta }_{\text{opt}})$. Knowing the value of β_opt allows to choose the optimal design value of f, since typically the working wavelength λ, fiber MFD and receiver diameter are constrained. Figure 5 shows ${\eta }_{0}^{(\mathrm{o}\mathrm{p}\mathrm{t})}$ and β_opt as a function of the obscuration ratio α. For α = 0, we have β_opt = 1.12, that allows to achieve a maximum optical coupling efficiency of 81.5% ≈ −0.89 dB.

Zoom In Zoom Out Reset image size

2.3.2. Effect of wavefront perturbations and AOs correction

As was first derived by [41, 42], the instantaneous wavefront aberration ${\Psi}(\vec{r},t)$ introduced by the turbulent channel at a point $\vec{r}$ on the receiver aperture can be decomposed in a superposition of Zernike polynomials defined over the normalised pupil coordinates (r, φ), with $r=2\vert \vec{r}\vert /{D}_{\text{Rx}}$. In the decomposition of ${\Psi}(\vec{r},t)$, each polynomial term ${Z}_{n}^{m}(r,\varphi )$ of radial degree n and azimuthal degree m is weighted by a time-dependent coefficient ${b}_{n}^{m}(t)$, yielding

$$\begin{equation}{\Psi}(r,\varphi ,t)=\sum\limits _{n,m}{b}_{n}^{m}(t){Z}_{n}^{m}(r,\varphi ).\end{equation} \tag{ 21 }$$

The Zernike coefficient variances $\left\langle {b}_{n}^{m2}\right\rangle$ represent the statistical strength of a particular aberration order, and depend on the ratio of the receiver aperture D_Rx to the Fried parameter r₀ with a modal term scaling with the radial order n [42, 43]:

$$\begin{equation}\left\langle {b}_{n}^{m2}\right\rangle ={\left(\frac{{D}_{\text{Rx}}}{{r}_{0}}\right)}^{\frac{5}{3}}\frac{n+1}{\pi }\frac{{\Gamma}\left(n-\frac{5}{6}\right){\Gamma}\left(\frac{23}{6}\right){\Gamma}\left(\frac{11}{6}\right)\mathrm{sin}\left(\frac{5}{6}\pi \right)}{{\Gamma}\left(n+\frac{23}{6}\right)}.\end{equation} \tag{ 22 }$$

An expression of the instantaneous coupling efficiency in the presence of wavefront perturbations was derived by [32, 44] directly in terms of the Zernike coefficients ${b}_{n}^{m}$:

$$\begin{equation}{\eta }_{\text{AO}}(t)=\mathrm{exp}\left[-\sum\limits _{n,m}{b}_{n}^{m}{(t)}^{2}\right].\end{equation} \tag{ 23 }$$

From (23), since the coefficients are independent, Gaussian-distributed random variables with zero mean and variance given by (22)—i.e. ${b}_{n}^{m}\sim \mathcal{N}(0,\left\langle {b}_{n}^{m2}\right\rangle )$—we derive the average coupling efficiency:

$$\begin{equation}\left\langle {\eta }_{\text{AO}}\right\rangle =\prod\limits _{n,m}\frac{1}{\sqrt{1+2\left\langle {b}_{n}^{m2}\right\rangle }}.\end{equation} \tag{ 24 }$$

Equation (24) makes the calculation of the average SMF coupling efficiency in the presence of a partial AO compensation of turbulence up to an order n_max straightforward. Assuming an ideal AO system with infinite control bandwidth, it is sufficient to completely suppress the coefficients in the productory corresponding to radial orders n ⩽ n_max. Figure 6 shows the average SMF coupling efficiency $\left\langle {\eta }_{\text{AO}}\right\rangle$ as a function of the ratio D_Rx/r₀ for increasing order of AO correction, assuming infinite control bandwidth.

Zoom In Zoom Out Reset image size

From (22), we have that the strength of the aberration orders decreases for higher orders. It is interesting then to view the results of figure 6 in a different form: figure 7 shows the required correction order for different target AO efficiencies, as a function of turbulence strength—parametrized by the ratio D_Rx/r₀. If the receiver diameter is fixed and r₀ is known from a previous characterization of the channel ${C}_{\text{n}}^{2}$ value, then using equation (24) it is possible to derive the number of aberration orders that need to be corrected in order to achieve a given target coupling efficiency.

Zoom In Zoom Out Reset image size

So far, we have assumed that the AO system is able to perfectly correct turbulence in real time. However, this is not the case, as various imperfections affect the efficiency, such as wavefront reconstruction errors and control bandwidth limitations. In this work, we focus on the latter term. As discussed in reference [45], the effect of control bandwidth limitations can be taken into account by introducing a mode attenuation factor ${\gamma }_{n}^{2}$ for the nth order aberration coefficients:

$$\begin{equation}{\gamma }_{n}^{2}=\frac{\int \vert {W}_{n}(\tilde{\nu }){\vert }^{2}\vert \varepsilon (\tilde{\nu }){\vert }^{2}\mathrm{d}\tilde{\nu }}{\int \vert {W}_{n}(\tilde{\nu }){\vert }^{2}\mathrm{d}\tilde{\nu }},\end{equation} \tag{ 25 }$$

where $\vert {W}_{n}(\tilde{\nu }){\vert }^{2}$ represents the power spectral density of the temporal spectrum of the nth order aberrations, $\tilde{\nu }$ is the frequency of the AO loop, and $\varepsilon (\tilde{\nu })$ represents the transfer function between the residual phase and the turbulent wavefront fluctuations, and depends on the AO system's open-loop transfer function $G(\tilde{\nu })$:

$$\begin{equation}\varepsilon (\tilde{\nu })=\frac{1}{1+G(\tilde{\nu })}.\end{equation} \tag{ 26 }$$

In this work, we consider a typical AO system with a pure-integrator control based on wavefront sensing with a Shack–Hartmann wavefront sensor, and correction with a deformable piezo-electric mirror. The open-loop transfer function is then given by

$$\begin{equation}G(\tilde{\nu })={K}_{i}\frac{{\text{e}}^{-\tau \tilde{\nu }}\left(1-{\text{e}}^{-T\tilde{\nu }}\right)}{{(T\tilde{\nu })}^{2}},\end{equation} \tag{ 27 }$$

where K_i is the gain of the integrator, τ is the overall latency of the control-actuator stage, and T is the wavefront sensor integration time.

The temporal power spectra of the different aberration orders have been studied in reference [46], where the authors find that the power spectral density scales polynomially with a cut-off frequency ${\tilde{\nu }}_{\text{c}}^{(n)}$ depending on the radial order n, average wind velocity $\bar{v}$ and receiver diameter D_Rx:

$$\begin{equation}\vert {W}_{n}(\tilde{\nu }){\vert }^{2}\sim \begin{cases}{\tilde{\nu }}^{-2/3}\quad \hfill & \tilde{\nu }\leqslant {\tilde{\nu }}_{\text{c}},n=1\hfill \\ {\tilde{\nu }}^{0}\quad \hfill & \tilde{\nu }\leqslant {\tilde{\nu }}_{\text{c}},n\ne 1\hfill \\ {\tilde{\nu }}^{-17/3}\quad \hfill & \tilde{\nu } > {\tilde{\nu }}_{\text{c}}\hfill \end{cases}\end{equation} \tag{ 28 }$$

with

$$\begin{equation}{\tilde{\nu }}_{\text{c}}^{(n)}=0.3(n+1)\bar{v}/{D}_{\text{Rx}}.\end{equation} \tag{ 29 }$$

The SMF efficiency of equation (24) is thus modified in the case of finite control bandwidth as:

$$\begin{equation}\left\langle {\eta }_{\text{AO}}\right\rangle =\prod\limits _{\begin{subarray}{c}n,m\\ n\leqslant {n}_{\mathrm{max}}\end{subarray}}\frac{1}{\sqrt{1+2{\gamma }_{n}^{2}\left\langle {b}_{n}^{m2}\right\rangle }}+\prod\limits _{\begin{subarray}{c}n,m\\ n > {n}_{\mathrm{max}}\end{subarray}}\frac{1}{\sqrt{1+2\left\langle {b}_{n}^{m2}\right\rangle }},\end{equation} \tag{ 30 }$$

with n_max the maximum aberration order corrected.

To conclude, in figure 8 we show $\left\langle {\eta }_{\text{AO}}\right\rangle$ as a function of the corrected aberration order and wavefront sensor integration time, in a scenario with D_Rx/r₀ = 8.5 and an average wind velocity corresponding to a light breeze, assuming in this instance a negligible latency of the control loop. We can see that to reach the desired coupling efficiency, represented by dashed blue lines in the figure, both the number of corrected orders and the bandwidth of the sensor must be increased, otherwise, if either one of the two is kept fixed, the SMF coupling has an upper bound that cannot be exceeded.

Zoom In Zoom Out Reset image size

2.3.3. Effect of atmospheric scintillation

In addition to phase perturbations, we also take into account the irradiance fluctuations introduced by atmospheric scintillation, which result in a random apodization of the pupil transmittance function, which in turn affects the maximum SMF coupling efficiency. A rigorous calculation of the scintillation contribution η_S to the SMF efficiency, which considers the optical system modulation transfer function and log-amplitude spatial covariance function C_χ (r) can be found in reference [32] and is based on the result of references [39, 41]. Nonetheless, a good approximation for the average value of η_S is only dependent on the on-axis pupil-plane scintillation index ${\sigma }_{\text{I}}^{2}$, which is given by equation (16) in the limit of an infinitesimal aperture i.e., d → 0:

$$\begin{equation}\left\langle {\eta }_{\text{S}}\right\rangle \approx \mathrm{exp}[-{C}_{\chi }(0)]=\mathrm{exp}[-{\sigma }_{\chi }^{2}]={(1+{\sigma }_{\text{I}}^{2})}^{-\frac{1}{4}},\end{equation} \tag{ 31 }$$

where we use the fact that the log-amplitude variance ${\sigma }_{\chi }^{2}$ is related to the scintillation index through ${\sigma }_{\text{I}}^{2}=\mathrm{exp}(4{\sigma }_{\chi }^{2})-1$.

Figure 9 shows a plot of $\left\langle {\eta }_{\text{S}}\right\rangle$ and ${\sigma }_{\text{I}}^{2}$ as a function of the Rytov variance of equation (14), highlighting the behavior of ${\sigma }_{\text{I}}^{2}$, which increases for the weak fluctuation regime, reaches a maximum value in the focusing regime, and then tends to ${\sigma }_{\text{I}}^{2}\sim 1$ in the strong fluctuation regime. The coupling efficiency reduction due to scintillation is contained to ∼1 dB, consistent with the fact that in the wave structure function the phase contribution dominates over the log-amplitude one [39]. This term should anyway be taken into account, as it cannot be corrected by the AO system.

Zoom In Zoom Out Reset image size

2.3.4. Optimum receiver diameter

We have seen that the average collection efficiency $\left\langle {\eta }_{{D}_{\text{Rx}}}\right\rangle$ of equation (11) increases as the receiver diameter increases. Conversely, for a fixed AO correction order n, the SMF coupling efficiency of $\left\langle {\eta }_{\text{AO}}\right\rangle$ (24) decreases with increasing receiver diameter. This leads to a trade-off between the collection efficiency and the fiber-coupling efficiency (as represented in figure 10). Once the wavelength, link distance, transmitter waist, and turbulence strength are fixed, if the maximum order of the AO correction is also fixed—for example by the available number of deformable mirror actuators, or the size of the wavefront sensor—then it is possible to find the optimum receiver diameter D_opt that maximizes the overall channel efficiency η_CH.

Zoom In Zoom Out Reset image size

Figure 11 shows the optimum receiver diameter as a function of link distance, for different orders of aberration correction and a moderate turbulence strength of ${C}_{\text{n}}^{2}=1{0}^{-14}\enspace {\text{m}}^{-2/3}$. For short link distances (shorter than the transmitter Rayleigh distance ${z}_{0}=\pi {W}_{0}^{2}/\lambda$) the beam size does not diverge much, and the fiber-coupling term dominates, causing the optimal diameter to decrease with increasing distance. For link distances z ∼ z₀, the beam size starts increasing, the collection efficiency term dominates, leading to an increasing optimum beam diameter. For longer propagation distances (z ≫ z₀), the decrease in spatial coherence of the beam is more severe and the turbulence term dominates again, leading to a smaller optimum diameter.

Zoom In Zoom Out Reset image size

2.3.5. Probability distribution of coupling efficiency

As anticipated, when the receiver rate approaches the saturation limit of the single-photon detectors, the full statistics of channel efficiency can no longer be ignored, and the whole probability distribution p_CH has to be considered for the expected SKR to be estimated correctly, as we will present in section 2.5.

The derivation of the probability distribution of the SMF coupling efficiency with partial AOs correction can be found in reference [32], which, however, does not include the effect of finite control bandwidth. Since the irradiance fluctuation statistical contribution is already taken into account for the collection efficiency, we restrict the calculation of the SMF coupling distribution to phase distortions only, and present some examples of probability distributions calculated as a function of the maximum corrected aberration order (figure 12), and varying the bandwidth of the AO control loop (figure 13).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given a set of Zernike coefficients {b_j } (where j indicates the OSA/ANSI index j = [n(n + 1) + m]/2) with variances $\langle {b}_{j}^{2}\rangle$, which may be corrected by a compensation factor ${\gamma }_{j}^{2}$—where ${\gamma }_{j}^{2}=0$ for perfect compensation, and ${\gamma }_{j}^{2}=1$ for uncorrected coefficients [32]—we define the quantity ξ(t) as the instantaneous sum of the squared Zernike coefficients:

$$\begin{equation}\xi (t)=\sum\limits _{j}{\,b}_{j}^{2}(t).\end{equation} \tag{ 32 }$$

The probability distribution of ξ, including the effect of finite AO control bandwidth is then, using the result of reference [47]:

$$\begin{equation}{p}_{\xi }(\xi )=\frac{1}{\pi }{\int }_{0}^{\infty }\frac{\mathrm{cos}\left[{\sum }_{j}\frac{1}{2\,}\mathrm{arctan}\left(2{\gamma }_{j}^{2}\langle {b}_{j}^{2}\rangle u\right)-\xi u\right]}{{\prod }_{j}{\left[1+4{u}^{2}{\left({\gamma }_{j}^{2}\langle {b}_{j}^{2}\rangle \right)}^{2}\right]}^{1/4}}\,\mathrm{d}u.\end{equation} \tag{ 33 }$$

The probability distribution of η_SMF is then

$$\begin{equation}{p}_{\text{SMF}}({\eta }_{\text{SMF}}\vert {\eta }_{\mathrm{max}})=\frac{1}{{\eta }_{\text{SMF}}}{p}_{z}\left[\mathrm{log}\left(\frac{{\eta }_{\mathrm{max}}}{{\eta }_{\text{SMF}}}\right)\right],\end{equation} \tag{ 34 }$$

where η_max = η₀ ⋅ η_S is the maximum normalized coupled flux, given by the product of the optical coupling efficiency of the system η₀, and the spatial scintillation term η_S.

Based on the results of reference [32], which calculated the channel loss term due to SMF coupling, and reference [29], which calculated the channel loss due to the finite receiver aperture in the presence of beam broadening, beam wandering and scintillation, we calculate the overall channel transmittance probability distribution considering both contributions (pupil-plane and focal-plane losses) and exploiting the law of total probability to write p_CH(η_CH) as:

$$\begin{equation}{p}_{\text{CH}}({\eta }_{\text{CH}})=\int {p}_{\text{SMF}}({\eta }_{\text{CH}}\vert {\eta }_{0}{\eta }_{\text{S}}{\eta }_{{D}_{\text{Rx}}}){p}_{{D}_{\text{Rx}}}({\eta }_{{D}_{\text{Rx}}})\mathrm{d}{\eta }_{{D}_{\text{Rx}}},\end{equation} \tag{ 35 }$$

where ${p}_{\text{SMF}}({\eta }_{\text{CH}}\vert {\eta }_{0}{\eta }_{\text{S}}{\eta }_{{D}_{\text{Rx}}})$ is the probability of obtaining a normalized flux η_CH in the SMF fiber, given a maximum input normalized flux ${\eta }_{0}{\eta }_{\text{S}}{\eta }_{{D}_{\text{Rx}}}$, through equation (34).

In figure 14 we show some examples of channel probability distributions for eight case studies, with average link losses between −7 dB and −48 dB and input parameters summarized in table 1. In all case studies, we assume an unobstructed receiver aperture, and maximum optical efficiency of η₀ = 81.5%. Cases 1 and 3 correspond to a scenario with a short urban link with strong turbulence, small aperture Tx/Rx telescopes, and mere tip/tilt correction. Cases 4, and 5 correspond to a scenario with moderate turbulence and longer link distance and highlight the effect of a smaller/larger receiver aperture, which leads to a trade-off between large collection efficiency and minimum aberration order corrected. Cases 2 and 6 show the effect of AO on longer, moderately turbulent links with large aperture receivers. Cases 7 and 8 show examples of highly lossy channels.

Zoom In Zoom Out Reset image size

Table 1. Input parameters for the simulation of the eight case studies.

	Case 1	Case 2	Case 3	Case 4	Case 5	Case 6	Case 7	Case 8
⟨ηCH⟩	−7	−15	−17	−23	−25	−38	−43	−48	(dB)
${C}_{\text{n}}^{2}$	10–13	10–14	10–13	10–14	10–14	10–14	10–14	10–14	(m−2/3)
W0	25	60	25	60	60	60	60	25	(mm)
DRx	50.8	200	50.8	200	50	200	400	200	(mm)
z	1	10	2	10	10	20	20	30	(km)
nmax	1	4	1	1	1	1	2	1

The channel efficiency η_CH that we have studied in the previous sections has a significant impact on the key performance indicators of QKD, such as the SKR and error rate. Clearly, a low efficiency decreases the signal rate and increases the relative weight of errors, which both have a negative effect on the SKR. These parameters are estimated using large samples of data acquired during a long experiment, averaging over the distributions p_CH. This might suggest that the fluctuations of η_CH can be neglected and only its mean value is relevant. Yet, this would be equivalent to calculating the expected value of a function by applying it to the expected value of its parameters. Since the functions that model the generation of a key are not all affine, i.e., they are not all compositions of a translation and a linear map, neglecting the fluctuations and using only the mean values is in general incorrect.

A minor non-affine effect is due to the threshold behavior of single-photon detectors, which produce an electrical signal whenever they observe one or more photons. Consequently, the output rate is:

$$\begin{equation}{R}_{0}=\left[R\cdot (1-{\text{e}}^{-\mu \cdot {\eta }_{\mathrm{det}}\cdot {\eta }_{\text{CH}}})+{R}_{\text{bkg}}\cdot {\eta }_{\mathrm{det}}+{R}_{\text{dark}}\right]\cdot (1+{p}_{\text{ap}}),\end{equation} \tag{ 36 }$$

where R is the repetition rate of the source, μ is the mean number of photons per pulse, η_det represents all the fixed efficiency terms at the receiver not accounted in η_CH, such as the detector quantum efficiency or insertion losses of fiber components, R_bkg is the background photon rate, R_dark is the dark rate, and p_ap is the probability that an afterpulse is triggered by another detection (higher-order afterpulses are neglected). This is not an affine function of η_CH, which enters the formula in the exponential term characteristic of the Poisson distribution of the number of photons in laser pulses. Yet, in typical conditions, in which μ ⋅ η_det ⋅ η_CH ≪ 1, it is legitimate to approximate this function to an affine one because $1-{\text{e}}^{-\mu \cdot {\eta }_{\mathrm{det}}\cdot {\eta }_{\text{CH}}}\approx \mu \cdot {\eta }_{\mathrm{det}}\cdot {\eta }_{\text{CH}}$.

However, there is another, much more relevant, non-affine effect, caused by the saturation of the detectors. These devices are blinded just after an event and might be further kept off to combat the phenomenon of afterpulses, which increases noise. This so-called dead time T_d implies that there is maximum rate of output signals R_sat = 1/T_d that the detectors can produce.

This means that the actual (output) detection rate is not R₀ but [48]

$$\begin{equation}{R}_{\mathrm{det}}=\frac{{R}_{0}\cdot {R}_{\text{sat}}}{{R}_{0}+{R}_{\text{sat}}}.\end{equation} \tag{ 37 }$$

Because this is not an affine function of η_CH, the fluctuations of the latter cannot be neglected. We note that, although equation (37) is derived for a continuous source, it is approximately valid also for a pulsed one if the repetition rate is much greater than R_sat.

To quantify the importance of these fluctuations, we estimate the detection rate of a QKD system for the eight cases of table 1, first considering the entire distributions of η_CH shown in figure 14 and then only their mean values. In figure 15 we show the overestimation factor (converted to dB to visualize it better) caused by neglecting fluctuations. This can reach a value of almost 9 dB for the distribution of case 8. The effect is larger when ⟨R₀⟩ ≈ R_sat and vanishes for ⟨R₀⟩ ≫ R_sat or ⟨R₀⟩ ≪ R_sat, when equation (37) is well approximated by its linearization. The distributions for which this error is greater are those which have stronger tails (i.e., a high kurtosis). Indeed, neglecting fluctuations means neglecting the suppression of the tails caused by the saturation of the detectors: the greater the tails are, the graver is the error caused by neglecting them.

Zoom In Zoom Out Reset image size

In figure 16, we consider an arbitrary QKD scenario with a 1 GHz source, a dead time T_d = 10 μs, and 15% detection efficiency, which are realistic features for a single-photon avalanche diode (SPAD) [49]. We show a simulation of the detection rate (as a function of the channel efficiency η_CH) which neglects fluctuations, and compare it with the more correct values which consider this effect. We can see a clear separation between the two methods of estimation, which grows larger when ⟨R₀⟩ ≈ R_sat and for distributions of greater kurtosis. This shows that performance predictions that consider only the mean value of the channel efficiency can be severely inaccurate.

Zoom In Zoom Out Reset image size

A crucial requirement for the realization of daylight free-space QKD is the successful filtering of the background radiation. Figure 17 shows the spectrum of the diffuse atmospheric radiance I_diff, extracted with LOWTRAN. As we can see, the spectrum peaks at blue wavelengths, and decreases for wavelengths in the infra-red.

Zoom In Zoom Out Reset image size

This datum considers only the sky brightness and neglects the fact that the transmitter partially blocks the FOV of the receiver, therefore it is overestimated in realistic scenarios, especially for short links. Nonetheless, we will use it in the following discussion, which focuses on the impact of background on the performance of QKD.

As a first-order approximation, we can consider the diffuse radiance to be uniform over the receiver FOV. We can estimate the detection rate of background photons per detection window $\left\langle {r}_{\text{sky}}\right\rangle$ at the quantum signal wavelength as a function of the receiver aperture D_Rx, solid-angle FOV Ω, and filtering bandwidth δλ (h is the Planck's constant and c the speed of light)

$$\begin{equation}\left\langle {r}_{\text{sky}}\right\rangle =\frac{{I}_{\text{diff}}\cdot \pi {\left(\frac{{D}_{\text{Rx}}}{2}\right)}^{2}\cdot {\Omega}\cdot \delta \lambda }{hc/\lambda }.\end{equation} \tag{ 38 }$$

For the typically small FOV characteristic of free-space communication systems, Ω can be approximated by

$$\begin{equation}{\Omega}=2\pi \left(1-\mathrm{cos}(\mathrm{F}\mathrm{O}\mathrm{V})\right)\approx \pi {\mathrm{F}\mathrm{O}\mathrm{V}}^{2},\end{equation} \tag{ 39 }$$

where FOV is the one-dimensional field-of-view.

The FOV of the receiver optical system depends on the optical design: if the optical fiber, or free-space detector, is placed on an image-plane of the entrance pupil, which is the case for free-space detectors with large active diameter $(\sim 150\enspace \mu \text{m})$ or large-core multi-mode fibers (MMF), the FOV is essentially a free parameter, limited only by the size of the optical elements (lenses, mirrors) used in the optical system, and may be as large as 400 μrad. In the case of SMFs or free-space detectors with small active diameter $< 10\enspace \mu$m, the optimal choice is to place them on the focal plane of the optical system. In this configuration the FOV is constrained by the size of the active area or fiber MFD. Since the receiver focal length f is chosen so that the optical system efficiency in equation (19) is maximized, we have

$$\begin{equation}\mathrm{M}\mathrm{F}\mathrm{D}={\beta }_{\text{opt}}\frac{4}{\pi }\frac{\lambda f}{{D}_{\text{Rx}}},\end{equation} \tag{ 40 }$$

where β_opt is shown in figure 5 and β_opt = 1.12 for unobstructed apertures. If we define the FOV as the pupil incidence angle at which the spot on the focal plane is deflected to a distance equal to half the MFD, then we have that:

$$\begin{equation}\mathrm{F}\mathrm{O}\mathrm{V}=\frac{\mathrm{M}\mathrm{F}\mathrm{D}}{2f}={\beta }_{\text{opt}}\frac{2}{\pi }\frac{\lambda }{{D}_{\text{Rx}}}.\end{equation} \tag{ 41 }$$

Another consequence of this constraint is that the detection rate of diffuse background photons coupled into the system becomes almost independent of the receiver optical system parameters. Indeed, combining equations (38) and (41), we find that the background light coupled into the SMF is given by

$$\begin{equation}{\left\langle {r}_{\text{sky}}\right\rangle }_{\text{SMF}}={\beta }_{\text{opt}}^{2}\frac{{I}_{\text{diff}}}{hc}{\lambda }^{3}\delta \lambda \end{equation} \tag{ 42 }$$

which does not depend on f nor D_Rx. Figure 18 shows the expected noise count rate for a SMF coupled receiver.

Zoom In Zoom Out Reset image size

Typical DV-QKD systems apply a temporal filter to all detected events, with the purpose of reducing the impact of noise. Indeed, the latter is uniformly distributed in time whereas signal photons, being emitted at regular intervals, have a predictable time of arrival. In post-processing, one can apply a T_gat-wide temporal window centered at this time and discard all events that fall outside of it, thus suppressing noise by a factor T_gat ⋅ R. However, there is a tradeoff, because the temporal distribution of the signal events is enlarged by the optical pulse width, by the jitter of the source, of the detectors, and of the time-digitizing hardware. A small value of T_gat, while strongly reducing noise, might discard too much of the signal, negatively impacting the final SKR. Assuming a normal distribution of standard deviation J for the signal time of arrival, the filter reduces the signal detection rate by a factor $\mathrm{e}\mathrm{r}\mathrm{f}\left(\frac{{T}_{\text{gat}}}{J2\sqrt{2}}\right)$.

A numerical study of the tradeoff can guide the choice of T_gat. The figure of merit to maximize is the final SKR, which includes the contribution of the detection and error rates. We focus on the ratio T_gat/J between it and the standard deviation J of the temporal distribution of the signal (including all the aforementioned jitter contributions). We can expect the tradeoff to be influenced by (i) the signal-to-noise ratio (SNR) when the noise is gated to a window as wide as the signal (±3 times the standard deviation J) and (ii) the coding error.

The SNR obtained before any gating is not sufficient to describe the situation, because it does not consider the width of the signal. Intuitively, for the same ungated SNR, a temporally wider signal favors smaller values of T_gat/J to eliminate more noise. Our definition of the SNR, by considering only the portion of noise that falls under the signal, effectively combines the ungated SNR and the width of the signal in a single parameter. Other protocol parameters such as decoy intensity levels and probabilities also have an importance, but we focus only on the two quantities above for simplicity. All the parameters of the model except the SNR and coding error are arbitrarily fixed to a realistic QKD scenario.

In figure 19, we can see the results of an optimization with the Nelder–Mead algorithm [52] of the T_gat/J ratio to maximize the SKR, for a grid of values of the SNR (noise gated at ±3J) and coding error. Predictably, the smaller the SNR, the smaller the gating window should be, in order to discard more noise. High values of the coding error decrease the optimal T_gat/J. Indeed, although temporal gating alone cannot change the coding error, the higher the total error contribution, the more important it is to reduce it, even at the expense of discarding part of the signal. This shifts the balance of the tradeoff toward smaller windows, as can be seen in the top-left corner of the figure.

Zoom In Zoom Out Reset image size

The small ridge in figure 19 indicates a sudden jump in the optimal value of T_gat/J. This is because there are several terms contributing to the SKR (see (44)), and each is computed with several methods [27], choosing the best for every configuration of the parameters. This leads to the presence of multiple local maxima: when one of them is promoted to global maximum, overcoming another, the optimal T_gat/J changes. The exact position and entity of the ridge is not universal and depends on the specific scenario that we chose in this simulation, but its presence is to be expected in general when the SKR is optimized.

Because actual experiments accumulate only a finite amount of data, the parameters mentioned above cannot be estimated with perfect precision, leaving an opening for potential attackers. To counter this, QKD uses a broad range of statistical analyses [53–55]. Following reference [27], we construct confidence intervals based on the Hoeffding inequality [56] around the parameters of interest, and then use the pessimistic extrema of the intervals as our estimates. This penalizes the performance of the system, but guarantees that the key is secure with a very high probability 1 − _sec.

Fortunately, we can mitigate this cost by optimizing some parameters of the protocol. These are: the probability p_Z that each of Alice and Bob choose the key basis Z for their preparations and measurements, the probability p_μ that Alice chooses the stronger intensity level, and the intensities μ, ν of each level. Their optimal values that maximize the SKR change depending on the amount of data (block length) that is used in the statistical analysis. In our study, we consider two example QKD scenarios: a high-end system with superconducting nanowire single-photon detectors, a GHz source, and low coding error (scenario A), and a less expensive one with SPADs, a slower source, and a higher coding error. The detector parameters are arbitrary, but similar to those of commercial devices [49, 57]. The most relevant parameters of each scenario are listed in table 2.

In figure 20 we show the results of the optimization (with a simulated annealing procedure [58]) when the channel is characterized by the distribution of case 1 (figure 14), with ⟨η_CH⟩ ≈ −7 dB. We can see an upward trend for p_Z and p_μ for growing block length, indeed the larger the total sample size, the easier it is to accumulate the needed statistics in the check-basis basis X (mutually unbiased with respect to Z) and intensity level ν even if they are chosen rarely. The mean photon numbers μ and ν show more stability, with a slight downward trend for ν in the right-hand side of the plot. This is justified considering that its limit for infinite block length is zero, as this makes the bounds of the decoy-state method tighter. For shorter block lengths, ν must grow to increase the detection rate and accumulate the needed statistics. The optimal value of μ is the one that strikes the right balance between a low multi-photon emission probability and a high detection rate. However, it is also influenced by other factors like afterpulses and detector saturation, which is why it is smaller in scenario B. For low and decreasing block length, as p_μ shrinks to increase the statistics accumulated in the low intensity level, μ grows to keep the signal rate high, and ν responds by slightly decreasing to tighten the bounds.

Zoom In Zoom Out Reset image size

The jump which can be seen in figure 20(b) is similar to the ridge of figure 19: a local maximum is promoted to global and the optimal parameters suddenly move. The exact position and entity of the jump depends strongly on the scenario, but its presence is to be expected in optimizations like these.

In figure 21 we show the cost of finite-key effects and of using a wrong set of parameters. We express it as the ratio between the SKR and SKR_∞, i.e., the SKR that would be obtained by optimizing the above parameters for infinite length and by accumulating an infinitely long key block. The dashed line corresponds to the choice of optimizing the parameters for a fixed block length of 10⁷ bits. Generating a key is possible but the SKR quickly drops to zero if block lengths shorter than 10⁷ are used. This happens more slowly in scenario B because the noisier detectors reduce the impact of tuning signal-related parameters.

Zoom In Zoom Out Reset image size

The solid line shows what happens if the parameters are optimized for each block length. The results improve drastically, underlining the importance of optimizing the protocol parameters for the predicted size that will be used operatively. However, even with a large block length of 10⁹ bits, the SKR is reduced by a sizeable portion with respect to SKR_∞. While the precise value of this portion depends on the chosen scenarios, the fact that finite-key effects should not be neglected even for large block lengths is general.

To conclude our analysis and give an example of the capabilities of our full model, we calculate the final SKR of a QKD system considering all the effects we have studied in this work. We do this for the two scenarios of table 2 and the eight channel efficiency distributions of figure 14. As explained in section 2.5, we calculate R_det using equations (36) and (37) for multiple values of the fluctuating η_CH, and then find its weighted average over the corresponding distributions. For each configuration, we use the simulated annealing algorithm to optimize the protocol parameters (the same of section 4.1) and the temporal gating, for a fixed block length of 10⁷ bits.

The results are shown in figure 22. We can observe the characteristic linear behavior of the SKR with transmittance, and glimpse a drop for strong losses caused by the prevalence of noise. However, thanks to the breadth of phenomena that our model includes, these values go beyond the simple verification of this typical trend, and are accurate estimates of the performance of the considered QKD systems. Scenario B is strongly penalized by its slower source, lower detection efficiency, afterpulses and dead time. Because of this, only the first four distributions yield a positive SKR.

Zoom In Zoom Out Reset image size