Classical-quantum dual encoding for laser communications in space

For the purpose of classical optical communications, the Intensity Modulation (IM) technique (which maps a classical information bit to one of the output intensities of a laser beam), together with the Direct Detection (DD) technique (which directly measures the intensity of the received signal using a photodetector), constitutes the widely-used IM-DD scheme. Despite its low implementational complexity, the IM-DD scheme is not desirable for our purpose mainly due to the IM technique's incompatibility with the quantum part of our protocol.

Instead, in this work we adopt a discrete modulation scheme where every classical information bit is mapped to one of the two real-valued displacements $\{\alpha,-\alpha\}$ applied to a continuous-wave single-mode laser beam with optical frequency Ω₀ and overall amplitude $\beta \gt\gt |\alpha|$. This scheme resembles the Binary Phase-Shift Keying (BPSK) modulation scheme since it only encodes information onto one quadrature (i.e. the X quadrature) of the laser beam. Our modulation at a specific sinusoidal modulation frequency, ω, creates an upper sideband and a lower sideband at $\Omega_0+\omega$ and $\Omega_0-\omega$, respectively. The component remaining at Ω₀ is known as the carrier, which can be considered as a strong coherent component [34, 35]. We will show in the following that although our protocol does not use IM for modulation, it can still benefit from a reduced implementational complexity provided by the use of DD for detection.

Traditionally, the extraction of quadrature-encoded information requires a coherent homodyne detection involving the mixing of the signal with an independent Local Oscillator (LO). By independent LO, we mean a bright laser beam that is phase-locked to the signal and typically sent along with the signal via multiplexing. In this work, we eliminate the need for such an independent LO by adopting the DD technique, which effectively measures the X quadrature (i.e. the only quadrature used for information encoding in our scheme) using a photodetector, for signal detection. To see how the DD technique works, we first consider a propagating continuous-wave laser beam in a single mode represented by the continuum operator $\hat{A}$:

$$\begin{align} \hat{A} & = \beta + \delta\hat{A}\left(t\right), \end{align} \tag{ 1 }$$

where β is the mean value of the carrier amplitude (taken to be real) and $\delta\hat{A}(t)$ represents the modulations (including both the discrete classical modulation and the quantum fluctuations) around the mean value. The notation δ emphasises that the modulations are small compared with β (note that $\langle \hat{A} \rangle = \beta$). The operation of DD using a photodetector is represented by the number operator, $\hat{N}$. Linearising the operator by assuming $\beta \gg 1$, we have

$$\begin{equation} \begin{aligned} \hat{{N}} & = \hat{{A}}^{\dagger}\left(t\right) \hat{{A}}\left(t\right) \\ & = \left(\beta+\delta \hat{{A}}^{\dagger}\left(t\right)\right)\left(\beta+\delta \hat{{A}}\left(t\right)\right) \\ & = \beta^2+\beta \delta \hat{{A}}^{\dagger}\left(t\right)+\beta \delta \hat{{A}}\left(t\right)+\delta \hat{{A}}^{\dagger}\left(t\right) \delta \hat{{A}}\left(t\right) \\ & = \beta^2+\beta \delta \hat{{X}}\left(t\right) \end{aligned} \end{equation} \tag{ 2 }$$

where $\delta \hat{X}(t) = \delta \hat{A}^\dagger (t) + \delta \hat{A}(t)$ is the X quadrature operator, and we have neglected the higher-than-first-order terms in the modulations due to the assumption that the fluctuations are small (i.e. the contribution to the expectation value from the higher order terms is negligible). The last line of equation (2) indicates that the use of DD effectively measures the X quadrature. Equation (2) also explains the counter-intuitive concept of using DD to measure the X quadrature—in the same beam, the small signals (i.e. the modulation-induced fluctuations) beat with the carrier (i.e. the strong coherent component) and can be read out at a particular sideband frequency ω. The technique for side-band readout proceeds as it would for homodyne detection and is described in, e.g. [36, 37].

According to our Classical-Quantum Dual Encoding protocol, we are interested in the modulations around side-band frequency ω, hence we focus our attention on the evolution of the Fourier space operator $\hat X_{\bar A}^{^{\prime}}(\omega)$ which is the Fourier transform of $\delta \hat{X}_{\bar A}^{^{\prime}}(t)$, representing the modulations of the field leaving Alice's station. We can write this Fourier space operator as:

