Atmospheric channel model for QKD over satellite megaconstellations

Introduction

The growing interest in quantum communication networks has driven significant research efforts toward establishing reliable ground-to-satellite quantum links. Unlike classical optical communication, quantum key distribution (QKD) and other quantum protocols are highly sensitive to photon loss, noise and decoherence introduced by the transmission medium. The atmospheric channel, in particular, imposes a fundamental challenge, as turbulence, absorption and scattering translates into severe constraints on the fidelity and efficiency of quantum states transmission.

Ground-to-satellite quantum links have emerged as a promising approach to extend quantum networks over global distances, bypassing the exponential loss limitations of fiber-based infrastructures. Several milestone experiments [1], most notably those carried out with the Micius satellite, have demonstrated the feasibility of satellite-based QKD, entanglement distribution and quantum teleportation over distances exceeding a thousand km. However, the performance of such systems depends critically on the atmospheric conditions encountered along the optical path, which vary with altitude, geographic location, season and time of day.

Accurate modeling of the atmospheric channel is therefore essential for the design, optimization and performance prediction of quantum communication systems operating along ground-to-satellite links. Key physical phenomena that must be accounted for include beam wandering and broadening induced by refractive index fluctuations, transmittance fluctuations, and altitude-dependent absorption due to molecular and aerosol species. Concepts borrowed from thermodynamics and optics, such as density-dependent refractive indices and Rayleigh scattering will provide the mathematical description to tackle these physical phenomena.

In this work, we present a comprehensive atmospheric channel model tailored to ground-to-satellite quantum optical links. The model integrates turbulence-induced index variations, beam propagation effects, and wavelength-dependent transmittance profiles derived from standard atmospheric databases. We evaluate the impact of these factors on key classical and quantum communication metrics, such as the single-photon transmission probability and losses, power balance and channel capacity. Furthermore, we benchmark our model against data extracted from real-world quantum satellite missions such as Micius and Hispasat. This comparative analysis allows us to identify the strengths and limitations of this approach, prove the feasibility of space Quantum Key Distribution and to explore future quantum satellite missions.

The remainder of the paper is organized as follows. We briefly review prior works on the feasibility of satellite QKD in Section II. Our physical modeling of the atmospheric optical channel is presented in Section III, and the technical specifications relevant for QKD satellite missions are thoroughly discussed in Section IV. The validation of our model under realistic parameterization and a comparison with experimental values is the focus of § 5. Finally, Section VI summarizes the contributions of this work.

Related Work

The feasibility of satellite-based quantum communication was first theoretically analyzed in [2], who identified the low photon loss of free-space links at high altitudes as a key advantage over fiber-based networks. Subsequent important works carried out in [3] and [4] laid the groundwork for satellite QKD, addressing photon loss budgets and link efficiency under realistic atmospheric conditions.

On the experimental front, the launch of the Micius quantum science satellite in 2016 marked a turning point in the field. [5] demonstrated intercontinental QKD over a distance exceeding 7,600 km using a trusted-node architecture, while [6] reported the distribution of entangled photon pairs between ground stations separated by over 1,200 km. These experiments confirmed the practical viability of satellite-mediated quantum links and established key benchmarks for atmospheric channel loss. The atmospheric channel has been extensively studied in the physical domain. The Kolmogorov turbulence theory [7] and the Hufnagel-Valley refractive index structure parameter profile [8] remain the dominant analytical frameworks for characterizing beam propagation through turbulent media. More sophisticated models that study the small-scale fading and interference statistics of a signal include generalizations of other well-known distributions such as Gaussian [9] or Rayleigh distributions, called the $\lambda-\kappa-\mu$ distributions [10], [11], [12], as well as other generalized distributions like $\eta-\mu$ [13], $\kappa-\mu$ [14], $\eta-\kappa$ [15] and $\alpha-\eta-\mu$/IGA [16] can accurately predict signal fading by properly tuning the different model parameters.

Physics-based models such as ITU-R recommended models (for instance: ITU-R P.676-13 [17]) or geometrical optics models like the Yudachev-Pinhasi model [18], although rigorous and accurate, rely on large amount of atmospheric data and variables, hindering their code implementation.

Wavelength-dependent atmospheric transmittance has been modeled using radiative transfer codes such as MODTRAN [19] and LOWTRAN [20], which account for molecular absorption by water vapor, oxygen, and carbon dioxide, as well as aerosol scattering.

More recent efforts have focused on the comparison and validation of atmospheric models under quantum communication-specific figures of merit. For instance, [21] analyzed the performance of a ground-to-satellite QKD system under varying atmospheric conditions. Similarly, [22] conducted a systematic study of the parameter space governing satellite QKD performance, identifying elevation angle, ground station aperture, and turbulence strength as the dominant factors.

