Airborne Quantum Key Distribution Performance Analysis under Supersonic Boundary Layer

Airborne quantum key distribution (QKD) that can synergize with terrestrial networks and quantum satellite nodes is expected to provide flexible and relay links for the large-scale integrated communication network. However, the photon transmission rate would be randomly reduced, owing to the random distributed boundary layer that surrounding to the surface of the aircraft when the flight speed larger than Mach 0.3. Here, we investigate the airborne QKD performance with the BL effects. Furthermore, we take experimental data of supersonic BL into the model and compare the airborne QKD performance under different conditions. Simulation results show that, owing to the complex small-scale turbulence structures in the supersonic boundary layer, the deflection angle and correspondingly drifted offset of the beam varied obviously and randomly, and the distribution probability of photons are redistributed. And the subsonic and supersonic boundary layer would decrease ~35.8% and ~62.5% of the secure key rate respectively. Our work provides a theoretical guidance towards a possible realization of high-speed airborne QKD.


Introduction
Quantum key distribution (QKD) is a promising technology to realize secure transmission, the fundamental principles of quantum mechanics guarantee the safety of QKD, which ensures the generation of information-theoretical-secure keys for distant users [1][2][3][4][5]. Furthermore, QKD can resist eavesdropping detection and tamper relying on physical principles, which has tremendous utilization potential in terms of finance, government, and military [6,7]. Up to now, remarkable progress has been achieved both in fiber links [8][9][10][11] and quantum satellites [12][13][14][15][16][17][18][19] QKD systems, and gradually transferred from the laboratory to realistic applications. However, with constant orbits, limited communication time window and night-only quantum satellites, airborne [20][21][22][23][24][25][26][27] platforms are ideal mobile nodes that can synergize with terrestrial links and quantum satellites to build the mobile, on-demand, and real-time coverage quantum network. For example, airborne quantum nodes could serve as temporary relays to solve the last-mile quantum key exchange for an inner-city or a field network benefiting from their rapid deployment capabilities.
Several groups have reported their research in airborne quantum communications in recent years, the first demonstration was implemented in 2013, with the platform flying at the speed of 290 km/h and a height of 1.1 km [20].
In 2020, Nanjing University accomplished the first drone-based entanglement distribution [23] and achieved 200 m coverage and duration of 40 min, and in the next year, In 2020, Nanjing University accomplished the first drone-based entanglement distribution [23] and achieved 200 m coverage and duration of 40 min, and in the next year, they expand the distance to 1000 m by means of one more drone [21] as an optical-relay. Compared to satellite quantum communication, airborne QKD features in high-speed maneuverability and suffers from the complicate atmosphere environment, including atmospheric turbulence, background noise and attitude disturbance, worse still, the boundary layer (BL) would also severely impact the quantum signal transmission. When the aircraft is flying at a high speed, usually larger than 0. 3 Ma, a very thin layer of air will stick over the surface of the aircraft with high velocity, resulting in the BL, and it would introduce random disturbance to the transmitted photons, which would reduce the coupling efficiency and fidelity of quantum states [28], and decrease the performance of aircraft-based QKD. Currently, most implementations of airborne QKD only analyzed the noise from atmospheric turbulence and molecular scattering [25,26], but ignored the boundary layer effects. In 2021, our group proposed an airborne QKD performance evaluation scheme that takes the boundary layer effect into account. The analysis results show that the aerodynamic optical effects caused by the boundary layer should be attended to, which will greatly reduce the final key rate [29].
In this article, an airborne QKD performance evaluation scheme with supersonic BL effect is proposed. We first simulate the density and refractive index field of the BL. Then, the photon propagation model in the boundary layer is established by performing the ray tracing method that uses the Adams linear multistep method. In the end, we present the analysis results about the key rate of photons transmitted through the measured supersonic BL. Owing to the complex small-scale turbulence structures in the supersonic boundary layer, the deflection angle and correspondingly drifted offset of the beam varied obviously and randomly, and the distribution probability of photons is redistributed. The result shows that the subsonic and supersonic boundary layer would decrease ~35.8% and ~62.5% of the secure key rate respectively. With the increase of speed v, the key-rate curve is obviously jittering, and the QKD performance is continuously reduced. Our work provides a theoretical guidance towards a possible realization of high-speed airborne QKD.