$$\begin{align} \hat{X}_{\bar A}^{^{\prime}} & = X_C + \hat{X}_{\bar A}, \end{align} \tag{ 3 }$$

where $\langle X_C \rangle = \pm 2 \alpha$ is our discrete classical modulation, and $\hat{X}_{\bar A}$ represents quantum fluctuations. After traversing the channel to Bob this operator evolves to

$$\begin{align} \hat{X}_B^{^{\prime}} & = \sqrt{\eta} \hat{X}_{\bar A}^{^{\prime}} + \sqrt{1- \eta} \hat{X}_v \nonumber \\ & = \sqrt{\eta}\left(X_C + \hat{X}_{\bar A}\right) + \sqrt{1- \eta} \hat{X}_v , \end{align} \tag{ 4 }$$

where η is the transmission of the channel and $\hat{X}_v$ represents vacuum fluctuations that couple in due to the loss. As will be described shortly it is arranged that the variance of the quantum fluctuations of the beam leaving Alice's station are at the QNL, i.e. $\langle \hat{X}_{\bar A}^2 \rangle = V_{\bar A} = 1$. As a result the signal to noise ratio (i.e. the signal power divided by the noise variance) of the classical signal received by Bob will be $SNR = \eta \langle X_C^2 \rangle/(\eta \langle \hat{X}_{\bar A}^2 \rangle+$ $(1- \eta) \langle \hat{X}_v^2 \rangle) = 4 \eta \alpha^2$, where we have used that the vacuum fluctuations also have unit variance. By choosing α sufficiently large the resulting bit error rate will be small and faithful transmission of the classical signals can be achieved.

To generate a key, we need to extract the quantum fluctuations from the combined classical and quantum signal that has been mixed down at our chosen side-band frequency ω. Notice that the carrier power detected by Bob over the signal length will be proportional to $\eta \beta^2$ where η is the channel transmission for that particular time-bin (assumed constant over the signal length). Thus η can be measured shot-by-shot. Also, provided α is sufficiently large, we know X_C shot-by-shot from the sign of the signal displacement. Thus we can extract the quantum signal by subtracting the weighted classical signal via

$$\begin{align} \hat{X}_B & = \sqrt{\eta} \hat{X}_{\bar A}^{^{\prime}} + \sqrt{1- \eta} \hat{X}_v - \sqrt{\eta} X_C\nonumber \\ & = \sqrt{\eta} \hat{X}_{\bar A} + \sqrt{1- \eta} \hat{X}_v, \end{align} \tag{ 5 }$$

where we have used equation (3) in going from the first to the second line. Note that the shot-by-shot carrier power also enables us to calibrate the QNL level for each signal. The time bin length should be chosen appropriately; short enough so that the fluctuating power due to the atmosphere is approximately constant over the time bin but long enough to acquire sufficient signal strength.

We now analyse the generation and evolution of the quantum signal as shown diagrammatically in figure 1. As shown, Alice begins with two phase-locked beams, one of which has super-Poissonian amplitude noise, $\langle \hat{X}_a^2 \rangle = V_a \gt 1$, and the other has sub-Poissonian amplitude noise, $\langle \hat{X}_s^2 \rangle = V_s \lt 1$. One of the beams has a bright carrier—it does not matter which but we will assume here that it is the super-Poissonian beam. The sub-Poissonian, or squeezed beam is intrinsically quantum, however the noise on the super-Poissonian beam could either be quantum, e.g. anti-squeezing or amplified spontaneous emission, or it could be induced via classical modulations. The most important thing is that the noise has a Gaussian distribution. Alice mixes her two beams on a beam-splitter, retaining the reflected beam which she directly measures and sending the transmitted beam (after imposing the classical signal) through the channel to Bob. Eve intercepts all the lost light. We can then write the evolution of the quantum parts of the amplitude quadrature operators as

$$\begin{align} \hat X_{A} & = \sqrt{1-\varepsilon} \hat X_a+\sqrt{\varepsilon}\hat X_s \end{align} \tag{ 6 }$$

$$\begin{align} \hat X_E & = \sqrt{1-\eta} \left(\sqrt{\varepsilon} \hat X_a + \sqrt{1-\varepsilon} \hat X_s \right) + \sqrt{\eta} \hat X_v \end{align} \tag{ 7 }$$