Our model, although less rigorous than the models found in the literature, conforms a simple and straightforward method to obtain the atmospheric refractive index and probability of transmission, that is, a compromise between accuracy and simplicity, focusing on the latter.

Theoretical Framework and Model Description

Composition of the atmosphere

The main obstacle between Earth’s surface and satellites is the atmosphere. The atmosphere is a gas ensemble that surrounds the Earth composed of Nitrogen (78%), Oxygen (21%) and many other gasses such as Argon, Carbon dioxide, Ozone, Water vapor and Helium. We can better describe the atmosphere by identifying five different layers: Troposphere, Stratosphere, Mesosphere, Thermosphere and Exosphere. 1) Troposphere: the lowest layer of the atmosphere and the place where the Earth’s weather takes place. It extends to $\sim10$ km of altitude and the temperature decreases the deeper we reach. 2) Stratosphere: sits just above the Troposphere. It

extends from the previous layer to around 50 km and is the home of the Ozone layer. Due to the absorption of UV rays, the temperature increases with height. 3) Mesosphere: third highest atmospheric layer reaching a height of 80 km. Some nitrogen-based molecules get weakly-ionized as a result of the presence of the Sun, causing a de-ionization at night. The temperature decreases with height. 4) Thermosphere: characterized by an increase of the temperature with height and a high density of ionized gases. It extends to $\sim1000$ km. The line of 100 km of height is named Karman’s line and marks the starting point of the space. It is comprised of extremely low-pressure gasses such as Hydrogen and Helium. 5) Exosphere: the furthest atmospheric layer. It is comprised of extremely low-pressure gasses such as Hydrogen and Helium.

We need to resort to the refractive index of each layer to study how light interacts with the atmospheric layers and how much light gets absorbed by each layer, but first we need to know the pressure profile $P(z)$ and temperature profile $T(z)$ with respect to altitude $z$. There are many low-altitude atmospheric models [23], [24] that model pressure and temperature as piecewise functions, so the natural step is, starting from these models, to be able to model $P(z)$ and $T(z)$ as piecewise functions that properly describe the measured pressures and temperatures. We sampled the real, measured values from the National Oceanic and Atmospheric Administration (NOAA) and approximated the values to follow simple mathematical relations. For the temperature profile $T(z)$, we get

and, for the pressure profile $P(z)$, we arrive at

where $T(z)$ is expressed in degrees Celsius. Figure 1 plots these functions, where we can clearly see the behavior of pressure and temperature in the different atmospheric layers. From this, we may compute now the densities of the gasses in the atmosphere. Assuming that the atmospheric gasses behave as ideal gasses, then it is straightforward to check that their density can be computed as [25]

where $M_i$ stands for the molecular mass of each specie, $\chi_r$ stands for the molar fraction of the gasses. For water vapor, we have a similar expression [26]

where $h$ measures humidity and $P_s\equiv P_c\cdot\exp\left[\tau\left(\alpha_1\theta+\alpha_2\theta^{1.5}+\alpha_3\theta^3+\alpha_4\theta^{3.5}+\alpha_5\theta^4+\alpha_6\theta^{7.5}\right)\right]$ with $\tau\equiv T_c/T$, $\theta\equiv1-T/T_c$ and $T_c$, $P_c$ the temperature and pressure at the critical point of H$_2$O.

Since pressure $P(z)$ behaves as a decaying exponential, then the various densities must also obey an exponential law with altitude due to the contribution of $P(z)$ on (3)–(4). The densities computed will allow us to define how the atmosphere absorbs light on the following sections.

Atmospheric refractive index

It is known that, for a gas mixture, its refraction index must depend on the densities of the gasses that make up the total gas and on the wavelength of the incident light. As we previously noticed, the densities of the atmospheric gasses inherit the exponential behavior of atmospheric pressure, allowing us to describe the atmospheric refractive index as a function dependent on $P(z)$. We also know that, for high altitude $z$ ($P(z)\to0$), the atmospheric refractive index must tend to the value of 1 since, for high altitudes, we can approximate the atmosphere to the spatial vacuum.

Taking these considerations into account, we propose the following expression for the refractive index

where $P$ is the atmospheric pressure dependent on altitude, $P_0\equiv1\ \mathrm{atm}$ and $n_0(\lambda)=1+\Lambda_0/\lambda$. The parameter $\Lambda_0$ encodes the strength the wavelength’s contribution has on the total refractive index, its value will be set to $\Lambda_0=15\ \mathrm{nm}$ since it reproduces the correct earth surface’s refractive index for the wavelength range close to $\lambda=1000\ \mathrm{nm}$.