Background
In the air-to-ground communication scenario of airborne QKD, we assume is the shortest projection distance to the horizon plane, between the aircraft (Alice) and the receiving ground station (Bob). Alice flies at a constant speed in a certain path obliquely above Bob, as shown in Figure 1. ℎ is the relative height of Alice and Bob which is constant. is the distance between Alice and Bob. is zenith angle. During the entire flight, the azimuth angle of the aircraft is varying with a range of [−80°, 80°].

Photon Scattering under the Boundary Layer
Given the density field distribution of the boundary layer, the refractive index distribution of the flow field can be calculated by where ρ is the density, and n is the refractive index of BL. K GD is the Gladstone-Dale constent with unit of m 3 /kg. K GD is usually calculated by the empirical formula where λ is the wavelengthof photons with unit of micron. With the settings in Table 1, the G-D coefficient is 2.237 × 10 −4 m 3 /kg. The scattered photon path through the boundary layer can be calculated by ray tracing methods [30]. With the varying refractive index, the boundary layer can be divided into several sufficiently small cube units. And the internal refractive index of each unit (i, j, k) is uniform which is n(i, j, k). Then, the photon will be refracted at the boundary of the two grids, as shown in Figure 2, and the new direction of communication can be calculated by Snell law where n(x 1 , y 1 , z 1 ) is the refractive index of unit (x 1 , y 1 , z 1 ), θ 1 is the incident angle, θ 2 is the refraction angle. In the air-to-ground communication scenario, the incident angle θ 1 is equal to zenith angle ϕ which can be calculated as When the azimuths angle is 0 • , the incident angle θ 1 = ϕ = arctan(d/h). It is worth noting that when the azimuth angle is 0 • , the incident angle is not 0 • .  And the new incident angle can be calculated by 2 when the photon propagated to the next grid. Then, recording the whole location ( 1 , 1 , 1 ), ( 2 , 2 , 2 ) … ( , , ) by reciprocating the Snell law, the scattered photon path X via the boundary layer can be calculated.
Then, the wavefront information could be calculated. Optical path length (OPL) is calculated by integrating the refractive index n along the propagation path X [30].
And the optical path difference (OPD) is defined as follows, Afterward, assumed that the photon source conforms to the distribution of the Gaussian beam, the normalized intensity of the Gaussian beam can be expressed as where r is the radial distance from the center axis of the beam, l is the axial distance from the beam's focus, and is the effective beam waist of the downlink signal at the ground station [31] Here, we assumed that the Transmitter pointing precision = 150 μrad that is given in [21].
is the beam waist at the Ground station prior to pointing errors: where r0 = 0.4 m is the Fried parameter in zenith, and is zenith angle as shown in Figure  1. ω0 is the waist radius of the Gaussian beam The angled brackets denote the spatial average over the optical aperture.

Transmission Efficiency Analysis
Assuming that the ATP (Acquisition, Tracking, and Pointing) technique is perfect, the center axis of the beam can be aligned with the center of the receiver telescope. And then, the transmission efficiency from aircraft to the ground station can be calculated when the light is via the boundary layer or not. And the new incident angle can be calculated by θ 2 when the photon propagated to the next grid. Then, recording the whole location (x 1 , y 1 , z 1 ), (x 2 , y 2 , z 2 ) . . . (x n , y n , z n ) by reciprocating the Snell law, the scattered photon path X via the boundary layer can be calculated.
Then, the wavefront information could be calculated. Optical path length (OPL) is calculated by integrating the refractive index n along the propagation path X [30].
And the optical path difference (OPD) is defined as follows, Afterward, assumed that the photon source conforms to the distribution of the Gaussian beam, the normalized intensity of the Gaussian beam can be expressed as where r is the radial distance from the center axis of the beam, l is the axial distance from the beam's focus, and ω LP is the effective beam waist of the downlink signal at the ground station [31] Here, we assumed that the Transmitter pointing precision σ T = 150 µrad that is given in [21]. ω L is the beam waist at the Ground station prior to pointing errors: where r 0 = 0.4 m is the Fried parameter in zenith, and ϕ is zenith angle as shown in Figure 1. ω 0 is the waist radius of the Gaussian beam The angled brackets denote the spatial average over the optical aperture.

