Based on quantum mechanics, quantum key distribution (QKD) can realize security key sharing between remote communication parties, which plays an important role in our modern information society^[1–5]. In order to build the integrated quantum communication network, substantial progress has been achieved in both theoretical and experimental aspects, especially those that are fiber-based^[6–9] and satellite-based^[10–16] QKD. In 2021, the world’s first large-scale quantum communication network was constructed that integrated more than 700 terrestrial optical fiber QKD links and two high-speed QKD links in satellite-terrestrial free space, which enabled any user in the network to communicate in a distance of 4600 km^[16]. It shows that the quantum satellites can effectively provide QKD service even at the intercontinental scale and construct ultralong-distance global quantum networks. However, there are still some places where terrestrial fiber and ground stations cannot be constructed, like harsh mountainous areas, and air space above the sea. So airborne quantum nodes are expected to replace ground stations and provide flexible and relay links for large-scale integrated communication networks^[17]. Compared with satellite-to-station QKD, satellite-to-aircraft QKD features low atmospheric loss and long transmission distance, owing to the photon loss and turbulence predominantly occurring in the lower $\sim 10\mathrm{km}$ of the atmosphere^[18].

Since the first aircraft-based QKD experiment was verified successfully in 2013^[19], in the past decade, numerous studies have been focusing on the challenges of airborne QKD links^[19–25]. Compared with ground stations, airborne QKDs feature high-speed maneuverability and suffer complicated atmosphere conditions that include atmospheric turbulence^[26–30], background noise^[31–33], and attitude disturbance^[14]. Furthermore, a very thin layer of air will stick over the surface of the aircraft with high velocity, resulting in the boundary layer (BL)^[34,35]. It has been proved that the BL effect will seriously affect the airborne QKD performance when the aircraft speed is higher than 0.3 Ma^[36]. However, previous research of airborne QKD with BL effects is mainly under the air-to-ground scenario, while the satellite-to-aircraft scenario, which is one of the essential components in constructing a globe-wide quantum-secure communication network, has rarely been reported.

In this article, we propose a performance evaluation scheme of satellite-to-aircraft QKD with BL effects. We first propose a satellite-to-aircraft QKD scenario with decoy BB84 protocol. Then, the wavefront aberration of quantum signal states is evaluated by estimating the reflection index distribution of the surrounded BL and performing the ray tracing method by the Adams linear multistep method, which starts with the satellite ephemeris and aircraft trajectory. Afterward, the photon transmission efficiency caused by wavefront aberration is evaluated by the Strehl ratio. Finally, the overall photon quantum bit error rate (QBER) and final secure key rate can be estimated. The analyzed photon transmission loss in different incident angles shows that the effects of the BL are more serious when the aircraft moves towards the satellite. With common experimental settings, the BL would introduce a 31 dB loss to the transmitted photons, decrease $\sim 47\%$ of the quantum communication time, and decrease $\sim 51\%$ of the secure key rate. Our detailed satellite-to-aircraft QKD performance evaluation study can be performed on future airborne quantum communication designs.

The coordinates of satellite and aircraft in the WGS-84 coordinate system can be obtained from the satellite ephemeris and aircraft trajectory, as shown in Fig. 1. The points P and ${\mathit{P}}^{\prime }$ are the conventional terrestrial pole, the blue curve is the Greenwich meridian, $r$ is the mean radius of the Earth, $l_A$ is the longitude of the satellite in degrees, $L_A$ is the latitude of the satellite in degrees, $h_A$ is the satellite altitude, $l_B$ is the longitude of aircraft in degrees, $L_B$ is the latitude of aircraft in degrees, $h_B$ is aircraft altitude, and $\gamma$ is the angle between OA and OB.