Earth’s weather takes place in the lowest layers of the atmosphere, we could incorporate a small, random coefficient $\delta_r$ as a correction $n_R(P,\lambda)=n(P,z)+\delta_r$, $\delta_r\ll n$. The parameter $\delta_r$ will exponentially decrease with altitude, since in the upper layers of the atmosphere there are no weather effects. Figure 2 shows the variation for wavelength $\lambda=1000\ \mathrm{nm}$.

Propagation, scattering and absorption

When a wave propagates through a medium, its intensity decreases. Theoretically, we can explain this phenomenon if we promote the refractive index to be a complex number: $\bar{n}=n_R+in_I$. To illustrate how the imaginary part of the refractive index affects the attenuation and absorption, let us consider an electric field $\vec{E}(t,\vec{x})$ in the one-dimensional case ($\vec{x}\equiv z$)

There exist a plethora of phenomena that affect the intensity loss of an electromagnetic wave propagating through the atmosphere such as heavy weather conditions like storms, foggy skies, cumulonimbus clouds forming high columns, which would represent a high loss source, rendering the satellite communication nonviable, and require precise weather forecast models; high-energy absorption phenomena like Compton scattering, photoelectric absorption and pair production (light-matter interactions), which require high energy ($\sim60$ keV) [27], small wavelength ($\sim10^{-11}$ m) electromagnetic waves, and pointing and tracking errors due to tiny fluctuations in the refractive index in small regions of the atmosphere. We may overcome the problems posed by these effects by assuming that we will work in the low-energy domain ($\lambda\sim10^{-9}$ m) to avoid energetic absorption effects, that the line of sight between satellite and ground station will not be blocked by weather conditions, and that the pointing and tracking errors will be either negligible or corrected via correction techniques [28]. Therefore, the two main sources of intensity loss in the atmosphere we will consider are the energy absorption by the different atmospheric compounds and the scattering of light rays due to the non-negligible size of the atoms and molecules. Let us split the attenuation component of the intensity into two parts

The absorption coefficient can be sampled from NASA measures on the opacity of the atmosphere [29]. To further simplify the empirical values, we approximated them by employing step functions for the values where the opacity dropped to $\sim0\%$, and negative Lorentzian distributions centered on the absorption wavelengths of the considered gas species within a relevant segment of the electromagnetic spectrum to model relevant atmospheric absorptions. The Rayleigh absorption coefficient can be approximated by first calculating the total cross-section $\sigma_{\mathrm{Rayleigh}}$ using the following equation [30]

where $\lambda$ is lights’ wavelength, $n_R$ the real part of the refractive index and $r_i$ the molecular radius of each gas species. Note that the term Rayleigh absorption coefficient is not correct since Rayleigh scattering is elastic, but a tag to emphasize that, due to the presence of this coefficient, yet another contribution to the beam intensity’s loss is considered. The relation between the cross-section and the absorption coefficient is

where $\rho_i(z)$ are the densities of the gasses, $M_i$ the molar mass of each component, and $N_A$ Avogadro constant. The complete factor that modulates the amplitude of the electric field $\vec{E}_0$ will be called the transparency coefficient. If the transparency is plotted against the wavelength we arrive at Figure 3, where we can clearly see at which wavelengths the atmosphere absorbs the most.

Transmittance. Reflection and Transmission Coefficients

The amount of light that travels through a medium is referred to as transmittance (or transmissivity), and if we want to know how much light gets reflected or passes through a

medium, we need to calculate the Reflection and Transmission coefficients. To calculate them, we resort to the Fresnel equations

where $J_{\perp/\parallel}^{r/t}$ is the reflected/transmitted energy flux of the perpendicular/parallel component of the wave, defined as the dot product of the Poynting vector $\vec{S}$ with the direction of propagation of the wave. The closer the transmission coefficients are to the value of 1.0, the more intensity passes through a medium. In this work, the computation of the transmission coefficients will enable us to analyze what percentage of photons passes through the atmosphere under any conditions at a given altitude $z$ and wavelength $\lambda$ by resorting to equation (5). The next task is addressing light polarization, which is a key part in QKD protocols, and the laser’s behavior.

QKD protocol E91 and depolarization