Transmission Efficiency Analysis
Assuming that the ATP (Acquisition, Tracking, and Pointing) technique is perfect, the center axis of the beam can be aligned with the center of the receiver telescope. And then, the transmission efficiency from aircraft to the ground station can be calculated when the light is via the boundary layer or not.
Without the boundary layer effects, the detected events fulfil the Gaussian distribution. When the beam illuminates the receiving telescope, the transmission efficiency η 0 can be calculated by geometric optics as [31] where D R is the diameter of the receiving telescope, and β is the extinction optical thickness between sea level and altitude. After passing through the boundary layer, the intensity of the beam will be redistributed. Then, by weighting each beam with the Gaussian distribution in Equation (7), the photons distribution probability eventually reaches the receiver telescope will be calculated, when the light is via the boundary layer or not. Then, the photons distribution probability ratio E boundary can be calculated as where ε boundary is the statistics of photons distribution probability within the range of receiving telescope per unit time with the effect of BL. While ε 0 is the value without the effect of BL. After passing through the boundary layer, the transmission efficiency η boundary can be calculated as where SR is the Strehl ratio [32] SR ≈ exp − 2πOPD rms where OPD rms is the RMS of the OPD on the optical aperture.

Secure Key Rate Estimation
The decoy state method [33], as an important weapon to combat photon number splitting attack, is proposed to use a weakly coherent light source in the QKD protocol to replace the ideal single-photon source that cannot be achieved at present. Therefore, the Vacuum + weak BB84 protocol is selected in the airborne QKD system, and a formula for secure key rate is where Q 1 is the gain of single-photon states, e 1 is the error rate of single-photon states, f (x) is the bidirectional error correction efficiency as a function of error rate, H 2 (x) is the binary Shannon information function and µ is the intensity of the signal state. Qµ and Eµ respectively represent the gain of signal states and the overall quantum bit error rate.
In free-space quantum communication, it is necessary to consider the reduction of efficiency caused by the diffusion of the light spot at the receiving terminal. With the influence of the boundary layer, the total transmission efficiency is where η s is receiving optical module efficiency, and η d is detector efficiency. Substituting Equations (16) into the formula about Vacuum + weak BB84 protocol in literature [33], and the secure key rate R can be calculated.

Evaluation Result and Discussion
Here, the typical airfoil that named "NACA0015" is chosen for the performance analysis of our specified air-to-ground QKD system. The specific parameter settings of the aircraft, source, ground station and protocols are shown in Table 1. The boundary layer will be generated around the airfoil and its density field distribution can be simulated by the computational fluid dynamics software (Ansys Fluent, Canonsburg, PA, USA) as shown in Figure 3a, with the parameters shown in Table 1. Usually, the boundary layer thickness based on the velocity boundary layer concept that the region in which flow adjusts from zero velocity at the wall to a maximum in the main stream of the flow is termed the boundary layer. According to the simulation results, the airfoil velocity boundary layer thickness is about 10 mm. Here, in order to ensure that the path lengths in the boundary layer are approximately equal and the results are more accurate under different incident angles, we expend the concept of the boundary layer as the field that caused by the motion of aircraft. So, the generalized boundary layer thickness is 400 mm.