In the ray-tracing methods, the satellite azimuth angle, the satellite elevation angle, and the distance between the satellite and aircraft are taken as the input parameters. Therefore, it is necessary to transform the WGS-84 coordinate system into the spherical coordinate system based on the aircraft, as shown in Fig. 2. The spherical coordinate system is established with the origin at the geometric center of the airborne receiver telescope. The $x$ axis is parallel to the tangent direction of the longitude and points north. The $y$ axis is parallel to the tangent direction of the longitude and points east. The $z$ axis merges with the local vertical and points to the zenith. The positions of satellite and aircraft in the coordinate system are S and O. The azimuth angle $\alpha$, the elevation angle $\beta$, and the distance $d$ between satellite and aircraft can be calculated by the following formula: $$\alpha =\arctan \left[\frac{\sin (|l_A-l_B|)\cos (L_A)}{\cos (L_B)\sin (L_A)-\cos (L_A)\sin (L_B)\cos (L_B-L_A)}\right],$$(1)$$d=(r+h_A){\left[1+{\left(\frac{r+h_B}{r+h_A}\right)}^2-2\left(\frac{r+h_B}{r+h_A}\right)\cos (\gamma )\right]}^{1/2},$$(2)$$\cos (\gamma )=\cos (L_B)\cos (L_A)\cos (l_A-l_B)+\sin (L_B)\sin (L_A),$$(3)$$\beta =\arccos \left[\frac{(r+h_A)\sin (\gamma )}{d}\right].$$(4)

The aero-optical effects are fundamentally caused by the gradient refractive index $n$ due to the variable-density flow field, which is expressed by the Gladstone–Dale equation^[37], $$n=1+\rho K_{\mathrm{GD}},$$(5)where $\rho$ is the density of the flow field. $K_{\mathrm{GD}}$ is the Gladstone–Dale constant decided only by the wavelength $\lambda$ (µm) of photons^[37], $$K_{\mathrm{GD}}=2.23\times {10}^{-4}\times (1+\frac{7.52\times {10}^3}{{\lambda }^2}).$$(6)