There are a plethora of quantum applications and protocols that could be tested in order to validate the model and to create a framework with which to evaluate viability. Our simulator had already implemented a working E91 [31] protocol. This presented an opportunity to compare against real world data taken from satellite entanglement distribution, as exposed in the paper by Yin et al. [32]. Specifically, this work employs the Micius satellite as a middle point to distribute entanglement between two ground stations 1,120 kilometers apart. Even though in the aforementioned work BBM92 [33] is used, owing to the fact that both protocols perform entanglement distribution to generate the key, we consider that the comparison should yield similar metrics in terms of QBER and Secret Key Rate (SKR). Because of this shared core principle and provided that our experiments were designed to match the exact parameters of this reference setup, with as much precision as Netsquid [34] (our simulation tool) allows, we argue that simulation results in this context should adjust to the empirical data and as such, should validate our model. In regards to the specific Netsquid requirements, we used the depolarizing channel since depolarization is a major impediment to E91, the reason lying in the nature of the protocol. When measuring with their respective bases (of which only two coincide), depolarization produces a loss of correlation or anti-correlation that tampers with expected QBER values, preventing adequate identification of eavesdroppers in the channel. For this depolarizing channel model Netsquid requires a depolarization rate ($R_{\mathrm{pol}}$) that introduces a probability of depolarization whose expression is deduced from the Lindblad master equation. The complete Lindblad Master equation reads ($\hbar= 1$)

where $R_{\mathrm{pol}}$ is the depolarization rate, $\rho$ is the density matrix, $L_i$ are the jump operators, $H$ is the time-independent closed system’s Hamiltonian and $[~ ,~]$, $\{~ ,~\}$ stand for the commutator and anticommutator, respectively. By choosing Pauli matrices as our jump operators $L_i=\sigma_i$, the interaction term on (13) becomes

since $\sigma_i^\dagger=\sigma_i$ and $\sigma_i^2=\mathbb{I}$ where $\mathbb{I}$ represents the $2\times2$ identity matrix, the anticommutator becomes $\{\mathbb{I},\rho\}\equiv2\rho$. Now, we employ the relation $2\mathbb{I}-\rho=\sum_i\sigma_i\rho\sigma_i$, and (14) becomes

Equation (15) is a first order, ordinary, differential matrix equation which can be solved by employing the following change of variables to reduce it to the standard von Neumann’s form for closed systems

and, by applying the properties of commutators on the new variables $[H,\rho]=\exp(-R_{\mathrm{pol}}\cdot t)\cdot[H,\bar{\rho}]$, we arrive at

which can be easily solved on the Schrödinger picture as

where $U(t)=\exp(-iHt)$ stands for the time-evolution operator. Finally, to solve for $\rho(t)$, we must apply the inverse change of variables noting that $\bar{\rho}(0)=\rho(0)-\mathbb{I}/2$, yielding the final expression

Therefore, the depolarization will trigger under the influence of a decaying exponential which depends on $R_{\mathrm{pol}}$. This exponential is the method Netsquid uses to evaluate the depolarization probability. Notwithstanding this effort, we shall as well justify which depolarization probability we chose. This depolarization probability can be deducted by the polarization extinction ratio (PER), this figure is defined as the ratio of powers of different polarizations. It characterizes the degree of polarization of the devices. To match the PER that should be expected in the Micius experiment, we referred to the work authored by Wu et al. [35] analyzing the optical elements pertaining to the Micius experiment. Defining an expected QBER due to depolarization that essentially matches the depolarization probability ($[\mathrm{QBER}_{\mathrm{pol}}\sim1-\exp(-R_{\mathrm{pol}}\cdot t)]$), allows us to, from a reference PER, obtain the depolarization rate. Essentially $\mathrm{QBER}_{\mathrm{pol}}\approx1/\mathrm{PER}$, this simplification is explained in their work and relies on the fact that the minimum required PER is high valued with respect to 1, and that deviations between intended and real polarization values are negligible with the equipment that is used. Depending on the delay the beam experiences the actual probability and thus, depolarization rate varies. The variations of the value of $R_{\mathrm{pol}}$ near the visual range of the electromagnetic spectrum is small enough to consider $R_{\mathrm{pol}}$ as a constant and, as a general figure and due to time constraints on the code implementation, we took $R_{\mathrm{pol}}=8$ Hz, which corresponds to higher depolarization probabilities than those experienced in the travel distances in the orbit points in between ground stations (corresponding to $5.2$ Hz or $2\%$ depolarization probability approximately) which, in our opinion, balances the optimistic predictions the model asserts and allows us to treat transmissions within an expected conservative depolarization range throughout all simulations. Additionally, since the value of $\mathrm{QBER}_{\mathrm{pol}}$ exponentially depends on $R_{\mathrm{pol}}$, small changes on the value of the depolarization rate will not trigger significant changes on the QBER values. We consider that this choice balances the optimistic predictions the model asserts and allows us to treat transmissions within an expected conservative depolarization range throughout all simulations.

Theoretical model for the laser: Gaussian Beams

The easiest way to model a real light beam is using the Gaussian Beam description. A Gaussian beam is described by an electric (magnetic) field whose amplitude envelope in the transverse plane is given by a Gaussian curve. One of the most important consequences of this profile is the fact that the cross-section of the laser does not increase in a linear fashion, instead the cross-section decreases until a minimum section called waist $w_0$ and then starts increasing again. The equations for the Gaussian beam traveling in the $z$ direction are [36]