$$\begin{align} \hat X_B & = \sqrt{\eta} \left(\sqrt{\varepsilon} \hat X_a + \sqrt{1-\varepsilon} \hat X_s \right) + \sqrt{1-\eta} \hat X_v, \end{align} \tag{ 8 }$$

where ε is the transmissivity of Alice's beamsplitter, η is the transmissivity of the channel, $\hat X_{A}$, $\hat X_B$, and $\hat X_E$ are the X-quadrature operators for Alice, Bob, and Eve, respectively, and $\hat X_s$ and $\hat X_a$ are the X-quadrature operators for the squeezed (sub-Poissonian) and antisqueezed (super-Poissonian) beams, respectively.

For zero information leakage we require that the modulation Alice prepares to send to Bob has the same variance as the vacuum state (i.e. the QNL). That is, we require

$$\begin{align} V_{\bar A} & = \varepsilon V_a+\left(1-\varepsilon\right) V_s = 1. \end{align} \tag{ 9 }$$

For the modulation to be QNL we need

$$\begin{align} \varepsilon & = \frac{1-V_s}{V_a-V_s},\; 1-\varepsilon = \frac{V_a-1}{V_a-V_s}. \end{align} \tag{ 10 }$$

Then,

$$\begin{align} \langle \hat X_A \hat X_B \rangle & = \sqrt{\eta}\sqrt{V_a-1}\sqrt{1-V_s} \end{align} \tag{ 11 }$$

$$\begin{align} \langle \hat X_E \hat X_B \rangle & = 0. \end{align} \tag{ 12 }$$

That is, there are correlations in the q quadrature between Alice and Bob, but there are no correlations in the q quadrature between Eve and Bob. This means the Holevo information $\chi_{EB} = 0$. The Holevo information, χ_EB , represents the maximum information that Eve can extract from her signal about Bob's measurements. Since $\langle \hat X_E \hat X_B \rangle = 0$, Eve's system is not correlated with Bob's system in the X quadrature.

Given no information leakage to Eve about Bob's measurement results we will find we only need Alice and Bob's covariance matrix. Alice and Bob's covariance matrix is as follows:

$$\begin{align} \Gamma_{AB} & = \left[ \begin{matrix} \langle \hat X_A^2 \rangle & \langle \hat X_A \hat P_A \rangle & \langle \hat X_A \hat X_B \rangle & \langle \hat X_A \hat P_B \rangle \\ \langle \hat P_A \hat X_A \rangle & \langle \hat P_A^2 \rangle & \langle \hat P_A \hat X_B \rangle & \langle \hat P_A \hat P_B \rangle \\ \langle \hat X_B \hat X_A \rangle & \langle \hat X_B \hat P_A \rangle & \langle \hat X_B^2 \rangle & \langle \hat X_B \hat P_B \rangle \\ \langle \hat P_B \hat X_A \rangle & \langle \hat P_B \hat P_A \rangle & \langle \hat P_B \hat X_B \rangle & \langle \hat P_B^2 \rangle \end{matrix} \right] \nonumber \\ & = \left[ \begin{matrix} a_q & 0 & c_q & 0 \\ 0 & a_p & 0 & c_p \\ c_q & 0 & b_q & 0 \\ 0 & c_p & 0 & b_p \end{matrix} \right], \end{align} \tag{ 13 }$$

where

$$\begin{align} a_q & = V_s\varepsilon + V_a\left(1{-}\varepsilon\right) \end{align} \tag{ 14 }$$

$$\begin{align} a_p & = \frac{\varepsilon}{V_s} + \frac{1{-}\varepsilon}{V_a} \end{align} \tag{ 15 }$$

$$\begin{align} b_q & = \eta_f\left[V_a\varepsilon + V_s\left(1{-}\varepsilon\right)\right] - \eta_f + 1 + \text{Var}\left(\sqrt{\eta}\right)\left(V_a\varepsilon + V_s\left(1{-}\varepsilon\right) -1\right) \end{align} \tag{ 16 }$$

$$\begin{align} b_p & = \eta_f\left( \frac{\varepsilon}{V_a} + \frac{1{-}\varepsilon}{V_s}\right) - \eta_f+ 1 + \text{Var}\left(\sqrt{\eta}\right)\left( \frac{\varepsilon}{V_a} + \frac{1{-}\varepsilon}{V_s} - 1 \right) \end{align} \tag{ 17 }$$