Evaluation Result and Discussion
Here, the typical airfoil that named "NACA0015" is chosen for the performance analysis of our specified air-to-ground QKD system. The specific parameter settings of the aircraft, source, ground station and protocols are shown in Table 1. The boundary layer will be generated around the airfoil and its density field distribution can be simulated by the computational fluid dynamics software (Ansys Fluent, Canonsburg, PA, USA) as shown in Figure 3a, with the parameters shown in Table 1. Usually, the boundary layer thickness based on the velocity boundary layer concept that the region in which flow adjusts from zero velocity at the wall to a maximum in the main stream of the flow is termed the boundary layer. According to the simulation results, the airfoil velocity boundary layer thickness is about 10 mm. Here, in order to ensure that the path lengths in the boundary layer are approximately equal and the results are more accurate under different incident angles, we expend the concept of the boundary layer as the field that caused by the motion of aircraft. So, the generalized boundary layer thickness is 400 mm. Literature [34] used the nano-particle-based planar laser scattering technique to measure the density distribution of the supersonic (Ma = 3.0) turbulent boundary layers, as shown in Figure 3b. The velocity boundary layer thickness is about 10 mm [33], and the generalized boundary layer thickness is about 13 mm. Although this kind of boundary layer is not the actual aircraft boundary layer, the turbulence structure is similar to the actual situation. Due to the limited amount of data, we only analyze the results at a certain moment, without considering an average over enough runs. It is also meaningful to research the influence of this kind of boundary layer on QKD performance.
When the photon trajectories of different incident angles pass through the boundary layer and reached the ground station, the evaluated deflection angle and the drifted offset of the beam are shown in Figure 2. The deflection angle can be calculated as | 1 − | which is absolute value of the angle difference between the incident angle 1 and the last refracted angle when the photon propagated via the boundary layer. In both cases, the deflection angle and correspondingly drifted offset of the beam were varied obviously in Figures 4 and  5. Because of the complex small-scale turbulence structures in the supersonic boundary layer, the spatial distribution of the density field is anisotropic and random, and the results vary randomly as shown in Figure 5. When the azimuth angle is near zero, the incident photon will pass through the region with a large refractive index gradient, as shown in Figure 3b, which could influence the deflection angle obviously. In a practical airborne quantum communication system, owing to the ATP technique, the center axis of the beam can always be aligned with the center of the receiver telescope. Therefore, the deflection angle and the drifted offset were not analyzed emphatically in this paper, which could been ignored. Literature [34] used the nano-particle-based planar laser scattering technique to measure the density distribution of the supersonic (Ma = 3.0) turbulent boundary layers, as shown in Figure 3b. The velocity boundary layer thickness is about 10 mm [33], and the generalized boundary layer thickness is about 13 mm. Although this kind of boundary layer is not the actual aircraft boundary layer, the turbulence structure is similar to the actual situation. Due to the limited amount of data, we only analyze the results at a certain moment, without considering an average over enough runs. It is also meaningful to research the influence of this kind of boundary layer on QKD performance.
When the photon trajectories of different incident angles pass through the boundary layer and reached the ground station, the evaluated deflection angle and the drifted offset of the beam are shown in Figure 2. The deflection angle can be calculated as |θ 1 − θ n | which is absolute value of the angle difference between the incident angle θ 1 and the last refracted angle θ n when the photon propagated via the boundary layer. In both cases, the deflection angle and correspondingly drifted offset of the beam were varied obviously in Figures 4 and 5. Because of the complex small-scale turbulence structures in the supersonic boundary layer, the spatial distribution of the density field is anisotropic and random, and the results vary randomly as shown in Figure 5. When the azimuth angle is near zero, the incident photon will pass through the region with a large refractive index gradient, as shown in Figure 3b, which could influence the deflection angle obviously. In a practical airborne quantum communication system, owing to the ATP technique, the center axis of the beam can always be aligned with the center of the receiver telescope. Therefore, the  Afterward, the photons distribution probability that eventually reaches the receiver telescope at azimuth angle is 10°,was calculated as shown in Figure 6, when the photon via the supersonic BL or not. It shows that the photons distribution probability is redistributed when the scattered photon passes through the supersonic BL. But the photons distribution probability is fundamentally invariant when the photon via the BL when flight speed is 0.7 Ma, which is not displayed.   Afterward, the photons distribution probability that eventually reaches the receiver telescope at azimuth angle is 10°,was calculated as shown in Figure 6, when the photon via the supersonic BL or not. It shows that the photons distribution probability is redistributed when the scattered photon passes through the supersonic BL. But the photons distribution probability is fundamentally invariant when the photon via the BL when flight speed is 0.7 Ma, which is not displayed.  Afterward, the photons distribution probability that eventually reaches the receiver telescope at azimuth angle is 10 • ,was calculated as shown in Figure 6, when the photon via the supersonic BL or not. It shows that the photons distribution probability is redistributed when the scattered photon passes through the supersonic BL. But the photons distribution probability is fundamentally invariant when the photon via the BL when flight speed is 0.7 Ma, which is not displayed.  Afterward, the photons distribution probability that eventually reaches the receiver telescope at azimuth angle is 10°,was calculated as shown in Figure 6, when the photon via the supersonic BL or not. It shows that the photons distribution probability is redistributed when the scattered photon passes through the supersonic BL. But the photons distribution probability is fundamentally invariant when the photon via the BL when flight speed is 0.7 Ma, which is not displayed.  Then, the receiving photons distribution probability ratio E boundary in Equation (12) over the azimuth angle can be calculated, as shown in Figure 7, when the photon via the supersonic BL. The curve random jitter is more obvious, and E boundary > 1 at some azimuth angle which shows that photons bunching has happened. Moreover, when the refractive index gradient is less than zero, the refraction angle is larger than incident angle. In this case, the photon spot dispersion caused by the small incident angle is more serious due to the limited number of grids. As shown in Figure 4, The small incidence angle will cause the beam to separate to the two sides, intensify the diffusion of light spots, and reduce the receiving photons distribution probability ratio E boundary . The case of the small incident angle with the negative refractive index gradient is indicated by the black box in Figure 3b. So, it is possible that significant dispersion phenomenon occurs around the zero azimuthal angle. Then, the receiving photons distribution probability ratio in Equation (12) over the azimuth angle can be calculated, as shown in Figure 7, when the photon via the supersonic BL. The curve random jitter is more obvious, and > 1 at some azimuth angle which shows that photons bunching has happened. Moreover, when the refractive index gradient is less than zero, the refraction angle is larger than incident angle. In this case, the photon spot dispersion caused by the small incident angle is more serious due to the limited number of grids. As shown in Figure 4, The small incidence angle will cause the beam to separate to the two sides, intensify the diffusion of light spots, and reduce the receiving photons distribution probability ratio . The case of the small incident angle with the negative refractive index gradient is indicated by the black box in Figure 3b. So, it is possible that significant dispersion phenomenon occurs around the zero azimuthal angle. When Figure 6 is converted into a top view, as shown in Figure 8, the distribution of photons can be well displayed with and without the effect of supersonic BL. By comparing the two cases where the azimuth is 10° and 20°, the probability of photon appearing in the range of the receiving telescope after passing through the boundary layer is different. When the azimuth is 10°, the photons passe through the boundary layer and are partially refracted out of the scope of the telescope. When the azimuth is 20°, some photons that are not in the scope of the telescope are refracted into the telescope after passing through the boundary layer. By contrast, the case of > 1 can be explained. When Figure 6 is converted into a top view, as shown in Figure 8, the distribution of photons can be well displayed with and without the effect of supersonic BL. By comparing the two cases where the azimuth is 10 • and 20 • , the probability of photon appearing in the range of the receiving telescope after passing through the boundary layer is different. When the azimuth is 10 • , the photons passe through the boundary layer and are partially refracted out of the scope of the telescope. When the azimuth is 20 • , some photons that are not in the scope of the telescope are refracted into the telescope after passing through the boundary layer. By contrast, the case of E boundary > 1 can be explained. Then, the receiving photons distribution probability ratio in Equation (12) over the azimuth angle can be calculated, as shown in Figure 7, when the photon via the supersonic BL. The curve random jitter is more obvious, and > 1 at some azimuth angle which shows that photons bunching has happened. Moreover, when the refractive index gradient is less than zero, the refraction angle is larger than incident angle. In this case, the photon spot dispersion caused by the small incident angle is more serious due to the limited number of grids. As shown in Figure 4, The small incidence angle will cause the beam to separate to the two sides, intensify the diffusion of light spots, and reduce the receiving photons distribution probability ratio . The case of the small incident angle with the negative refractive index gradient is indicated by the black box in Figure 3b. So, it is possible that significant dispersion phenomenon occurs around the zero azimuthal angle. When Figure 6 is converted into a top view, as shown in Figure 8, the distribution of photons can be well displayed with and without the effect of supersonic BL. By comparing the two cases where the azimuth is 10° and 20°, the probability of photon appearing in the range of the receiving telescope after passing through the boundary layer is different. When the azimuth is 10°, the photons passe through the boundary layer and are partially refracted out of the scope of the telescope. When the azimuth is 20°, some photons that are not in the scope of the telescope are refracted into the telescope after passing through the boundary layer. By contrast, the case of > 1 can be explained.  And the photons distribution probability with several different azimuth angles can be calculated, when the photons via the supersonic BL or not. As shown in Figure 9, different from the subsonic boundary layer, the supersonic boundary layer can diffuse the probability of photons distribution. It can cause the photon of receiving terminal to be so irregular that the transmission efficiency and secure key rate decrease.
Entropy 2023, 25, x FOR PEER REVIEW 9 of 12 receiving photons distribution when the azimuth angle is 20° with and without the effect of supersonic BL. The green dotted line indicates the aperture of the receiving telescope.
And the photons distribution probability with several different azimuth angles can be calculated, when the photons via the supersonic BL or not. As shown in Figure 9, different from the subsonic boundary layer, the supersonic boundary layer can diffuse the probability of photons distribution. It can cause the photon of receiving terminal to be so irregular that the transmission efficiency and secure key rate decrease. By taking the total transmission efficiency into the formula of the key rate, the QKD performance with the BL effects is evaluated and the result is shown in Figure 10. When the photon via the BL, the key rate curve drops obviously. Compared with the subsonic surface layer, the effect of the supersonic surface layer is more obvious. The curve varies even more dramatically, and the effect of supersonic BL at individual points is very obvious. Therefore, the estimated average secure key rate is around 943 bit/s and 551 bit/s respectively when the photon via the subsonic and supersonic surface layer. If there's no boundary layer surrounding the aircraft, the estimated average secure key rate would be around 1468 bit/s. It shows that the actual supersonic BL can also have an impact on the QKD performance, and the effect is more obvious. In addition, due to the higher speed, the effective communication time between the vehicle and the ground station will be shorter, and the total key rates will be smaller. Further, although the impact of the deflection angle and the drifted offset could be ignored, some potential factors will also affect the performance of QKD which from ATP system. Therefore, it is feasible to implement QKD by supersonic aircraft, but the results would be unsatisfactory. By taking the total transmission efficiency into the formula of the key rate, the QKD performance with the BL effects is evaluated and the result is shown in Figure 10. When the photon via the BL, the key rate curve drops obviously. Compared with the subsonic surface layer, the effect of the supersonic surface layer is more obvious. The curve varies even more dramatically, and the effect of supersonic BL at individual points is very obvious. Therefore, the estimated average secure key rate is around 943 bit/s and 551 bit/s respectively when the photon via the subsonic and supersonic surface layer. If there's no boundary layer surrounding the aircraft, the estimated average secure key rate would be around 1468 bit/s. It shows that the actual supersonic BL can also have an impact on the QKD performance, and the effect is more obvious. In addition, due to the higher speed, the effective communication time between the vehicle and the ground station will be shorter, and the total key rates will be smaller. Further, although the impact of the deflection angle and the drifted offset could be ignored, some potential factors will also affect the performance of QKD which from ATP system. Therefore, it is feasible to implement QKD by supersonic aircraft, but the results would be unsatisfactory. Entropy 2023, 25, x FOR PEER REVIEW 10 of 12 Figure 10. The secure key rate over the azimuth angle when the photon via subsonic BL, supersonic BL, or not.

Conclusions
Airborne quantum key distribution (QKD) that can synergize with terrestrial networks and quantum satellite nodes is expected to provide flexible and relay links for the large-scale integrated communication network. However, the photon transmission rate would be randomly reduced, owing to the random distributed boundary layer that surrounding to the surface of the aircraft when the flight speed larger than Mach 0.3, which would change the local refractive index and energy flux density drastically. Different from satellite based implementations, airborne platforms need to consider the influence of boundary layer due to the atmosphere. In this article, an airborne QKD performance evaluation scheme with supersonic BL effect is proposed. Through modeling and analysis, owing to the complex small-scale turbulence structures in the supersonic boundary layer, the deflection angle and correspondingly drifted offset of the beam varied obviously and randomly, and the distribution probability of photons are redistributed. The result shows that the subsonic and supersonic boundary layer would decrease ~35.8% and ~62.5% of the secure key rate respectively. Our work provides a theoretical guidance towards a possible realization of high-speed airborne QKD.