where $w(z) = w_0(1 + (z/z_0)^2)^{1/2}$, $z_0 = \pi n_R w_0^2/\lambda$, $R(z) = z(1 + (z_0/z)^2)$, and $\xi(z) = \arctan(z/z_0)$, and, the intensity of the beam is given by

$\eta_w$ is the wave impedance of the medium in which the beam is propagating.

Optical link power model

We turn now our attention to the satellites used in the communication. The usual optical link can be modeled by an equation, a product of efficiencies and losses that relates transmitted power to the received power, this is a modification of the Friis equation for the link budget, the specific modifications applied for the case of optical communications is extracted from Hemmati’s book on Deep Space Optical Communication [37], as for all expressions in this section

where $\eta_{(\cdot)}$ stands for the efficiency of the transmitter’s optics, aperture illumination and receiver’s collector respectively, $A_{(\cdot)}$ is the aperture are for transmission and reception and $L_{(\cdot)}$ stands for the different sources of power loss. The transmitter and receiver gains can respectively be defined as $G_T = 4\pi A_T/\lambda^2$ and $G_R = 4\pi A_R/\lambda^2$; the terms $A_T, A_R$ are the apertures of both the transmitter and receiver.

In terms of space communications, augmenting the transmitted power could be a good idea to increase the received power, but it would be restricted to the power budget requirements. Another idea could be increasing the aperture area, which would reduce the beam’s width but, in return, increases the probability and magnitude of the PAT error.

Since we are working with the gaussian beam formalism, the expression for the aperture illumination efficiency is

where $\gamma$ refers to the obscuration ratio, $S$ is the RMS optical path difference (for smooth optics, its value is 28% of the peak to valley value), $\alpha\equiv A_d/2w(z)$ where $A_d$ corresponds to the aperture diameter and $w(z)$ is the beam’s radius.

Background interference sources

Most of the interference can be boiled down to diffused light from the atmosphere ($P_{\mathrm{diffuse}}$), or sourced from the planets ($P_{\mathrm{planetary}}$), scattered light in the optics ($P_{\mathrm{stray}}$). These can be expressed as

where $H(\lambda)$ is the spectral irradiance of the background, $L_\lambda$ refers to the sky irradiance, $I_\lambda$ is the exo-atmospheric solar constant (0.074 W/cm$^2 \cdot\mu$m), $BSDF(\theta)$ is the Directional Scatter Distribution, which follows a power law distribution as a function of $\theta$ and $T(\theta)$ being the atmospheric absorption.

Specifications and Technical Details

To correctly address design goals and measure accurate results that validate the model, it is important to create a reference framework. Testing specific cases and results would only be accurate if the model tried to replicate the conditions of a certain set of experimental results. Of course, a certain abstraction level is needed to reduce or neglect some variables or take into account others. In terms of satellite quantum communications we have to take into account some limitations when designing a framework. First, technical standards do not meet established legal or commercial regulations, neither is there a specific standard that is commonly and internationally accepted. This is unlike other telecom services where specific power, precision and wavelengths would exist. Second, each country with existing quantum communication projects may have different providers, goals and requirements, so two completely different sets of parameters can be completely valid because they aim to fulfill different goals or different conditions. Thirdly, there are some countries that do have such a significant presence in the space sector, that their technologies and standards can get imposed as a side effect, with no guarantee that other future projects will adopt it. Finally, the actual constraints of the different subsystems of a space system vastly differ in the case of a ground station and the satellite by itself. This asymmetry is key to correctly understand the nature of the communication links.

Optical Space Segment Specifications

Space segment specifications represent the maximum constraint that the link reaches, as the power budget in satellites and deployed systems are strict. Apertures and efficiencies in the satellite are actually the most important parameters to take into account because the modules have to fulfil a given set of requirements, which conditions shapes, mass, volume, materials, and many other variables. As such, one cannot simply assume an aperture of an arbitrary number of meters, because for that arbitrarily large telescope to be orbiting there would need to exist appropriate equipment to actually deploy it. Looking at QKD projects, in the case of the Micius satellite a 300-mm-aperture Cassegrain telescope was used and a 5 $\mu$W laser was used for tracking, the QKD link was operating at 850 nm. The diffraction loss was estimated to be 22 dB for a link of 1200 km, the tracking precision was 1.2 micro radians [38]. The Micius satellite is actually a relatively large (over 500 kilograms) satellite in comparison to other initiatives. In contrast, the QUBE-II [39], which is based on the CubeL mission (2021), for CubeSIL and the precursor mission QUBE-I (2024) which follows the CUBESAT standard [40], it represents a radical shift in paradigm compared to Micius. Firstly, CubeSats operate in different standardised categories that are internationally respected, from U1 (10 × 10 × 10 cm for 1 to 1.33 kg) to U6 or U12 (20 × 20 × 30 cm for 20 to 50 kg). The difference between categories is remarkable, and actually, U6 U12 CubeSats are in fact deployed for commercial use with specialized modules. The QUBE-I operates at 850 nm and has an aperture of 50 mm. Larger satellites can host more equipment and are relatively less constrained in terms of power. As a result, the non-standardized satellites are not only more technologically advanced and thus less restricted, but this develops in such an advantage over CubeSats, that large satellites dedicated to QKD can be deployed to geostationary orbits.