$$\begin{align} c_q & = \sqrt{\eta_f}\sqrt{\varepsilon}\sqrt{1{-}\varepsilon}\left(V_a{-} V_s \right) \end{align} \tag{ 18 }$$

$$\begin{align} c_p & = \sqrt{\eta_f}\sqrt{\varepsilon}\sqrt{1{-}\varepsilon}\left( \frac{1}{V_a} {-} \frac{1}{V_s}\right). \end{align} \tag{ 19 }$$

Here we have allowed for a fluctuating channel transmission as expected for an atmospheric channel, where $\eta_f = \langle \sqrt{\eta} \rangle ^2$ and Var$(\sqrt{\eta}) = \langle {\eta} \rangle - \eta_f$, and where $\langle \eta \rangle = \int_0^{\eta_\text{max}} \eta p(\eta) \; \text{d}\eta$ and $\langle \sqrt{\eta} \rangle = \int_0^{\eta_\text{max}} \sqrt{\eta} p(\eta) \; \text{d}\eta$. Notice that given the no information leakage condition (equation (10)). we find $b_q = 1$, so the fluctuating channel transmission does not add noise to the quantum amplitude correlations.

The secret key protocol we propose is an adaption of CV-QKD with reverse reconciliation [33], using the no leakage protocol [31] and a passive eavesdropper [26]. In the reverse reconciliation protocol Alice guesses the value of Bob's quadrature measurement outcomes based on her own outcomes. They then perform a reconciliation protocol which corrects errors in Alice's guesses at the cost of shortening the data string. In the limit of an infinitely long string this process does not leak any information to the eavesdropper—this is called the asymptotic limit. The asymptotic secret key rate with reverse reconciliation is [27]

$$\begin{align} K_{\text{asym.}} & = \beta I_{AB} - \chi_{EB}, \end{align} \tag{ 20 }$$

where β is the reverse reconciliation efficiency, I_AB is Alice and Bob's mutual information, and χ_EB is a bound on Eve's maximal information given by the Holevo Bound. For the no leakage protocol with a passive eavesdropper we have $\chi_{EB} = 0$ and so $K_{\text{asym.}} = \beta I_{AB}$.

The best possible performance of our key protocol will be in this asymptotic limit, with ideal reconciliation (β = 1), strong squeezing ($V_s \lt\lt 1$) and ideal detectors. Under these conditions we have $I_{AB} = \frac{1}{2}\log_2\frac{a_q}{a_q-{c_q^2}/{b_q}}$. Substituting the values in this limit we find the highest achievable rate of the zero-leakage protocol is

$$\begin{align} K_{\text{ideal}} & = -\frac{1}{2}\log_2{\left(1-\eta_f\right)}. \end{align} \tag{ 21 }$$

In order to include the effects of finite data block lengths a more sophisticated analysis is required, as will be introduced in section 4.

As a specific example, we consider the satellite-to-Earth communication setup shown in figure 2—it should be noted that our system is not restricted to such an arrangement. In this work, the ground-station altitude is denoted as h₀, the satellite zenith angle at the ground station is denoted as θ_z, and the satellite altitude at $\theta_{\text{z}} = 0$ is denoted as H. The channel distance can then be given by $L = (H-h_0)/\text{cos}\theta_{\text{z}}$. We assume a restricted eavesdropping scenario where the channel is monitored [25], and hence Eve is passive and not able to perform an active attack [38]. However, as indicated in figure 2, we allow Eve to capture all the lost light.

Zoom In Zoom Out Reset image size

The random inhomogeneities (i.e. turbulent eddies) within the Earth's turbulent atmosphere (which can be seen as a random medium) give rise to small random refractive index fluctuations (i.e. optical turbulence) that cause continuous phase modulations on an optical beam, imposing amplitude and phase distortions on the optical beam as it propagates through an atmospheric channel. Specifically, the family of turbulent eddies whose sizes are bounded above by the outer scale L₀ and below by the inner scale l₀ constitutes the inertial subrange [39]. It is commonly assumed that atmospheric turbulence satisfies

$$\begin{equation} \langle n\left(\mathbf{R}\right) \rangle\! = \!1, \,\, \delta n\left(\mathbf{R}\right)\!\ll\!1, \,\, \langle\delta n\left(\mathbf{R}\right)\rangle = 0, \end{equation} \tag{ 22 }$$