The trajectory of a ray in inhomogeneous media is determined by solving the ray equation^[38], $$\frac{\mathrm{d}}{\mathrm{d}p}\left(n\frac{\mathrm{d}s}{\mathrm{d}p}\right)=\nabla n,$$(7)where $s$ is the position vector of a typical point on the ray, $p$ is the path length of the ray, and $\nabla n$ is the gradient of the refractive index. Equation (7) can be written as a set of first-order differential equations, $$\{\begin{array}{l}\frac{\mathrm{d}F}{\mathrm{d}p}=\nabla n\hfill \\ \frac{\mathrm{d}s}{\mathrm{d}p}=\frac{1}{n}F\hfill \end{array},$$(8)where $F$ is the vector of light. The vectors $s$ and $F$ can be written as $$s=\left(\begin{array}{l}x\hfill \\ y\hfill \\ z\hfill \end{array}\right),F=\left(\begin{array}{l}{F}_x\hfill \\ F_y\hfill \\ F_z\hfill \end{array}\right)=n\left(\begin{array}{l}\frac{\mathrm{d}x}{\mathrm{d}p}\hfill \\ \frac{\mathrm{d}y}{\mathrm{d}p}\hfill \\ \frac{\mathrm{d}z}{\mathrm{d}p}\hfill \end{array}\right).$$(9)

The ray path can be calculated by numerical solution^[39,40].

In the airborne QKD procedure, the BL effect is also called the aero-optical effect in classical optics. Aero-optical effects will be introduced to the photons, which are propagated through the density-varying flow field of the BL. Typical aero-optical effects mainly include wavefront aberration, jitter, intensity attenuation, and so on. Relevant parameters of aero-optical effects are the optical path length (OPL), the optical path difference (OPD), and the Strehl ratio (SR)^[41].

The refractive index field of the airborne BL can be calculated by dividing the density field $\rho$ into sufficiently small squares and performing the Gladstone–Dale equation. The scattered photon path $P$ through the BL can be calculated by performing the ray-tracing methods.

The OPL of the photons is calculated by integrating the refractive index $n$ along the propagation path $P$^[42,43], $$\mathrm{OPL}(x,y,t)={\int }_{P}n(x,y,t)\mathrm{d}p.$$(10)

OPD shows the configuration of the wavefront and is defined as $$\mathrm{OPD}(x,y,t)=\mathrm{OPL}(x,y,t)-\overline{\mathrm{OPL}}.$$(11)

The overline denotes the spatial average over the optical aperture. The phase difference of the photons can be defined by $$\varphi =\frac{2\pi ·\mathrm{OPD}}{\lambda }.$$(12)

There is a distance between the BL and receiver telescope, as shown in Fig. 3, which depends on whether the aircraft is a transmitter or receiver in the airborne QKD scenario. When the transmission distance is similar to the communication distance, all aero-optical effects introduced by the BL need to be considered. However, when the distance is far less than communication distance, even if the effect of the BL is introduced into the divergence angle and the deflection angle, the deflection and divergent effect could be ignored, but the wavefront aberration should be taken into consideration. As the distance is far less than communication distance in the satellite-to-aircraft downlink QKD scenario, only the wavefront aberration would have been taken into consideration.

Previous research on airborne QKD with BL effects is mainly under the air-to-ground scenario, with the aircraft as the transmitter and the ground station as the receiver. In this paper, the satellite is the transmitter and the aircraft is the receiver. The primary difference is that the distance between the BL and the receiving telescope is long in the air-to-ground scenario, whereas the distance is tight in the satellite-to-aircraft scenario. As mentioned in Section 2.3, due to the different distance, the performance evaluation schemes are different in the two scenarios. The performance evaluation procedure for satellite-to-aircraft QKD scheme is shown in Fig. 4.

The satellite-to-aircraft downlink QKD scenario is shown in Fig. 2, the quantum photon source is located at the satellite (Alice), and the QKD receiving module with a spatial single-mode receiver is fixed in the upper fuselage of the aircraft (Bob). The receiver telescope position is temporarily set in a reasonable range. Assume that Alice is flying with a constant velocity, direction, and altitude. Given the aircraft specification, speed $v$, flying altitude $h_B$, and air density ${\rho }_h$, the density field distribution of the BL can be simulated by computational fluid dynamics software (such as CFX, Fluent, Star-CD, and COMSOL).

According to Section 2.1, the time-varying data of azimuth angle $\alpha$ ($0^{°}\le \alpha {<360}^{°}$), the elevation angle $\beta$ ($0^{°}\le \beta \le {90}^{°}$), and the distance $d$ between satellite and aircraft can be calculated by the satellite ephemeris and aircraft trajectory. Generally, the interval of $0^{°}\le \beta {<10}^{°}$ is used for links calibration.

The performance evaluation procedure for satellite-to-aircraft QKD scheme mainly contains three steps: photon aberrations evaluation, transmission efficiency calculation, and key rate estimation, as shown in Fig. 4.

According to the downlink satellite-to-aircraft QKD scenario, the distance between the BL and the receiving telescope is tiny. The deflection and divergent effect from aero-optics can be ignored, and the effective beam waist of the downlink photon at the receiving telescope is constant, no matter whether the BL exists or not.

When the Gaussian mode beam is propagating through the BL to the aircraft, the effective beam waist ${\omega }_{\mathrm{DP}}$ of the downlink photon at the receiving telescope can be expressed as $${\omega }_{\mathrm{DP}}=\sqrt{{\omega }_D^2+{({\sigma }_T·d)}^2},$$(13)where ${\sigma }_T$ is the pointing error of the transmitter telescope.

${\omega }_D$ is the beam waist at the ground station prior to pointing errors, $${\omega }_D=d\frac{\lambda }{\pi ·{\omega }_0}{\left[1+\frac{0.83}{\sin \beta }{\left(\frac{D_T}{r_0}\right)}^{5/3}\right]}^{3/5},$$(14)where $r_0$ is the Fried parameter in zenith^[44], and $D_T$ is the diameter of the transmitter telescope. ${\omega }_0=0.316D_T$ is the waist radius of the transmitted Gaussian beam^[44].

According to Section 2.1, we establish the spherical coordinate system with the origin at the geometric center of the airborne receiver telescope, as shown in Fig. 5. We take sufficient incident points uniformly in the light spot range with the radius of the effective beam waist ${\omega }_{\mathrm{DP}}$, and the initial value of $s$ in Section 2.2 can be denoted as $$s_0=\left[\begin{array}{ccc}{x}_1& y_1& z_1\\ x_2& y_2& z_2\\ ⋮& ⋮& ⋮\\ x_m& y_m& z_m\end{array}\right],$$(15)where $m$ is the number of incident points. The initial value of $F$ can be denoted as $$F_0=n\left(\begin{array}{c}\cos \alpha \cos \beta \\ \sin \alpha \cos \beta \\ \sin \beta \end{array}\right),$$(16)where $n$ is the refractive index field of the BL, which can be calculated by the density field distribution in Eq. (5). And, using fourth-order Runge–Kutta method, Eq. (8) can be expanded as $$\{\begin{array}{l}\begin{array}{l}{r}_{n+1}=r_n+\frac{1}{6}h\{T_n+2[T_n+\frac{h}{2}D(r_n)]+2[T_n+\frac{h}{2}D(r_n+\frac{h}{2}{T}_n)]\}\hfill \\ +\frac{1}{6}h\{T_n+hD[r_n+\frac{h}{2}{T}_n+\frac{h^2}{4}D(r_n\left)\right]\}\hfill \end{array}\hfill \\ \begin{array}{l}{T}_{n+1}=T_n+\frac{1}{6}h\{2D(r_n+\frac{h}{2}{T}_n)+2D[r_n+\frac{h}{2}{T}_n+\frac{h^2}{4}D(r_n)]\}\hfill \\ +\frac{1}{6}h\{D(r_n)+D[r_n+h{T}_n+\frac{h^2}{2}D(r_n)]\}\hfill \end{array}\hfill \end{array},$$(17)where $h$ is the step size of the ray-tracing method along the negative $z$ axis, and $D$ is the partial differential of $n$, $$D=n\left(\begin{array}{l}\frac{\partial n}{\partial x}\hfill \\ \frac{\partial n}{\partial y}\hfill \\ \frac{\partial n}{\partial z}\hfill \end{array}\right).$$(18)

The initial values $s_0$ and $F_0$ are substituted into Eq. (17) to calculate the next position $s_1$, and vector $F_1$. And then, the entire ray-tracing process can be calculated. But, to speed up the computer calculation, the values of $s_0$, $F_0$, $s_1$, $F_1$, and the refractive index field $n$ are taken as the input parameters for the Adams linear multistep method that is combined with the interpolation method based on polynomial fitting.^[45] Thus, when the Gaussian mode beam is propagating through the BL to the aircraft, the OPL of the photons can be evaluated with the ray-tracing method. Finally, the wavefront phase difference of photons can be calculated by the equation in Section 2.3.

When the beam propagates through the BL and illuminates the receiving telescope, the transmission efficiency ${\eta }_0$ can be calculated as^[44]$${\eta }_0=\mathrm{SR}·\exp (-\frac{\tau }{\sin \beta })·\{1-\exp [-0.5{\left(\frac{D_R}{{\omega }_{\mathrm{DP}}}\right)}^2]\},$$(19)where $D_R$ is the diameter of the receiving telescope, and $\tau$ is the extinction optical thickness from altitude to infinity,^[46] which is the reduction coefficient in brightness of stellar objects as their photons pass through the atmosphere. The effects of extinction depend on transparency, the elevation of the observer, the zenith angle, and the angle from the zenith to one’s line of sight. SR is the Strehl ratio^[47], $$\mathrm{SR}\approx \exp (-{\varphi }_{\mathrm{rms}}^2).$$(20)

In the satellite-to-aircraft QKD system, the photon transmission efficiency $\eta$ will be decreased, with the aero-optical effects of the aircraft BL, which can be calculated as $$\eta ={\eta }_0{\eta }_s{\eta }_d,$$(21)where ${\eta }_s$ is the system receiving efficiency caused by constant optical components and ${\eta }_d$ is the detector efficiency.

The decoy state is a common method in implemented experiments that combine with the QKD protocols like BB84 and measurement device independent (MDI) QKD, which can efficiently defend against the photon number splitting attacks and can perform the weak coherent photon source to replace the single-photon source in the implementations. The decoy state QKD protocol has been widely performed in fiber-based, satellite-based, and airborne-based QKD systems. Thus, in the satellite-to-aircraft QKD scheme, we perform weak-vacuum decoy BB84^[48] protocol with signal photon intensity $\mu$ and decoy photon intensity $\nu$. The modulating probability of signal (decoy) states is $P_s$ and $P_d$. The final secure key rate can be calculated as $$R\ge q\{Q_1[1-H_2(e_1)]-Q_{\mu }f(E_{\mu })H_2(E_{\mu })\},$$(22)where $Q_1$ is the gain of the received single-photon states, $e_1$ is the error rate of single-photon states, and $f(x)$ is the information reconciliation efficiency for correcting error bits. $Q_{\mu }$ and $E_{\mu }$ represent the gain of signal states and the overall QBER, respectively. $H_2(x)$ is the binary Shannon entropy, which can be calculated as $$H_2(x)=-x\log (x)-(1-x)\log (1-x).$$(23)

Given the photon transmission efficiency $\eta$, $Q_{\mu }$ is calculated as $$Q_{\mu }=Y_0+1-e^{-\eta \mu },$$(24)where $Y_0$ is the dark count rate of QKD systems. Thus, the error gain of signal quantum states can be given by $$E_{\mu }{Q}_{\mu }=e_0{Y}_0+e_d(1-e^{-\eta \mu }),$$(25)where $e_0$ is the error rate of dark counts, usually $e_0=0.50$. $e_d$ is the misalignment error rate of QKD systems.

Thus, the QBER $E_{\mu }$ can be calculated as $$E_{\mu }=E_{\mu }{Q}_{\mu }/Q_{\mu }.$$(26)

The gain of single-photon states $Q_1$ can be calculated as $$Q_1\ge Q_1^{L,\nu ,0}=\frac{{\mu }^2{e}^{-\mu }}{\mu \nu -{\nu }^2}(Q_{\nu }{e}^{\nu }-Q_{\mu }{e}^{\mu }\frac{{\nu }^2}{{\mu }^2}-\frac{{\mu }^2-{\nu }^2}{{\mu }^2}{Y}_0),$$(27)where $L$ denotes the lower bound value, and $Q_{\nu }$ is the gain of decoy states.

The error rate of single-photon states $e_1$ can be calculated as $$e_1\le e_1^{U,\nu ,0}=\frac{E_{\nu }{Q}_{\nu }{e}^{\nu }-e_0{Y}_0}{Y_1^{L,\nu ,0}\nu },$$(28)where $Y_1^{L,\nu ,0}$ is the yield of single-photon states, $$Y_1^{L,\nu ,0}=\frac{Q_1}{\mu e^{-\mu }}.$$(29)

The error gain of decoy states $E_{\nu }{Q}_{\nu }$ can be calculated as $$E_{\nu }{Q}_{\nu }=e_0{Y}_0+e_d(1-e^{-\eta \nu }).$$(30)

Afterward, the QBER $E_{\mu }$ and secure key rate $R$ can be obtained by taking the photon transmission efficiency $\eta$ into weak-vacuum decoy BB84 protocol.

The specific parameters of the aircraft, quantum photon source payload, and optical receiving module are shown in Table 1. Some parameters refer to the Micius quantum experiment science satellite in the literature^[10].

Here, the typical wing-body configuration DLR-F6 is chosen for the performance analysis of our specified satellite-to-aircraft QKD system. Given the detailed aircraft description from Table 1, the BL will be generated around the DLR-F6^[49] and its density field distribution can be simulated by the computational fluid dynamics software (Ansys Fluent), as shown in Fig. 6. Assuming that the receiver telescope that conforms with the aircraft is fixed in the symmetry axis of the upper fuselage surface, the $500\mathrm{mm}\times 500\mathrm{mm}\times 50\mathrm{mm}$ BL above the fuselage has been taken, as shown in Fig. 6. Afterward, the refractive index distribution can be calculated by Eq. (5). For convenience, the projecting distances between the geometric center of the receiver telescope and the tip of the nose in the flight direction are denoted as $X_0$. The value of $X_0$ should be greater than 37 mm to avoid the flight deck, which is temporarily set as 66 mm.

In the coordinate system shown in Fig. 6, assume that the satellite orbits the earth at an altitude of 500 km that is directly over the aircraft at various azimuths. Setting azimuth angle α = 0°, 90°, 180°, 270° and elevation angle ${10}^{°}\le \beta \le {90}^{°}$, the total loss at different incident angles can be calculated as shown in Fig. 7. As the elevation angle increases, the total loss declines sharply, which indicates that the loss is generally high when the elevation angle $\beta \le {30}^{°}$. Moreover, when $\alpha {=180}^{°}$ and ${30}^{°}\le \beta \le {90}^{°}$, the total loss also increases on average about 3 dB. This indicates that when the aircraft moves toward the satellite, the effects of the BL are more serious. So, to reduce the influence of the BL, with permission to adjust the direction of flight, it is possible to change the direction of flight and turn the aircraft away from the satellite.

Assume that the aircraft is flying at a constant speed $v=0.7\mathrm{Ma}$, altitude $h_B=10\mathrm{km}$, and heads due south or east, and the initial position of aircraft is (${34}^{°}{15}^{\prime }{56}^{\prime \prime }\mathrm{N}$, ${108}^{°}{57}^{\prime }{13}^{\prime \prime }\mathrm{E}$). It is assumed that the aircraft returns to the initial position after each orbit of communication is completed. We import the satellite ephemeris of Micius and the aircraft trajectory into the Satellite Tool Kit (STK) from 12:00 on May 29, 2022, to 12:00 on June 5, 2022. The schematic diagram of satellite-to-aircraft QKD is shown in Fig. 8; all the trajectories of Micius that establish available communication links of satellite-to-aircraft QKD are marked.

Then, during a week, the performance of the whole satellite-to-aircraft QKD session is evaluated. When the aircraft is heading south, the result is shown in Fig. 9. Significantly, the abscissa represents the link time of all established links in a week, and the total link time is 6800 s. The additional channel loss to the transmitted photons is around 30 dB during the total link time, as shown in Fig. 9(a). The time when the secure key rate is more than zero is denoted as the quantum communication time. Therefore, the total quantum communication time is 3625 s, as shown in Fig. 9(c), and the estimated total final key size is around $6.718\times {10}^6\mathrm{bits}$. When the aircraft is heading east, the result is shown in Fig. 10. The additional channel loss to the transmitted photons is around 32 dB during the total link time, as shown in Fig. 10(a). The total quantum communication time is 3685 s, as shown in Fig. 10(c), and the estimated total final key size is around $6.322\times {10}^6\mathrm{bits}$. If there is no BL surrounding the aircraft, the estimated secure key rate would be around $1.326\times {10}^7\mathrm{bits}$. In summary, the BL effects cannot be ignored in the satellite-to-aircraft QKD scenario and heavily decrease the final secure key rate.

Airborne QKD will be a flexible bond between the terrestrial fiber QKD network and quantum satellites, which can establish a mobile, on-demand, and real-time coverage quantum network. However, the randomly distributed BL always surrounds the surface of the aircraft, which would introduce random wavefront aberration, jitter, and extra intensity attenuation to the transmitted photons between the aircraft and the ground station. Previous research of airborne QKD with BL effects is mainly under the air-to-ground scenario, while the satellite-to-aircraft scenario is rarely reported. In this article, we proposed a detailed performance evaluation scheme of satellite-to-aircraft QKD with BL effects. The analyzed photon transmission loss shows that the effects of the BL are more serious when the aircraft moves towards the satellite. In our proposed satellite-to-aircraft QKD scenario, the BL would introduce $\sim 31\mathrm{dB}$ loss to the transmitted photons, decrease $\sim 47\%$ of the quantum communication time, and decrease $\sim 51\%$ of the secure key rate. This indicates that the aero-optical effects caused by the BL cannot be ignored. Our detailed satellite-to-aircraft QKD performance evaluation study can be used in future airborne quantum communication designs.