SYSTEM PARAMETERS
Orbit	GEO
Position	30° W – 61° W
Power	258 W (TBC)
Mass	212 kg (TBC)
QKD Tech.	Discrete Variable
Protocol	BB84 + Decoy
Link Budget	47 dB (TBC)
Wavelength	1550 nm
Telescope aperture	50 cm
Key Rate	30–500 bps
Key Volume	1.5–3 Gb
Pointing Accuracy	$1 \,\mu\mathrm{rad}$
Services	Commercial

This is a great advantage in terms of precision but it is not attainable for all latitudes. This is because after a certain lower bound elevation angle, the distance that the photons need to travel over the atmosphere is so large that the transmitivity vanishes. This ”atmospheric wall” effect supposes a great disadvantage for GEO satellites because they inherently need to orbit in the Earth’s plane. This means that for many regions of the planet direct line of sight (LOS) is heavily impeded for GEO quantum communication. For example, a 30 mm aperture transmitter/receiver in a CubeSat supposed to operate in GEO is not realistic because it would not be able to compensate the ”atmospheric wall” effect. This could be easily overcome by running the model and measuring the received power. For LEO and MEO satellites, using apertures like 300 mm or more can be unrealistic, especially for LEO. While for bigger apertures, better yields are expected, the satellites that are usually deployed in these conditions are restricted in terms of space, shape, and volume. In other words, when working with a CubeSat of a given category, then we know that the aperture does not only need to be proportional to the size of the CubeSat specification, but that it needs to fit its size. For satellites developed and operated by government actors, these constrains are less of an issue. For instance, the Hispasat QKD-GEO satellite project aims for GEO QKD with the requirements seen in Table I.

Terrestrial Network Integration

Integration with terrestrial networks needs to balance the use of the different wavelengths that are usable in the context of QKD. As the wavelength increases in value, the required precision to operate also increases. This means that the greater the wavelength (lower frequencies) the higher the precision the equipment must attain. This is caused by the dispersion of the laser. Recalling the expression of the Gaussian Beam, dispersion depends on the value of $z_0$ which constrains the waist of the beam and how it propagates onward. To better understand this, Figure 4 depicts the two main approaches to quantum communications. Basically, smaller wavelengths have a greater divergence angle, needing more precise optical equipment, Equipment becomes precise after the implementation of advanced ATP systems (Advanced Tracking and Pointing System), to levels below the micro radian.

Equipment produced to work in the 800 nm (although not this exact wavelength but approximated) is far more efficient and cheap. The 1550 nm wavelength is however, much better suited to work in existing fiber optic systems. Thus, it is a matter of reconciliation between network domains. In summary, free space quantum communications works in the 800 nm mark because of the efficiency and low power consumption of the emitters compensating the higher dispersion rate with technology such as ATP systems, this tradeoff is beneficial but has the side effect of isolating free space networks from ground networks operating with fiber optics in terms of wavelength. This interfacing problem is not an issue. Ground Stations can exchange information with both ground and space network (at the cost of slight increase in complexity).

Our model reflects this tradeoff between the wavelength and the precision. For each actual scenario the designer would have to take into account the expected values of dispersion of the beam and the exact wavelength the equipment is working on.

Model Results

Atmospheric test for Micius and Hispasat missions

Based on the previously discussed conditions for the different projects we decided to test the model results by simulating mainly two specific cases, the Hispasat and the Micius satellite. On one hand, the Micius satellite, with abundant information of technical details and results is the ideal case to reproduce for our model. Comparing the results of the Micius experiment and the result of our model provides a unique opportunity for validation, as previously discussed. The Hispasat satellite does not have such an enormous base of scientific literature and measurements, but poses a similar initiative yet at the same time radically distinct approach (due to the changes in altitude, wavelength and receiver and transmitter specifications). Based on this premise, we have produced a series of figures as the output of the simulation based on the developed model, for which we have taken into account both the downlink and the uplink for each of the cases (Micius and Hispasat).