where $\mathbf{R} = [x,y,z]^{T}$ is the three-dimensional position vector, $n(\mathbf{R})$ denotes the refractive index at R, and $\delta n(\mathbf{R})\! = \!n(\mathbf{R}) - \langle n(\mathbf{R})\rangle$ denotes the small refractive index fluctuations [39]. Under the paraxial approximation and the assumptions in equation (22), the propagation of a monochromatic optical beam $\psi(\mathbf{R})$ through the turbulent atmosphere is governed by the stochastic parabolic equation [39]

$$\begin{equation} \nabla_{\text{T}}^{2} \psi\left(\mathbf{R}\right)+i 2 k \frac{\partial\psi\left(\mathbf{R}\right)}{\partial {z}} +2 \delta n\left(\mathbf{R}\right) k^{2} \psi\left(\mathbf{R}\right) = 0, \end{equation} \tag{ 23 }$$

where $\nabla_{\text{T}}^{2} = \partial^{2} / \partial x^{2}+\partial^{2} / \partial y^{2}$ is the transverse Laplacian operator.

In this work, we assume that all the transmitted optical beams are in the fundamental Gaussian mode. The general form of a Gaussian beam propagating along the z axis is given by

$$\begin{equation} \begin{aligned} \psi_{\text{G}}\left(\mathbf{r},z\right)& = \sqrt{\frac{2}{\pi}} \frac {1 } {w \left( z \right) } \text{exp} \left[ \frac {- |\mathbf{r}| ^ {2 } } {w ^ {2 } \left( z \right) } \right]\text{exp} \left[i\phi\left(\mathbf{r},z\right)\right], \end{aligned} \end{equation} \tag{ 24 }$$

where $\mathbf{r} = [x,y]^T$ denotes the two-dimensional position vector in the transverse plane, $w(z) = w_{0} \sqrt{1+(z / z_{\text{R}})^{2}}$, w₀ is the beam-waist radius, $z_{\mathrm{R}} = \pi w_{0}^{2}/\lambda$ is the Rayleigh range, λ is the (central) optical wavelength, and $\phi(\mathbf{r},z) = \frac { k |\mathbf{r}| ^ {2 } z } {2 \left( z ^ {2 } + z _ {\text{R} } ^ {2 } \right) } - \text{arctan} ( \frac {z } {z _ {\text{R} } } )$ with $k = 2\pi/\lambda$ being the optical wavenumber.

Although our scheme is, in principle, robust against the random phase distortions (recall that no local oscillator is required in our direct detection setup—the wavefront matching between the local oscillator and the signal beam is not required), turbulence-induced effects related to amplitude distortions can still negatively affect the performance of our scheme. Indeed, these effects (e.g. beam wandering, beam-shape deformation, beam broadening, and intensity scintillation) can introduce substantial losses in real-world scenarios where a finite-sized aperture is used at the receiver [39]. Due to the random nature of atmospheric turbulence, such losses are also randomly fluctuating, giving rise to non-Gaussian excess noise and thus further degrading the system performance (see discussions in e.g. [40, 41]).

We adopt the standard notation 'loss' for visualization purposes; however, in our calculations, we always use the transmissivity η (with the relation $\text{Loss}\,\text{[dB]} = -10 \log_{10} \eta$) to maintain consistency with our analyses in the previous sections. Assuming the transmitted beam is denoted by $\psi_{\text{G}}(\mathbf{r},0)$ and the received beam is denoted by $\psi_{\text{Rx}}(\mathbf{r},L)$, the (fluctuating) transmissivity is given by

$$\begin{equation} \eta = \int_{\mathcal{A}} \left|\psi_{\text{Rx}}\left(\mathbf{r},L\right)\right |^2\,\text{d}^2\mathbf{r}, \end{equation} \tag{ 25 }$$

where L denotes the channel distance, and $\mathcal{A}$ denotes the aperture area at the receiver. Note that we assume a circular receiver aperture whose radius is r_a in this work.

3.2.1. Turbulence characterization

The optical turbulence within the Earth's atmosphere must first be appropriately characterized. In this work, we adopt the widely used Hufnagel–Valley (HV) model to determine the structure parameter $C_n^2(h)$ that describes the strength of the optical turbulence within a satellite-based atmospheric channel. The HV model expresses $C_n^2(h)$ as the following function of altitude h [39].