For the Hispasat case a wavelength of 1550 nm was used as indicated by the requirements of the project, in the case of the Micius case, a 810 nm wavelength, in both cases the same ground station was employed (1 m in aperture). Despite both GS having different specifications for this parameter we wanted to use a 1 m aperture as a compromise between both cases, by not using bigger apertures we aimed to produce more conservative results. The satellites’ apertures are sampled from the original works regarding the Hispasat and Micius missions [41, 1], whose values are 50 cm and 30 cm, respectively.

Figure 5 shows the previously mentioned atmospheric wall effect, by which after a certain elevation angle in the uplink transmission, the initial segment of the communication experiments huge losses as the beams approach the horizon. Downlink losses are not nearly as severe for these extremes angles because relative to the total distance of the link, the final segment (the atmosphere in that case) is minuscule. Still, as the angles reach 90 degrees, the losses rise asymptotically.

Beyond basic power losses we can separate this result into other variables to see exactly why the model behaves like this. Figure 6 shows how transmissivity evolves for both Micius and Hispasat, the wall effect is more remarkable here for the Micius case while Hispasat manages to hold better transmissivity values across the collapsing region which supports the idea that the requirements and the concept of the mission conveyed, Hispasat’s QKD-GEO aimed for GEO quantum communications, if we take into account the latitude of a country such as Spain (which is much closer to the equator compared to other Northern European ones) and the fact that the transmissivity holds for longer in this case we can deduce that such a communication scheme is valid , even though Hispasat shows a lower transmissivity in general to that of the Micius satellite. In contrast, downlink does not show much variation for the same reasons as the loss of power did.

In the case of the losses that the atmosphere shows, in Figure 6 (right), these are a perfect negative image to that of the power loss. Specially so because in our model’s simulations we have not taken into account statistical losses for other phenomena not relative to the atmosphere, however, the difference in atmospheric losses between the cases are not as remarkable as in the power losses. This could be explained by the fact that in both cases we are operating at wavelengths at which the atmosphere behaves very similarly. The wall effect is also visible here and the uplink/downlink asymmetry is also present.

These results can be contrasted to the predictions of the Micius experiments. In terms of diffraction, approximately 22 dB were calculated, and absorption in atmosphere losses ranging from 3 dB to 8 dB. Our predictions for atmospheric losses in the elevation angles that Micius operated in (10° to 15°) our model predicts around −3 dB, that is coherent, specially taking into consideration that in our model there are other phenomena that were not taken into account, many contributions that despite being small individually, may amount to compensate with this offset, which despite existing, is within a reasonable margin of error. Total losses were calculated to be in the range −31 dB to −39 dB aproximately. Our model’s prediction stays at −29 dB. Again, taking into account that the model lacks certain phenomena, and that the conditions for the simulation were similar but not identical so as to be compared with the hypothetical Hispasat case, we think that the predictions are in fact reasonable. Figure ?? specifically targets the range that the Micius satellite works in, with the exact parameters that the Micius has. The difference is now much higher (the model predicts a −27 dB to −30 dB loss), and this is due to this lack of modeled phenomena. Results still approximate to the correct range, so while concrete and precise values are not possible with the actual state of the model, we are certain that with more refined additions to the code these could be achieved. In terms of the transmissivity according to the atmosphere losses, with absorption losses of −3 dB to −8 dB, that implies a transmissivity in the range of 0.5 to 0.16, our prediction stays at the 0.5 to 0.46 range approximately, so in terms of transmissivity we conclude that the predictions are close to experimental results.

Link capacity based on model results

Power budgets and transmissivity are a helpful approach towards the exploration of viable quantum communication link, but the analysis beyond viability pertains more expressions and relationships beyond this. For instance, capacity establishes a limit to the size of the possible key established by quantum key distribution protocols (QKD). For the capacity we can use the Pirandola–Laurenza–Ottaviani–Banchi (PLOB) bound [42] $C(T) = -\log_{2}(1-T)$, where $T$ is the total channel transmissivity (Figure 5, right). In our case this means we will be dealing with the transmissivity of the atmospheric link, which is generated by the simulation model. In terms of longer chains of repeaters, such as multi-satellite communication networks this can simply be reduced to the minimum capacity link, meaning that low yield links will produce a bottleneck in the network reducing the upper bound of the key significantly. This, together with our observations of the model results in the Uplink capacity being the bottleneck (so long as there are not worse links once in space, which is unlikely).

Satellite QKD comparison with Micius

To assess the practical viability of E91-based quantum key distribution (QKD) under realistic atmospheric conditions, we performed a series of numerical simulations to reproduce the scenario entanglement distribution scenario between the two ground stations and satellite involved in the aforementioned Micius experiment using our atmospheric channel model discussed on Section III.

Micius experiment using our atmospheric channel model discussed on Section III.

Simulations were carried out for elevation angles ranging from ($\varepsilon_1 \in[120^\circ, 10^\circ]$) measured from one of the ground stations, with a satellite altitude of 500 Km and a wavelength of $\lambda= 810$ nm, consistent with the LEO configuration employed in the Micius satellite experiments [1]. For each elevation angle we computed the the expected channel transmittance, the resulting quantum bit error rate (QBER) and the sifted key rate following the standard E91 security analysis. Table II presents the corresponding Secure Key Rate in bits per second and QBER of said angles.

Elevation (°)	SKR	SKR Std.Dev.	QBER	QBER Std.Dev.
117	0.30	0.10	0.010756	0.004940
85	0.35	0.14	0.007855	0.004941
54	0.20	0.06	0.008672	0.005117
36	0.24	0.09	0.007482	0.004962
24	0.32	0.14	0.008358	0.004359
17	0.28	0.09	0.009408	0.005044

The QBER obtained from our simulations remains below the 10% threshold, allowing for sufficient preemptive margin for eavesdropper detection. In contrast with real Micius data, our results predict lower QBER (at an average of 1% error rate), but for the reference real values that were measured to have that level of QBER the SKR relies close below the values, thus being more pessimistic. We performed 200 repetitions for each point.

The disagreements between the experimental results and our simulated values are small in absolute value and originate from our assumption that the satellite remains in the same position at the instant of time of the entanglement distribution for a given transmission of a key block (which is needed in order to use Netsquid, which does not incorporate dynamic systems), and consequently pointing and tracking errors are not realistic, which is an unavoidable limitation of static simulators. A noteworthy observation is that the channel shows more SKR when the satellite is between ground stations, in terms of transmissivity both uplink and downlink have a similar value in this region, amounting to a higher probability that both photons arrive at their destination. However it is also true that because of the slower fall in transmissivity as a function of the elevation angle in the downlink channel, that also has the effect of increasing SKR in positions that are better for the uplink, as that does not lower the downlink’s transmissivity as much. Overall, these results support the conclusion that E91-based QKD over a satellite link is numerically consistent with the performance demonstrated by the Micius experiment, and validate our atmospheric model as a reliable and open tool for predicting the feasibility of future satellite-based QKD missions.

Starlink exploratory assessment

The Starlink constellation comprises LEO satellites at altitudes ranging from 340 km to 570 km, with orbital inclinations designed to provide near-global coverage. Equipped with the results provided by the atmospheric channel model, we will study how feasible quantum communications employing the Starlink constellation are. We must have in mind for the purposes of simulation, that we considered similar satellite specifications as the Micius mission within significant changes on altitude (550 km from Earth’s surface), wavelength ($\lambda= 1550$ nm) and elevation angles ($\theta\in[20^\circ, 70^\circ]$) this due to the lack of existing reference infrastructure for quantum communications in the Starlink constellation, which means: there is no existing quantum communication infrastructure, and existing infrastructure for classical communication cannot be used as a reference. We performed simulations to obtain the atmospheric transmissivity and power balance to be able to plot the computed values against the different elevation angles $\theta$ to study the optimal window of operation, similar to previous cases. The simulated results can be observed in Figure 7 where we can clearly notice that, if such a Starlink network is developed, it would need to respect this technological

restrictions imposed by the specifications of the Micius ground station and the computed received power. In other words, mission requirements for a quantum communication network deployment would need to take into account how the received power decreases with the inclination angles.

The restriction in terms of throughput for these networks and megaconstellations could be severely limited depending on the particular protocol, but also in terms of the endpoints of the communication. Such low throughput (inferior to 1 bit/s) as seen in the validation phase of this work, corresponds to entanglement distribution networks where satellites and intermediaries between ground stations. As such, the effective probability of such a communication to be successful can be roughly described as the probability of success of both photons to arrive correctly at their corresponding ground stations. If one is lost during transmission, the pair is effectively lost. Notwithstanding these effects for ground to space quantum communications, it is true that performing ground to space key generation, or in-space quantum applications, would allow for far greater throughput, moreover, constellations and megaconstellations would allow for near perfect channels where only pointing and tracking errors worsen the rate of transmission. With highly populated megaconstellations these errors can be greatly reduced as the distance between satellites decreases. Providing a promising framework for quantum communications in space.

Conclusions

In this work, we have developed a theoretical framework to predict results regarding transmissivity, atmospheric and total power losses for satellite quantum communication scenarios and compared against the Micius satellite for validation and HıspaSat’s QKD-GEO as a way to analyze the viability of the established requirements. We used publicly available information to program our model. Our results match the observations from Micius and Hispasat, demonstrating that the model adjusts to experimental outcomes.

Acknowledgments

This work was supported in part by MCIU/AEI/10.13039/501100011033 (DISCOVERY Project) under Grant PID2023-1487160B-C31 and by the Galician Regional Government under Project ED431B 2024/41 (GPC).