$$\begin{align} C_{n}^{2}\left(h\right) = 0.00594\left(v_{\text{rms}}/27\right)^{2}\left(h \times 10^{-5}\right)^{10} \text{exp}{\left(-h/1000\right)}+2.7\!\times\!10^{-16} \text{exp}{\left(-h/1500\right)} +A \text{exp}{\left(-h/100\right)}, \end{align} \tag{ 26 }$$

where A is the ground-level turbulence strength in $\text{m}^{-2 / 3}$, and v_rms is the root-mean-square (rms) wind speed in m s⁻¹. The rms wind speed is given by

$$\begin{equation} v_{\text{rms}} = \left[\frac{1}{15 \times 10^{3}} \int_{5 \times 10^{3}}^{20 \times 10^{3}} V^{2}\left(h\right) \,\text{d}h\right]^{1 / 2}, \end{equation} \tag{ 27 }$$

where V(h) is the altitude-dependent wind speed profile. In this paper we adopt the Bufton wind profile [39]

$$\begin{equation} V\left(h\right) = V_{\text{g}}+30 \text{exp} \left[-\left(\frac{h-9400}{4800}\right)^{2}\right], \end{equation} \tag{ 28 }$$

where V_g is the ground-level wind speed.

Utilizing the structure parameter profile of the Earth's atmosphere, we adopt the scintillation index $\sigma_{\text{I}}^2$ (which is the normalized variance of intensity) and the Fried parameter r₀ (which quantifies the coherence length of turbulence-induced phase errors in the transverse plane) to quantify the effect of atmospheric turbulence on a propagating beam. For satellite-to-Earth channels, the scintillation index is given by [39]

$$\begin{equation} \sigma_{\text{I}}^{2} = \text{exp}\!\left[\frac{0.49 \sigma_{\text{R}}^{2}}{\left(1+1.11 \sigma_{\text{R}}^{12 / 5}\right)^{7 / 6}}\!+\!\frac{0.51 \sigma_{\text{R}}^{2}}{\left(1+0.69 \sigma_{\text{R}}^{12 / 5}\right)^{5 / 6}}\right]-1, \end{equation} \tag{ 29 }$$

with $\sigma_{\text{R}}^2$ being the Rytov variance,

$$\begin{equation} \sigma_{\text{R}}^{2} = 2.25 k^{7 / 6} \text{sec} ^{11 / 6}\left(\theta_{\text{z}}\right) \int_{h_{0}}^{H} C_{n}^{2}\left(h\right)\left(h-h_{0}\right)^{5 / 6} \,\text{d} h. \end{equation} \tag{ 30 }$$

For satellite-to-Earth channels, the Fried parameter is given by [39]

$$\begin{equation} r_0 = \left[0.423 k^{2} \text{sec} \theta_{\text{z}} \int_{h_0}^{H} C_{n}^{2}\left(h\right) \,\text{d} h\right]^{-3 / 5}. \end{equation} \tag{ 31 }$$

The rate at which atmospheric turbulence changes with time is commonly quantified by the Greenwood frequency [42, 43]

$$\begin{equation} f_{\text{G}} = 2.31 \lambda^{-6 / 5}\left[\text{sec} \theta_{\text{z}} \int C_n^2\left(h\right) V^{5 / 3}\left(h\right)\,\mathrm{d} h\right]^{3 / 5}. \end{equation} \tag{ 32 }$$

Correspondingly, the time interval over which atmospheric turbulence remains essentially unchanged can be quantified by the atmospheric coherence time [43]

$$\begin{equation} \tau_0 = \frac{0.134} {f_{\text{G}}}. \end{equation} \tag{ 33 }$$

3.2.2. Method of numerical simulations

In this work, the implementation of our numerical simulation largely follows our previous work [44] (see, e.g. the appendix of [44] for more discussions). Specifically, we numerically solve equation (23) using the split-step method [45] (also referred to as phase screen simulations) in order to faithfully capture the effects imposed by the optical turbulence on an optical beam propagating within the Earth's atmosphere. The main idea of this method is to model an atmospheric channel using a set of slabs with a random phase screen located in the midway of each slab. Two vacuum propagations with one random phase modulation in between are repeatedly performed until the beam reaches the receiver.

The implementation of our phase screen simulations requires careful configuration and optimization based on the evaluation of the scintillation index (recall equation (29)) and the Fried parameter (recall equation (31)) for both the whole channel and each individual slab. In order to characterize the atmospheric turbulence within the jth slab, both the scintillation index (denoted as $\sigma_{\text{I}_j}^2$) and the Fried parameters (denoted as $r_{0_j}$) have to be evaluated locally for the turbulent atmosphere within that slab. For the actuate modelling of the atmospheric channel, we set the width of each slab based on the two conditions (i.e. $\sigma_{\text{I}_j}^2\lt 0.1$ and $\sigma_{\text{I}_j}^2\lt 0.1\sigma_\text{I}^2$) described in [46].

After determining the widths of the slabs, we use a well-known Fast-Fourier-Transform (FFT)-based spectral-domain algorithm [47] to generate (random realizations of) the phase screen for each slab using its corresponding phase Power Spectral Density (PSD) function. In this algorithm, the generation of each individual phase screen starts with the spectral-domain generation of a uniform square grid of random complex numbers sampled from a complex Gaussian distribution whose real and imaginary parts each have zero mean and equal variances with zero cross-covariances. This grid of random numbers is further weighted by the corresponding phase PSD function of the atmospheric turbulence within the corresponding slab. A FFT is then performed to transform this spectral-domain grid of weighted random complex numbers into a spatial-domain grid of random phase values (corresponding to an array of sample points of the phase screen) with the same statistics as the turbulence-induced phase variations within the corresponding slab. In this work, we adopt the widely-used modified von Karman model for the atmospheric turbulence, giving the following phase PSD function for the jth slab

$$\begin{equation} \Phi_{\phi_j}^{\mathrm{mvK}}\left(f\,\right) = 0.023 r_{0_j}^{-5 / 3} \frac{\text{exp} \left(-f^{\,\,2} / f_{\text{m}}^{\,\,2}\right)}{\left(f^{\,\,2}+f_0^{\,\,2}\right)^{11 / 6}}, \end{equation} \tag{ 34 }$$

where f is the magnitude of the two-dimensional spatial frequency vector in the transverse plane in cycles/m, $f_0 = 1/L_{\text{outer}}$, and $f_{\text{m}} = 0.9422/l_\text{inner}$ [45].

For the free-space propagation, we utilize a physical optical propagation library named proper, whose routines implement FFT-based angular-spectrum and Fresnel-approximation methods to propagate a wavefront in near-field and far-field conditions, respectively [48]. As a fully numerical approach to solving the stochastic parabolic equation, despite its computationally intensive nature, the split-step method can take into account most of the important effects (e.g. beam wandering, beam broadening, and scintillation) in atmospheric optical propagation. Although this method can also be adopted to characterize the temporal profile of the turbulence (i.e. the temporal turbulent fluctuations within one coherence time τ₀), we notice that such a temporal profile characterization is not required within the scope of this work since an atmospheric channel can be considered constant for a large number of sequential signal pulses (more discussions will be provided in section 3.4). The split-step method has been widely used to study the problem of atmospheric optical propagation under various channel conditions, providing quantitative agreement with both analytical results (e.g. [49, 50]) and real-world experimental data (e.g. [51]). More technical details regarding the implementation of our numerical simulations can be found in [44, 50–52] and the references therein.

We restrict ourselves to the case of a low-Earth-orbit (LEO) satellite with a satellite altitude $H = 500\,\text{km}$, and we consider zenith angles θ_z ranging from 0^∘ to 60^∘. For the atmospheric parameters, we follow previous works (e.g. [44, 52]) and set $A = 9.6\times 10^{-14}\, \text{m}^{-2 / 3}$, $L_{\text{outer}} = 5\,\text{m}$, and $l_{\text{inner}} = 1\,\text{cm}$. For the wind speed profile, we set $V_{\text{g}} = 3\,\mathrm{m\,s}^{-1}$, giving a value of $v_{\text{rms}} = 21\,\mathrm{m\,s}^{-1}$. For the optical parameters, we set $\lambda = 1064\,\text{nm}$, and we set the beam-waist radius to $w_0 = 15\,\text{cm}$. Particularly, our setting of the $1064\,\text{nm}$ wavelength is due to the existence of high-quality sources of squeezed light at this wavelength (see, e.g. [53]). This wavelength also falls within the atmospheric transmission window and is currently one of the most widely used operating wavelengths in real-world space-based optical communication systems (see comprehensive review papers, e.g. [54]). For real-world implementations, we consider three reasonable receiver aperture radii, namely, $r_{\text{a}} = 15\,\text{cm}$, $r_{\text{a}} = 30\,\text{cm}$, and $r_{\text{a}} = 50\,\text{cm}$, at the ground station. We further fix the ground-station altitude to $h_0 = 0\,\text{m}$. In our numerical simulations, we generate 10 000 independent realizations of the satellite-to-Earth channel under each simulation setting. Our simulation parameters are listed in table 1.

Table 1. Simulation parameters.

A	Louter	linner	Vg	vrms	λ	w0	h0	H	θz	ra
$9.6\times 10^{-14}\,\text{m}^{-2/3}$	$5\,\text{m}$	$1\,\text{cm}$	$3\,\mathrm{m\,s}^{-1}$	$21\,\mathrm{m\,s}^{-1}$	$1064\,\text{nm}$	$15\,\text{cm}$	$0\,\text{m}$	$500\,\text{km}$	0∘ to 60∘	$\{15\,\text{cm},\,30\,\text{cm},\,50\,\text{cm}\}$

Although the existence of pointing errors can lead to an additional loss penalty in a real-world system, the resulting loss penalty can be significantly reduced by the use of state-of-the-art tracking techniques. Indeed, via the use of feedback-based tracking techniques, the loss penalty resulting from pointing errors over a $1200\,\text{km}$ satellite-to-Earth channel is estimated to be smaller than $3\,\text{dB}$ (which is lower than the turbulence-induced loss and much lower than the diffraction loss) in a well-known real-world demonstration of satellite-to-Earth QKD [55]. Therefore, in this work we ignore the pointing-error-induced loss as we expect that including such a loss will not significantly change any of our observations.

Before presenting our main results, we first provide a visualization of the channel loss statistics (predicted from our numerical simulations) over a satellite-to-Earth channel.

To illustrate the time-dependent fluctuation of channel loss, in figure 3 we first plot the fluctuating loss over a satellite-to-Earth channel ($\theta_{\text{z}} = 60^{\circ}$ and $r_{\text{a}} = 15\,\text{cm}$) within a $0.5\,\text{s}$ time period. It should be noted that figure 3 presents the channel loss fluctuation as a step function, with step size being the atmospheric coherence time τ₀—this is because the channel loss can be considered constant within the time interval specified by τ₀. Specifically, the atmospheric coherence time under this setting is found to be $\tau_0 = 2.29\,\text{ms}$. Given that the coherence timescale of an atmospheric channel is on the order of milliseconds (see, e.g. [39]) and that a pulsed source with a $100\,\text{MHz}$ repetition rate (which corresponds to a pulse period of $0.01\,\mu\text{s}$) is readily available, the channel can be considered constant for a large number ($\sim\!10^{5}$) of sequential signal pulses.

Zoom In Zoom Out Reset image size

In figure 4, we plot the Probability Density Function (PDF) of the channel loss, achieved with different r_a values under different θ_z values, over a satellite-to-Earth channel. From this figure, we can observe that both the average channel loss (indicated by the centre of a PDF curve) and the fluctuation of channel loss (indicated by the spread of a PDF curve) decrease as the receiver aperture radius r_a increases. The former observation is intuitive since a larger receiver aperture can capture a larger portion of the beam on the ground. The latter observation is due to the well-known effect of aperture averaging in free-space optical communications—the level of power fluctuation decreases as the receiver aperture size increases (more details can be found in many good textbooks, e.g. [39]). It can be also seen that under each setting, the distribution of channel loss is approximately centred at diffraction loss, and the spread of distribution (which indicates the fluctuation of channel loss) is generally minimal. Such observations indicate that diffraction loss (a deterministic loss) is the dominant source of loss over a satellite-to-Earth channel. Indeed, it is well known that beam wandering, which is the major source of turbulence-induced channel loss fluctuations, is negligible in a satellite-to-Earth channel (see, e.g. [39, 50, 55, 56]). This phenomenon is because an optical beam transmitted from a satellite becomes very large when it enters the Earth's atmosphere, and the turbulence-induced effects are imposed on the optical beam only in the last segment of its propagation (note that the effective thickness of the atmosphere is only around $6\,\text{km}$—see, e.g. [57]).

Zoom In Zoom Out Reset image size

In figure 5, we plot the mean channel loss, achieved with different r_a values under different θ_z values, over a satellite-to-Earth channel. The results in this figure are consistent with those in figure 3 and provide a quantitative representation of the mean channel loss and the channel loss fluctuation over a satellite-to-Earth channel.

Zoom In Zoom Out Reset image size