# Numerical simulations of atmospheric quantum channels

M. Klen  Bogolyubov Institute for Theoretical Physics, NAS of Ukraine,
Vulytsya Metrologichna 14b, 03143 Kyiv, Ukraine  Department of Theoretical and
Mathematical Physics, Kyiv Academic University, Boulevard Vernadskogo 36,
03142 Kyiv, Ukraine  A. A. Semenov  Bogolyubov Institute for Theoretical
Physics, NAS of Ukraine, Vulytsya Metrologichna 14b, 03143 Kyiv, Ukraine
Department of Theoretical and Mathematical Physics, Kyiv Academic University,
Boulevard Vernadskogo 36, 03142 Kyiv, Ukraine  Institute of Physics, NAS of
Ukraine, Prospect Nauky 46, 03028 Kyiv, Ukraine

###### Abstract

Atmospheric turbulence is one of the lead disturbance factors for free-space
quantum communication. The quantum states of light in such channels are
affected by fluctuating losses characterized by the probability distribution
of transmittance (PDT). We obtain the PDT for different horizontal links via
numerical simulations of light transmission through the atmosphere. The
results are compared with analytical models: the truncated log-normal
distribution, the beam-wandering model, the elliptic-beam approximation, and
the model based on the law of total probability. Their applicability is shown
to be strongly dependent on the receiver aperture radius. We introduce an
empirical model based on the Beta distribution, which is in good agreement
with numerical simulations for a wide range of channel parameters. However,
there are still scenarios where none of the above analytical models fits the
numerically simulated data. The numerical simulation is then used to analyze
the transmission of quadrature-squeezed light through free-space channels.

##  I Introduction

Quantum communication is a rapidly developing field with a wide range of
applications. One of its main advantages is the ability to establish secure
links between remote parties (for a review see Refs. [1, 2, 3, 4]). It can
also be used in many other applications, such as quantum digital signature
[5], connecting quantum devices with quantum teleportation [6, 7], and
entanglement swapping protocols [8], to name just a few.

Optical radiation fields serve as natural carriers of quantum information for
communication channels. Quantum light can be distributed through optical
fibers or free space. The latter offers several practical advantages,
including the ability to establish satellite-mediated global quantum
communication, communication through hard-to-access regions, and communication
with moving parties. Various implementations of free-space channels have been
reported in the literature for horizontal (see, e.g., Refs. [9, 10, 11, 12,
13, 14, 15, 16, 17, 18, 19]) and vertical (see, e.g., Refs. [20, 21, 22])
links, including experiments involving satellites [23, 24, 25, 26, 27, 28, 29,
30, 31, 32, 33, 34, 35, 36].

Atmospheric turbulence is one of the major obstacles for free-space quantum
communication. Its effect on the optical fields is described by methods [37,
38, 39] developed in classical atmospheric optics. The corresponding
techniques differ depending on the degrees of freedom of light and the types
of measurements involved in the protocol. For example, protocols dealing with
optical angular momentum [40, 41, 42] require a proper description of the
spatial structure of the modes [43, 44, 45, 46, 47, 48, 49]. Other types of
protocols and experiments involving photocounting [50, 51], polarization
analysis [9, 12, 52, 53], homodyne detection [10, 11, 54, 55, 16, 56], and so
on, demand a proper description of the fluctuating losses in atmospheric
channels. In this case, the light mode at the transmitter can be prepared, for
example, as a quasi monochromatic mode in the form of a pulsed Gaussian beam.
For discrete-variable protocols, using propagation of single photons in each
direction and/or measurements with long detection times, fluctuations of
losses (scintillation) can be neglected [57, 58, 59].

In this paper we consider communication protocols with a fixed spatial
structure of light modes at the transmitter. As mentioned, in this case the
quantum states of a quasi monochromatic light mode with a fluctuating shape at
the receiver are affected by fluctuating losses. They change the quantum
states of the quasi monochromatic mode according to the input-output relation
(cf. Refs. [60, 55, 61] and Appendix A),

$\displaystyle
P_{\mathrm{out}}\left(\alpha\right)=\int\limits_{0}^{1}d\eta\frac{1}{\eta}P_{\mathrm{in}}\left(\frac{\alpha}{\sqrt{\eta}}\right)\mathcal{P}(\eta).$
(1)

Here $P_{\mathrm{out}}(\alpha)$ and $P_{\mathrm{in}}(\alpha)$ are the Glauber-
Sudarshan $P$ functions [62, 63] at the channel output and input,
respectively, $\eta$ is the channel transmittance, and $\mathcal{P}(\eta)$ is
the probability distribution of transmittance (PDT). Obviously, the PDT is the
main characteristic of quantum channels. Equation (1) has a similar form in
both quantum and classical optics and thus the PDT can be obtained from purely
classical considerations.

The mentioned approach has been successfully applied to analyze quantum-state
transfer through free-space channels (see Ref. [64] for a brief summary). For
example, quadrature squeezing [55, 61], Bell nonlocality [52, 53], Gaussian
entanglement [65, 66], and higher-order nonclassical and entanglement
phenomena [67, 68] have been analyzed in the context of their transfer through
atmospheric channels. Also the input-output relation (1) has been applied to
analyze quantum communication protocols such as the decoy state protocol [69,
70], continuous-variable quantum key distribution protocols [54, 71, 72, 73,
74, 75, 76, 77, 78, 79, 80, 81], and continuous-variable quantum teleportation
[82, 83, 84].

A rigorous description of free-space quantum channels requires knowledge of
the PDT. Random fluctuations of the refractive index due to atmospheric
turbulence lead to random wandering of the beam centroid, fluctuations of the
beam-spot shape, etc. These effects in turn lead to fluctuations in the
intensity of the beam passing through the receiver aperture. The PDT
characterizes the quantum channel associated with a quasimonochromatic mode.
This characteristic can be obtained from purely classical considerations,
corresponding to coherent states at the transmitter.

Several analytical models of the PDT have been proposed. One of them focuses
on the effect of beam wandering [61]. Furthermore, the elliptic-beam model
[85] involves beam-spot distortions approximated by Gaussian elliptic beams
with randomly oriented semiaxes. In particular, this model shows that beam-
spot distortions can be effectively described by the truncated log-normal
distribution considered earlier in Refs. [50, 51]. Assuming that beam
wandering and broadening are statistically independent, the law of total
probability can be employed to derive the PDT [70].

All analytical models have a common drawback: They are applicable for limited
values of channel parameters such as turbulence strength, propagation
distance, receiver aperture size, and beam parameters. It is usually difficult
to determine which model to use for given values of the channel parameters. In
this paper we address this problem by providing numerical simulations of the
PDT for horizontal links and comparing the results with analytical models. For
this purpose, we use the sparse-spectrum model [86, 87, 88] of the phase-
screen method (see Refs. [89, 90, 91, 92] and Refs. [93, 48, 49, 94, 95] for
its recent applications). We also numerically check statistical properties of
beams after passing through the atmosphere: Gaussianity of the distribution of
the beam-centroid position, correlation between beam deflection and beam
shape, etc. These properties are directly related to the applicability of
analytical models.

The phase-screen method cannot be considered a substitute for experimental
investigations, since it involves approximations. Nevertheless, the sparse-
spectrum model yields the statistics of phase disturbances that are in good
agreement with the analytical expression (cf. Appendix B). Most importantly,
the method is free of assumptions about the beam shape and its statistical
characteristics made in the analytical models of the PDT. Thus, our results
can be useful to determine their applicability. Numerically simulated data can
also be used directly to analyze quantum light transmission in atmospheric
channels for a variety of scenarios. Taking the quadrature squeezing as an
instance, we demonstrate how numerical analysis can be employed to study the
transmission of nonclassical phenomena through free-space channels.

The rest of the paper is organized as follows. In Sec. II we introduce
preliminary information on the propagation of Gaussian beams through the
turbulent atmosphere, the sparse-spectrum formulation of the phase-screen
method, and analytical models. An empirical analytical model based on the Beta
distribution, which shows good agreement with the numerically simulated data,
is introduced in Sec. III. Numerical simulations of atmospheric channels and
comparison of the obtained numerically-simulated data with analytical models
are considered in Sec. IV. In Sec. V the statistical properties of beams after
passing through the turbulent atmosphere are analyzed. An application of the
discussed method to a description of the transmission of quadrature squeezing
through free-space channels is considered in Sec. VI. A summary and concluding
remarks are given in Sec. VII. In the Supplemental Material [96] one can find
numerically simulated data for different channels, an interactive
visualization tool, and the corresponding Python 3 codes.

##  II Preliminaries

In this section we briefly introduce preliminary results needed for our
consideration. All of them are related to channel characterization by the PDT
and are thus derived from classical considerations. First, we discuss the
fluctuating channel transmittance $\eta$. Then we briefly review the sparse-
spectrum formulation of the phase-screen method used for our numerical
simulations. Finally, we review the existing analytical models of the PDT.

###  II.1 Channel transmittance

The distribution of a light beam through the turbulent atmosphere along the
$z$ axis can be described by the paraxial equation [38] for the field
amplitude $u(\mathbf{r};z)$,

$\displaystyle 2ik\frac{\partial u(\mathbf{r};z)}{\partial
z}+\Delta_{\mathbf{r}}u(\mathbf{r};z)+2k^{2}\delta
n(\mathbf{r},z)u(\mathbf{r};z)=0,$ (2)

where $k$ is the wave number, $\delta n(\mathbf{r},z)$ is a fluctuating part
of the index of refraction of air, and
$\mathbf{r}=\begin{pmatrix}x&y\end{pmatrix}^{T}$ is the vector of transverse
coordinates. The Gaussian beams are defined by the boundary conditions of this
equation,

$\displaystyle u(\mathbf{r};0)=\sqrt{\frac{2}{\pi
W_{0}^{2}}}\exp\Bigl{[}-\frac{\mathbf{r}^{2}}{W_{0}^{2}}{-}\frac{ik}{2F_{0}}\mathbf{r}^{2}\Bigr{]}.$
(3)

Here $W_{0}$ and $F_{0}$ are the beam-spot radius and the wave-front radius at
the transmitter, respectively. For our purposes we will use two types of
initial conditions. The first corresponds to the case of $F_{0}=+\infty$. In
the absence of atmosphere, this beam is collimated i.e. it is minimally
divergent, for a propagation distance much shorter than the Rayleigh length,
$z\ll z_{\mathrm{R}}=kW_{0}^{2}/2$. The second corresponds to
$F_{0}=z_{\mathrm{ap}}$, where $z_{\mathrm{ap}}$ is the channel length. Such a
beam is considered to be geometrically focused (cf. Ref. [97]) because in the
absence of atmosphere it has a minimum beam spot radius in the plane of the
receiver aperture, $z=z_{\mathrm{ap}}$, given $W_{0}$ and $z_{\mathrm{ap}}$.

Random fluctuations of the index of refraction in the case of isotropic
homogeneous turbulence are described by the correlation function, the
Markovian approximation for which reads

$\displaystyle\langle\delta n$ $\displaystyle(\mathbf{r}_{1};z_{1})\delta
n(\mathbf{r}_{2};z_{2})\rangle$
$\displaystyle=\int_{\mathbb{R}^{2}}d^{2}\boldsymbol{\kappa}\Phi_{n}(\boldsymbol{\kappa};\kappa_{z}{=}0)e^{i\boldsymbol{\kappa}\cdot(\mathbf{r}_{1}-\mathbf{r}_{2})}\delta(z_{1}-z_{2}),$
(4)

where $\Phi_{n}(\boldsymbol{\kappa};\kappa_{z})$ is the turbulence spectrum.
We will use it in the modified von Kármán–Tatarskii form given by

$\displaystyle\Phi_{n}(\boldsymbol{\kappa};\kappa_{z})=\frac{0.033C_{n}^{2}\exp\left[-\left(\frac{\kappa\ell_{0}}{2\pi}\right)^{2}\right]}{\left(\kappa^{2}+L_{0}^{-2}\right)^{11/6}}.$
(5)

Here $C_{n}^{2}$ is the index-of-refraction structure constant characterizing
the local strength of the turbulence, $\ell_{0}$ and $L_{0}$ are inner and
outer scales of turbulence corresponding to minimum and maximum sizes of the
turbulence eddies, respectively, and
$\kappa=\sqrt{\boldsymbol{\kappa}^{2}+\kappa_{z}^{2}}$.

The beam intensity is defined as $I(\mathbf{r},z)=|u(\mathbf{r};z)|^{2}$. As
discussed in Refs. [61, 85, 70] (see also Appendix A), the transmittance is
given by integration of the field intensity over the aperture opening
$\mathcal{A}$ in the receiver aperture plane $z=z_{\mathrm{ap}}$,

$\displaystyle\eta=\int_{\mathcal{A}}d^{2}\mathbf{r}I(\mathbf{r};z_{\mathrm{ap}}).$
(6)

Fluctuations of $\delta n(\mathbf{r};z)$ lead to fluctuations of
$I(\mathbf{r};z)$ and consequently to randomization of $\eta$. Moreover, this
consideration does not take into account deterministic losses caused by
absorption and scattering, which can be easily included in the final results.

Next, we consider the second- and fourth-order correlation functions given by

$\displaystyle\Gamma_{2}(\mathbf{r};z)=\left\langle
I(\mathbf{r};z)\right\rangle,$ (7)

and

$\displaystyle\Gamma_{4}(\mathbf{r}_{1},\mathbf{r}_{2};z)=\left\langle
I(\mathbf{r}_{1};z)I(\mathbf{r}_{2};z)\right\rangle,$ (8)

respectively. These functions can be used to obtain various statistical
properties of the light beam passing through the turbulent atmosphere. For
instance, the first two moments of the channel transmittance can be evaluated
as

$\displaystyle\langle\eta\rangle=\int_{\mathcal{A}}d^{2}\mathbf{r}\Gamma_{2}(\mathbf{r};z)$
(9)

and

$\displaystyle\langle\eta^{2}\rangle=\int_{\mathcal{A}}d^{2}\mathbf{r}_{1}\int_{\mathcal{A}}d^{2}\mathbf{r}_{2}\Gamma_{4}(\mathbf{r}_{1},\mathbf{r}_{2};z).$
(10)

Other applications of these functions needed for analytical models of the PDT
will be discussed in Sec. II.3. In this context, it is also important to
mention the recent results in Refs. [94, 95] about spatial modes that are
different from Gaussian beams and which exhibit better transmittance for
particular realizations of turbulent channels.

###  II.2 Sparse-spectrum model for phase-screen method

According to Eq. (6), in order to sample the channel transmittance $\eta$, one
must first sample the field amplitude $u(\mathbf{r};z)$. This can be
accomplished using the sparse-spectrum model for the phase-screen method [86,
87, 88]. Its verification technique is discussed in Appendix B. The method
relies on the fact that given the field amplitude $u(\mathbf{r},z_{m-1})$ at
the point $z=z_{m-1}$, one can find an approximate solution to the paraxial
equation (2) at the point $z=z_{m}>z_{m-1}$ up to the second order of
$l=z_{m}-z_{m-1}$ (cf. Ref. [89]),

$\displaystyle u(\mathbf{r};z_{m})\approx
e^{-({il}/{4k})\Delta_{\mathbf{r}}}e^{-i\phi(\mathbf{r})}e^{-({il}/{4k})\Delta_{\mathbf{r}}}u(\mathbf{r};z_{m-1}).$
(11)

Here

$\displaystyle\phi(\mathbf{r})=k\int\limits_{z_{m-1}}^{z_{m}}dz\,\delta
n(\mathbf{r},z)$ (12)

is referred to as the phase screen. As follows from Eq. (II.1), the
correlation function of the latter is

$\displaystyle\langle\phi(\mathbf{r}_{1})\phi(\mathbf{r}_{2})\rangle=\int_{\mathbb{R}^{2}}d^{2}\boldsymbol{\kappa}\Phi_{\phi}(\boldsymbol{\kappa})e^{i\boldsymbol{\kappa}\cdot(\mathbf{r}_{1}-\mathbf{r}_{2})},$
(13)

where

$\displaystyle\Phi_{\phi}(\boldsymbol{\kappa})=2\pi
lk^{2}\Phi_{n}(\boldsymbol{\kappa};0)$ (14)

is the power spectral density of the phase screens. Equation (11) implies that
propagation between the points $z_{m-1}$ and $z_{m}$ is reduced to vacuum
propagation up to the midpoint of this interval, followed by the phase
increment by $\phi(\mathbf{r})$ and finally vacuum propagation again.

The phase-screen method can be outlined as follows. First, the propagation
distance is divided into $M$ slabs $[z_{m-1},z_{m}]$, where $m=1,\ldots,M$,
$z_{0}=0$, and $z_{M}=z_{\mathrm{ap}}$. Then the phase screens are sampled at
the midpoints of these slabs. Finally, the vacuum propagation of the field
amplitude $u(\mathbf{r},z)$ is simulated successively between these points and
the $\mathbf{r}$-dependent phase is incremented on them.

The vacuum propagation can be simulated using the Fourier transform of the
field amplitude. In the sparse-spectrum model [86, 87, 88], the phase screens
are expanded into $N$ harmonics such that

$\displaystyle\phi(\mathbf{r})\approx\mbox{\rm
Re}\sum_{n=1}^{N}a_{n}e^{i\boldsymbol{\kappa}_{n}\cdot\mathbf{r}}.$ (15)

The sampled values of $\phi(\mathbf{r})$ are obtained via sampling the
expansion coefficients $a_{n}$ and the random vectors
$\boldsymbol{\kappa}_{n}$. Their statistics should satisfy Eq. (13). This can
be achieved if two conditions are fulfilled. First, the real and imaginary
parts of the coefficients $a_{n}$ are normally distributed with

$\displaystyle\langle a_{n}\rangle=0,\quad\langle
a_{n}a_{m}\rangle=0,\quad\langle a_{n}a^{*}_{m}\rangle=s_{n}\delta_{nm}.$ (16)

Second, the values of $s_{n}$ and the random vectors $\boldsymbol{\kappa}_{n}$
are chosen such that

$\displaystyle\sum_{n=1}^{N}s_{n}p_{n}(\boldsymbol{\kappa})=2\Phi_{\phi}(\boldsymbol{\kappa}),$
(17)

where $p_{n}(\boldsymbol{\kappa}_{n})$ are the probability distributions for
the random vectors $\boldsymbol{\kappa}_{n}$.

Let us split the domain of $\boldsymbol{\kappa}$ into $N$ rings,
$K_{n-1}<\kappa<K_{n}$, where

$\displaystyle
K_{n}=K_{\mathrm{min}}\exp\left[\frac{n}{N}\ln\left(\frac{K_{\mathrm{max}}}{K_{\mathrm{min}}}\right)\right],$
(18)

$n=0,\ldots,N$, and $K_{\mathrm{min}}$ and $K_{\mathrm{max}}$ are the inner
and outer bounds of the spectral domain, respectively. We define
$p_{n}(\boldsymbol{\kappa}_{n})$ such that it has zero values outside the
corresponding ring. In this case, Eq. (17) is reduced to the form

$\displaystyle
s_{n}p_{n}(\boldsymbol{\kappa}_{n})=2\Phi_{\phi}(\boldsymbol{\kappa}_{n}).$
(19)

Integrating this expression inside the ring, we get

$\displaystyle
s_{n}=2\pi\int\limits_{K_{n-1}}^{K_{n}}d\kappa\kappa\Phi_{\phi}(\kappa).$ (20)

The probability distributions for the vector $\boldsymbol{\kappa}_{n}$ are
given by
$p_{n}(\boldsymbol{\kappa}_{n})=2s_{n}^{-1}\Phi_{\phi}(\boldsymbol{\kappa}_{n})$.
We approximate it by sampling the absolute values of $\boldsymbol{\kappa}_{n}$
as

$\displaystyle\kappa_{n}=\sqrt{K_{n-1}^{2}+\xi_{n}\left(K_{n}^{2}-K_{n-1}^{2}\right)},$
(21)

where $\xi_{n}\in[0,1]$ is the uniformly distributed random variable, and by
independently sampling the polar angle of $\boldsymbol{\kappa}_{n}$ as a
uniformly distributed random variable in the domain $[0,2\pi]$.

###  II.3 Analytical and semi-analytical models

####  II.3.1 Truncated log-normal model

The log-normal distribution is a widely used model that describes light
passing through atmospheric channels [50, 51, 13]. In the domain of its
applicability, this distribution should be truncated, since the channel
efficiency cannot exceed the value of one. The normalized truncated log-normal
PDT is given by

$\displaystyle\mathcal{P}(\eta;\mu,\sigma)=\frac{1}{\mathcal{F}(1)}\frac{1}{\sqrt{2\pi}\eta\sigma}\exp\Bigl{[}-\frac{(\ln\eta+\mu)^{2}}{2\sigma^{2}}\Bigr{]},$
(22)

where $\mathcal{F}(x)$ is the cumulative distribution function of the log-
normal distribution. The parameters $\mu$ and $\sigma$ are related to the
first and second moments of the transmittance [cf. Eqs. (9) and (10)],
respectively, as

$\displaystyle\mu=\mu\left(\langle\eta\rangle,\langle\eta\rangle^{2}\right)\approx-\ln\left[\frac{\langle\eta\rangle^{2}}{\sqrt{\langle\eta^{2}\rangle}}\right],$
(23)
$\displaystyle\sigma^{2}=\sigma^{2}\left(\langle\eta\rangle,\langle\eta\rangle^{2}\right)\approx\ln\left[\frac{\langle\eta^{2}\rangle}{\langle\eta\rangle^{2}}\right].$
(24)

Thus, the truncated log-normal model is completely characterized by two
moments of the transmittance, i.e., by the information contained in the field
correlation functions (7) and (8).

####  II.3.2 Beam-wandering model

Here we briefly sketch the beam-wandering model considered in Ref. [61]. In
this model it is assumed that the beam intensity at the receiver-aperture
plane has a Gaussian shape, but $W_{0}$ is replaced with the short-term beam-
spot radius $W_{\mathrm{ST}}$,

$\displaystyle
W_{\mathrm{ST}}^{2}=W_{\mathrm{LT}}^{2}-4\sigma_{\mathrm{bw}}^{2}.$ (25)

Here $W_{\mathrm{LT}}$, defined as

$\displaystyle
W_{\mathrm{LT}}^{2}=4\int_{\mathbb{R}^{2}}d\mathbf{r}x^{2}\Gamma_{2}(\mathbf{r};z),$
(26)

is the long-term beam-spot radius and

$\displaystyle\sigma_{\mathrm{bw}}^{2}=\int_{\mathbb{R}^{4}}d\mathbf{r}_{1}d\mathbf{r}_{2}x_{1}x_{2}\Gamma_{4}(\mathbf{r}_{1},\mathbf{r}_{2};z)$
(27)

is the variance of a beam-centroid coordinate. Assuming that the beam-centroid
position is normally distributed, an analytical form for the PDT is derived in
Ref. [61], (see also Appendix C). The basic assumption behind this model is
that atmospheric turbulence results mostly in beam wandering.

####  II.3.3 Elliptic-beam model

The main idea behind the elliptic-beam model [85] is to approximate the beam
intensity by the elliptic Gaussian form

$\displaystyle
I(\mathbf{r};z_{\mathrm{ap}})=\frac{2}{\pi\sqrt{\det\mathbf{S}}}\exp\Bigl{[}-2({\mathbf{r}}{-}{\mathbf{r}}_{0})^{\rm
T}{\mathbf{S}}^{-1}({\bf r}{-}{\bf r}_{0})\Bigr{]}.$ (28)

Here

$\displaystyle\mathbf{r}_{0}=\int_{\mathbb{R}^{2}}d^{2}\mathbf{r}\,\mathbf{r}\,I(\mathbf{r};z_{\mathrm{ap}})$
(29)

and

$\displaystyle\mathbf{S}=4\int_{\mathbb{R}^{2}}d^{2}\mathbf{r}\left[(\mathbf{r}-\mathbf{r}_{0})(\mathbf{r}-\mathbf{r}_{0})^{\mathrm{T}}\right]I(\mathbf{r};z_{\mathrm{ap}})$
(30)

are the beam-centroid position and the spot-shape matrix, respectively. As it
has been discussed in Ref. [85], all parameters of the elliptic-beam model can
be approximately related to the field correlation functions (7) and (8) (see
also Appendix D).

In this paper, along with the analytical approximation considered in Ref.
[85], we will use a semianalytical technique to obtain the elliptic-beam PDT.
This method is free from assumptions about the statistical characteristics of
the elliptic-beam model. In particular, this technique involves the following
steps: (i) The values of the field intensity $I(\mathbf{r};z_{\mathrm{ap}})$
are sampled using the sparse-spectrum model for the phase-screen method, (ii)
the sampled values of $\mathbf{r}_{0}$ and $\mathbf{S}$ are obtained using
Eqs. (29) and (30), (iii) these values are successively substituted into Eqs.
(28) and (6) to obtain the sampled values of the transmittance $\eta$, and
(iv) the elliptic-beam PDT is reconstructed from the sampled values of $\eta$.
This semi-analytical approach is not considered a practical tool for obtaining
the PDT. Rather, it is used solely to assess the general applicability of the
elliptic-beam model, without relying on analytical approximations of the model
parameters.

####  II.3.4 The model based on the law of total probability

As discussed in Ref. [85], the beam-spot distortion leads to the log-normal
shape of the PDT. At the same time, beam wandering leads to the log-negative
Weibull distribution [61]. The model based on the law of total probability
[70] employs these observations and combines them into a single PDT that
accounts for two effects simultaneously. The corresponding function is given
by

$\displaystyle\mathcal{P}(\eta)=\int_{\mathbb{R}^{2}}d^{2}\boldsymbol{r}_{0}\mathcal{P}(\eta|\boldsymbol{r}_{0})\rho(\boldsymbol{r}_{0}).$
(31)

Here $\rho(\boldsymbol{r}_{0})$ is the Gaussian probability distribution for a
beam-centroid coordinate with variance (27). The function

$\displaystyle\mathcal{P}(\eta|\boldsymbol{r}_{0})=\mathcal{P}(\eta;\mu_{r_{0}},\sigma_{r_{0}})$
(32)

is the truncated log-normal PDT [cf. Eq. (22)], where
$\mu_{r_{0}}=\mu\left(\langle\eta\rangle_{r_{0}},\langle\eta^{2}\rangle_{r_{0}}\right)$
and
$\sigma_{r_{0}}=\sigma\left(\langle\eta\rangle_{r_{0}},\langle\eta^{2}\rangle_{r_{0}}\right)$
are obtained from Eqs. (23) and (24), respectively. The conditional moments
$\langle\eta\rangle_{r_{0}}$ and $\langle\eta^{2}\rangle_{r_{0}}$ are defined
similarly to (9) and (10), respectively. However, in this case the field
correlation functions are replaced by
$\Gamma_{2}^{(c)}(\mathbf{r};z)=\left\langle
I^{(c)}(\mathbf{r};z)\right\rangle$ and
$\Gamma_{4}^{(c)}(\mathbf{r}_{1},\mathbf{r}_{2};z)=\left\langle
I^{(c)}(\mathbf{r}_{1};z)I^{(c)}(\mathbf{r}_{2};z)\right\rangle$, where

$\displaystyle
I^{(c)}(\boldsymbol{r};z)=I(\boldsymbol{r}+\boldsymbol{r}_{0};z)$ (33)

is the beam intensity function with the origin of the coordinate system placed
at the position $\boldsymbol{r}_{0}$ of the beam centroid.

There are two ways to obtain the conditional moments
$\langle\eta\rangle_{r_{0}}$ and $\langle\eta^{2}\rangle_{r_{0}}$. The first
requires direct integration of the field correlation functions
$\Gamma_{2}^{(c)}(\mathbf{r};z)$ and
$\Gamma_{4}^{(c)}(\mathbf{r}_{1},\mathbf{r}_{2};z)$. This leads to an involved
integration of the analytic expressions for these functions, which in turn are
often expressed in terms of multiple integrals. However, given the numerically
simulated data, these moments can be evaluated with standard numerical
integration. The PDT obtained with such a semianalytical technique can be used
for a general analysis of the model applicability. The second possibility is
the analytical approximation based on the assumption of weak beam wandering
(see Ref. [70] and Appendix E). An analysis of the corresponding PDT will give
us information about the suitability of this approximation.

##  III Beta-distribution model

In this section we introduce an empirical model based on the Beta distribution
[98]. The choice of such an ansatz is motivated by several reasons. First,
this distribution is naturally defined in the domain $\eta\in[0,1]$. Then, it
can be parametrized by the two moments $\langle\eta\rangle$ and
$\langle\eta^{2}\rangle$. Finally, this distribution has a highly variable
shape depending on the parameter values. We will show that the Beta-
distribution model fits the results of numerical simulations better than other
analytical models for a wide range of channel parameters.

The PDT for the Beta-distribution model reads

$\displaystyle\mathcal{P}(\eta;a,b)=\frac{1}{B(a,b)}\eta^{a-1}(1-\eta)^{b-1}.$
(34)

Here $B(a,b)$ is the Beta function. The parameters $a$ and $b$ are given by

$\displaystyle
a=a\left(\langle\eta\rangle,\langle\eta^{2}\rangle\right)=\frac{\langle\eta\rangle-\langle\eta^{2}\rangle}{\langle\eta^{2}\rangle-\langle\eta\rangle^{2}}\langle\eta\rangle,$
(35) $\displaystyle
b=b\left(\langle\eta\rangle,\langle\eta^{2}\rangle\right)=a\left(\langle\eta\rangle,\langle\eta^{2}\rangle\right)\left(\frac{1}{\langle\eta\rangle}-1\right).$
(36)

These expressions parametrize the PDT in terms of $\langle\eta\rangle$ and
$\langle\eta^{2}\rangle$.

In many practical scenarios, the Beta PDT can be considered as a reasonable
alternative to the truncated log-normal distribution. This implies that one
can reformulate the model based on the law of total probability such that
$\mathcal{P}(\eta|\boldsymbol{r}_{0})$ in Eq. (31) is replaced with the Beta
PDT as

$\displaystyle\mathcal{P}(\eta|\boldsymbol{r}_{0})=\mathcal{P}(\eta;a_{r_{0}},b_{r_{0}}),$
(37)

where
$a_{r_{0}}=a\left(\langle\eta\rangle_{r_{0}},\langle\eta^{2}\rangle_{r_{0}}\right)$
and
$b_{r_{0}}=b\left(\langle\eta\rangle_{r_{0}},\langle\eta^{2}\rangle_{r_{0}}\right)$.
In this case, the conditional moments $\langle\eta\rangle_{r_{0}}$ and
$\langle\eta^{2}\rangle_{r_{0}}$ are defined in the same way as discussed in
Sec. II.3.4.

##  IV Numerical simulations

In this section we will simulate the PDTs for various horizontal channels and
compare the results with the analytical models considered in Secs. II.3 and
III. We will focus on three channels classified by the value of the Rytov
parameter $\sigma_{R}^{2}=1.23C_{n}^{2}k^{7/6}z_{\mathrm{ap}}^{11/6}$ (see
Ref. [97]). This parameter characterizes the strength of scintillation and the
overall impact of turbulence on the entire channel, taking into account both
the local turbulence strength characterized by $C_{n}^{2}$ and the channel
length $z_{\mathrm{ap}}$. We will consider three scenarios: weak impact of
turbulence with $\sigma_{\mathrm{R}}^{2}<1$, moderate impact of turbulence
with $\sigma_{\mathrm{R}}^{2}=1,\ldots,10$, and strong impact of turbulence
with $\sigma_{\mathrm{R}}^{2}>10$, as shown in Table 1. The numerically
simulated data, an interactive tool for their visualization, and the
corresponding code can be accessed in the Supplemental Material [96]. It
should be noted that our simulations do not include deterministic losses
caused by absorption and scattering. For practical applications of our
results, these losses should be additionally taken into account. In a typical
scenario, they can be considered on the order of 0.1 dB/km (see Ref. [99]).

Table 1:  Parameter values for channels with weak, moderate, and strong impact of turbulence considered in this paper.  Parameter | Weak | Moderate | Strong
---|---|---|---
Rytov parameter $\sigma_{\mathrm{R}}^{2}$ | 0.2 | 1.5 | 33.3
structure constant $C_{n}^{2}$ (m-2/3) | $5{\times}10^{-15}$ | $1.5{\times}10^{-14}$ | $6{\times}10^{-16}$
outer scale $L_{0}$ (m) | $80$ | $80$ | $80$
inner scale $\ell_{0}$ (m) | $10^{-3}$ | $10^{-3}$ | $10^{-3}$
channel length $z_{\mathrm{ap}}$ (km) | 1 | 1.6 | 50
beam-spot radius at the transmitter, $W_{0}$ (cm) | 2 | 2 | 6
wavelength $\lambda=2\pi/k$ (nm) | 809 | 809 | 808

In order to quantify the differences between the analytical or semianalytical
and numerical PDTs, we will use the Kolmogorov-Smirnov (KS) statistic (see,
e.g., Ref. [100]), given by

$\displaystyle D_{M}=\sup_{\eta}\left|F_{M}(\eta)-F(\eta)\right|.$ (38)

Here $F(\eta)=\int_{0}^{\eta}d\eta^{\prime}\mathcal{P}(\eta^{\prime})$ is the
cumulative probability distribution of transmittance corresponding to the
analytical model, $F_{M}(\eta)=M^{-1}\sum_{i=1}^{M}\theta(\eta-\eta_{i})$ is
the empirical probability distribution obtained from $M$ sampled
transmittances $\eta_{i}$, and $\theta(\eta)$ is the Heaviside step function.
If an analytical or semi-analytical model presumes samples of transmittance,
$F(\eta)$ is replaced by the empirical probability distribution obtained from
the corresponding sampled values. The cases $D_{M}=0$ and $D_{M}=1$ correspond
to complete equivalence and maximum discrepancy, respectively, between the
numerically simulated data and analytical models.

We simulate data for three channels with the technique described in Sec. II.2.
For the channels with weak and moderate impact of turbulence, we use a spatial
grid with 512 points along one axis. For the channel with strong impact of
turbulence, we use a spatial grid with 4096 points. The spatial grid steps are
0.3, 0.4, and 1 mm for the channels with weak, moderate, and strong impact of
turbulence, respectively. The number of spectral rings is $N=1024$. The inner
and outer bounds of the spectrum are $K_{\mathrm{min}}=1/15L_{0}$ and
$K_{\mathrm{max}}=2/\ell_{0}$, respectively. The number of phase screens is
chosen from the condition that the Rytov parameter for the interscreen
distance does not exceed $0.1$ (cf. Refs. [92, 101]). It is equal to 10 for
the channels with weak and moderate impact of turbulence and 30 for the
channel with strong impact of turbulence. This implies that the values of the
Rytov parameter for the interscreen distances in these cases are $3\times
10^{-3}$, $2.2\times 10^{-2}$, and $6.5\times 10^{-2}$, respectively. Small
values of this parameter enable proper simulation of atmospheric channels. The
number of samples $M=10^{5}$ for all channels.

###  IV.1 Channel with weak impact of turbulence

In this section we consider a 1 km channel with weak local turbulence, as
detailed in the second column of Table 1. These parameter values result in a
weak impact of turbulence over the entire channel. For $F_{0}=+\infty$, the
Rayleigh length is $z_{\mathrm{R}}=kW_{0}^{2}/2\approx 1.553~{}\mathrm{km}$,
which is greater than the channel length $z_{\mathrm{ap}}{=}1~{}\mathrm{km}$.
Therefore, in the absence of the atmosphere, this beam can be approximately
considered collimated. When $F_{0}=z_{\mathrm{ap}}$, the beam-spot radius at
the receiver aperture plane is minimized in the absence of the atmosphere for
a given value of $W_{0}$ and $z_{\mathrm{ap}}$. Such a beam can be considered
geometrically focused, as discussed in Ref. [97].

First, consider the collimated beam. The dependence of the KS statistic on the
aperture radius is shown in Fig. 1. Here and in the following, the aperture
radius is normalized by $W_{\mathrm{LT}}$ obtained according to Eq. (26) from
numerically simulated $\Gamma_{2}(\mathbf{r};z_{\mathrm{ap}})$. Evidently, the
beam-wandering and elliptic-beam models show the worst agreement with the
numerically simulated data. This result is counterintuitive since the beam-
spot distortion is expected to be small for weak impact of turbulence.
However, our simulations show that the contribution of this effect is still
significant compared to beam wandering. The Beta-distribution model shows the
best agreement with the numerically simulated data for the considered channel.
This model has a simple analytical form [cf. Eq. (34)] and can be easily used
for analytical considerations of various nonclassical effects for quantum
light distributed in such channels. An example of the numerically simulated
PDT and different analytical PDTs can be seen in Fig. 2. This example clearly
demonstrates that the Beta-distribution model has the best agreement with the
numerically simulated PDT in this case. The truncated log-normal and beam-
wandering models show much worse agreement.

Figure 1:  The KS statistics, quantifying the differences between the
numerically-simulated and analytical PDTs, shown for the beam-wandering model
(W), the elliptic-beam model (E), the truncated log-normal model (L), the
model based on the law of total probability with the approximation of weak
beam wandering (T${}_{\textrm{L}}$), and the Beta-distribution model (B) for
the channel with weak impact of turbulence and $F_{0}=+\infty$.  Figure 2:
Numerically simulated (N), Beta-distribution (B), beam-wandering (W), and
truncated log-normal (L) PDTs shown for the channel with weak impact of
turbulence and $F_{0}=+\infty$.

Second, we consider the geometrically focused beam. The corresponding KS
statistic as a function of aperture radius is shown in Fig. 3. The best result
for small values of the aperture radius is demonstrated by the total
probability model based on the Beta distribution. For the remaining domain,
this approximation no longer works. However, the model based on the law of
total probability without the assumption of weak beam wandering still shows
the best agreement with the numerically simulated data, as can be seen from
the KS statistic for the corresponding semianalytical model. It is important
to note that the estimation of the PDT for such semianalytical models is time
consuming. For this reason, it is performed only with $5\times 10^{3}$
samples, leading to a significant increase in the statistical error for the KS
statistics. The elliptic-beam model, while showing good agreement with the
numerically simulated data in the global minimum domain for the KS statistic,
has a small discrepancy between the semianalytical and analytical models. The
best choice among the analytical models for large values of the aperture
radius is given by the Beta-distribution model.

Figure 3:  The KS statistics, quantifying the differences between the
numerically-simulated and analytical PDTs, shown for the beam-wandering model
(W), the elliptic-beam model (E), the truncated log-normal model (L), the
model based on the law of total probability with the truncated log-normal (TL)
and Beta (T${}_{\textrm{B}}$) distributions, and the Beta-distribution model
(B) for the channel with weak impact of turbulence and
$F_{0}=z_{\mathrm{ap}}$. Solid and dashed lines correspond to analytical and
semianalytical models, respectively. Figure 4:  Numerically simulated (N) PDT,
Beta (B) PDT, elliptic-beam (E) PDT, truncated log-normal (L) PDT, and the PDT
based on the law of total probability with Beta distribution
(T${}_{\textrm{B}}$), shown for the channel with weak impact of turbulence and
$F_{0}=z_{\mathrm{ap}}$.

The numerically simulated and analytical PDTs for the geometrically focused
beam are shown in Fig. 4. It is interesting to note that the shape of the
numerically simulated PDT is similar to the PDT for the elliptic-beam model,
but for different values of the parameters. Consequently, the experimental or
numerically simulated data can still be fitted by the elliptic-beam and beam-
wandering models. However, the estimated values of the channel parameters
obtained with such fitting are biased. In general, the modes of the elliptic-
beam PDTs are less(greater) than the modes of the numerically simulated PDTs
for aperture radii less(greater) than the value corresponding to the minimum
in the KS statistic.

It is worth noting that geometrical focusing significantly improves the
channel characteristics under the considered conditions. This means that the
effects of diffraction are still significant. In particular, the mean
transmittance $\langle\eta\rangle$ is larger for $F_{0}=z_{\mathrm{ap}}$ than
for $F_{0}=+\infty$.

###  IV.2 Channel with moderate impact of turbulence

The channel considered in this section (see the third column in Table 1) is
characterized by the same beam parameters as the channel considered in Sec.
IV.1. However, its length is 1.6 km. Such a channel has been implemented in
Erlangen, Germany (see Refs. [16, 54, 56]). In contrast to the previous case,
here we consider stronger local turbulence, resulting in a moderate impact on
the light beam. Since the propagation distance is approximately equal to the
Rayleigh length, the corresponding beam is significantly divergent for
$F_{0}=+\infty$. Therefore, the geometrically focused beam,
$F_{0}=z_{\mathrm{ap}}$, is the optimal choice for given $W_{0}$ and
$z_{\mathrm{ap}}$.

The corresponding KS statistics are shown in Fig. 5. The best agreement with
the numerically simulated data is achieved by the models based on the law of
total probability. This fact becomes clear from the KS statistic for the
corresponding semianalytical model. However, these KS statistics, similar to
the previously considered case, suffer from a significant statistical error.
The Beta-distribution model shows the best agreement among the analytical
models for almost the entire range of aperture radius values. Exceptions are
the beam-wandering and elliptic-beam models, which show the best result in a
range with aperture radius close to the long-term beam-spot radius
$W_{\mathrm{LT}}$. The latter models may still be capable of fitting
numerically simulated data for a wider range of aperture radii. However, this
may lead to a biased estimation of the turbulence parameters.

Figure 5:  The KS statistics, quantifying differences between the numerically
simulated and analytical or semianalytical PDTs, shown for the beam-wandering
model (W), the elliptic-beam model (E), the truncated log-normal model (L),
the model based on the law of total probability with the truncated log-normal
(TL) and Beta (T${}_{\textrm{B}}$) distribution, and the Beta-distribution
model (B) for the channel with moderate impact of turbulence and
$F_{0}=z_{\mathrm{ap}}$. Solid and dashed lines correspond to analytical and
semi-analytical models, respectively. Figure 6:  Numerically simulated PDTs
for the channel with the moderate impact of turbulence. The cases
$F_{0}=z_{\mathrm{ap}}$ and $F_{0}=+\infty$ are given for the untracked (solid
lines) and perfectly tracked (dot-dashed lines) beams.

Examples of the numerically simulated PDTs are shown in Fig. 6 for
$F_{0}=z_{\mathrm{ap}}$ and $F_{0}=+\infty$. For the considered aperture
radius $R_{\mathrm{ap}}=1.9~{}\mathrm{cm}$, the PDT shapes do not resemble the
beam-wandering and/or elliptic-beam PDTs. We also consider a scenario with the
beam-tracking procedure (see, e.g., Ref. [102]). The idea is to reduce the
signal noise caused by beam wandering. A typical setup consists of a beacon
laser at the receiver side and a position-sensitive sensor at the transmitter
side. The latter is used for fast control of the source to align the beam
centroid with the aperture center. The transformation of the PDT due to the
partial and the complete mitigation of the beam wandering has been considered
in Ref. [70]. The PDTs for the perfectly tracked beams are shown in Fig. 6.
These distributions describe only the beam-spot distortion. One can see that
in this case, contrary to the assumptions in the model based on the law of
total probability, these PDTs do not agree with the truncated log-normal
distributions. It is also clear that the beam tracking procedure leads to an
increase in the PDT modes. This implies that the channel characteristics are
significantly improved with beam tracking.

###  IV.3 Channel with strong impact of turbulence

In this section we consider a 50 km channel with weak strength of local
turbulence (see the fourth column in Table 1). The propagation distance plays
a crucial role in such a scenario; thus the overall impact of turbulence
becomes strong. The Rayleigh length in this case is $z_{\mathrm{R}}\approx 14$
km. In the absence of atmosphere, the beam-spot radius at the aperture plane
is almost the same for $F_{0}=z_{\mathrm{ap}}$ and $F_{0}=+\infty$. Therefore,
it is reasonable to consider only the latter case.

Figure 7:  The KS statistics for the beam-wandering model (W), the elliptic-
beam model (E), the truncated log-normal model (L), the model based on the law
of total probability with the truncated log-normal (TL) and Beta
(T${}_{\textrm{B}}$) distributions, and the Beta-distribution model (B) for
the channel with strong impact of turbulence and $F_{0}=z_{\mathrm{ap}}$.
Solid and dashed lines correspond to analytical and semi-analytical models,
respectively. Figure 8:  Numerically simulated (N), truncated log-normal (L),
Beta (B), and beam-wandering (W) PDTs for the channel with strong impact of
turbulence and $F_{0}=z_{\mathrm{ap}}$, shown for different values of the
aperture radius.

The KS statistics, quantifying the differences between the numerically
simulated data and the analytical or semianalytical models for the considered
channel, are shown in Fig. 7. The general conclusions are similar to the case
of moderate impact of turbulence: The beam-wandering and elliptic-beam models
show good agreement with the numerically simulated data for the aperture
radius close to the long-term beam-spot radius $W_{\mathrm{LT}}$; outside this
range, the best result is usually given by the Beta-distribution model.
However, there is an important exception: For small apertures, the best result
is given by the truncated log-normal model and the model based on the law of
total probability with the truncated log-normal distribution. The KS statistic
for the semi-analytical models based on the law of total probability
demonstrates the best result for large apertures.

The numerical PDT and the Beta PDT for different values of the aperture radius
are shown in Fig. 8. Obviously, the PDT skewness changes with increasing
aperture; this property is inherent in all considered channels. For small
apertures, the PDT mode is left centered such that its shape is close to the
truncated log-normal distribution or the model based on the law of total
probability. This can be explained by a leading role of the beam-spot
distortion. For larger apertures, the finite size of the beam becomes
important. As a consequence, the PDT resembles the distributions for the beam-
wandering and/or elliptic-beam models, which are right centered.

##  V Statistical characteristics of beam parameters

The analytical and semianalytical models considered in the previous sections
are based on assumptions about the statistical properties of the beam
parameters. First, the beam-wandering model, the elliptic-beam model, and the
model based on the law of total probability all assume that the beam-centroid
position follows a two-dimensional Gaussian distribution. Then the model based
on the law of total probability includes the assumption that the beam-centroid
position is statistically independent of the beam shape. Finally, the
elliptic-beam model supposes a particular shape of the probability
distribution for the semiaxes of the ellipse, forming the beam shape in the
Gaussian approximation. In this section we will verify these assumptions
through numerical simulations.

###  V.1 Gaussianity of the beam-centroid position

Here we verify that the beam-centroid position follows a two-dimensional
Gaussian distribution. To accomplish this, we calculate the beam centroid
position using Eq. (29) for each sampling event and analyze the corresponding
statistical data. Note that $x_{0}$ and $y_{0}$ are statistically independent.
Therefore, it is sufficient to check whether a single coordinate follows a
Gaussian distribution.

Table 2:  Skewness and excess kurtosis for the beam-centroid coordinate $x_{0}$ in the case of channels with weak, moderate, and strong impact of turbulence, as given in Table 1. Impact of turbulence |  Wave-front radius $F_{0}$ | Skewness |  Excess kurtosis
---|---|---|---
weak | $z_{\mathrm{ap}}$ | $-0.0043$ | $-0.0038$
$+\infty$ | $0.0058$ | $0.014$
moderate | $z_{\mathrm{ap}}$ | $-0.003$ | $-0.0279$
$+\infty$ | $0.0172$ | $-0.0046$
strong | $+\infty$ | $0.0008$ | $-0.1064$
Figure 9:  Histogram obtained from the numerically simulated data,
demonstrating the distribution of the beam-centroid coordinate $x_{0}$ in the
case of a channel with strong impact of turbulence (cf. Table 1), compared
with the corresponding Gaussian distribution (dashed line).

To characterize non-Gaussianity, we will use skewness and excess kurtosis
[103]. The corresponding values for different channels are given in Table 2.
The largest discrepancy with the assumption of Gaussian distribution for the
beam-centroid position is shown by the channel with strong impact of
turbulence. However, even in this case, the skewness and excess kurtosis are
small. This implies that the numerically simulated data demonstrate good
agreement with Gaussian distributions in all cases. In Fig. 9 we compare
numerically simulated data with the corresponding Gaussian distribution for
the case where we find maximal absolute values of the skewness and excess
kurtosis. As can be seen, the deviation from the Gaussian distribution is
negligible even for such a scenario.

###  V.2 Statistical contribution of beam wandering

In this section we analyze the contribution of beam wandering into the PDT
shape. This can be characterized by correlations between the beam-deflection
distance $r_{0}=|\mathbf{r}_{0}|$ and the transmittance $\eta$. Non-zero
values of the corresponding Pearson correlation coefficient

$\displaystyle S(r_{0},\eta)=\frac{\left\langle\Delta
r_{0}\Delta\eta\right\rangle}{\sqrt{\left\langle\Delta
r_{0}^{2}\right\rangle\left\langle\Delta\eta^{2}\right\rangle}}$ (39)

indicate a contribution of beam wandering into the PDT.

Figure 10:  Pearson correlation coefficient (39) as a function of aperture
radius in units of the long beam-spot width corresponding to each scenario.
Here $W$, $M$, and $S$ correspond to the channels with weak, moderate, and
strong impact of turbulence, respectively (see Table 1). The indices $+\infty$
and $z_{\mathrm{ap}}$ indicate the cases with $F_{0}=+\infty$ and
$F_{0}=z_{\mathrm{ap}}$, respectively.

We analyze each sampling event and find the corresponding values of $\eta$ and
$r_{0}$ for it by using Eqs. (6) and (29), respectively. With this set of
sampled data, we can estimate the Pearson correlation coefficient. The results
are shown in Fig. 10 as a function of the aperture radius in units of the
corresponding long beam-spot width $R_{\mathrm{ap}}/W_{\mathrm{LT}}$.
Obviously, the contribution of beam wandering in the PDT strictly depends on
the aperture radius. The maximum anticorrelations appear for
$R_{\mathrm{ap}}\approx 0.5W_{\mathrm{LT}}$ for all channels. Also channels
with weak impact of turbulence and with $F_{0}=+\infty$ show relatively small
anticorrelations compared to other cases. Strong anticorrelations shown by the
Pearson coefficient (39) indicate the need to explicitly account for the
contribution of beam wandering such as it is done in the elliptic-beam model
and in the model based on the law of total probability.

###  V.3 Correlations between beam shape and beam deflection

The model based on the law of total probability [70] includes the assumption
that the beam shape is independent of the beam deflection. This implies that
$I^{(c)}(\mathbf{r})$ in Eq. (33) is statistically independent on
$\mathbf{r}_{0}$. This assumption will be verified in this section.

A straightforward way to do this is to study the correlations between the beam
deflection $r_{0}$ and the transmittance of the perfectly tracked beam
obtained via integration similar to Eq. (6) but for $I^{(c)}(\mathbf{r})$. For
each sampling event, we find the corresponding beam-centroid position
$\mathbf{r}_{0}$ using Eq. (29), place the aperture center there, and
calculate the corresponding transmittance $\eta$ using Eq. (6). The result is
the sample set of two random variables, $r_{0}=|\mathbf{r}_{0}|$ and $\eta$,
which we use to estimate the Pearson correlation coefficient similar to Eq.
(39). If it differs from zero, then the assumption of statistical independence
between the beam shape and the beam deflection fails.

Figure 11:  Pearson correlation coefficient, defined similarly to Eq. (39) but
with the transmittance of the perfectly tracked beam, as a function of
aperture radius in units of the long beam-spot width corresponding to each
scenario. Lines are marked as in Fig. 10.

Values of this coefficient as a function of aperture radius are shown in Fig.
11. In most cases the transmittance of the perfectly tracked beam and the
beam-deflection distance are only weakly anticorrelated at least in the
Gaussian approximation. Relatively strong anticorrelations occur only in the
case of channels with strong impact of turbulence.

Figure 12:  Tangential beam-spot width $W_{r}$. The beam centroid is deflected to the point $\mathbf{r}_{0}=(x_{0},y_{0})$. The angle $\chi$ defines the direction of the vector $\mathbf{r}_{0}$ to the $x$ axis. Table 3:  Pearson correlation coefficient between the beam-deflection distance $r_{0}$ and the radial beam width $W_{r}$ [cf. Eq. (42)], for the channels with weak, moderate, and strong impact of turbulence, as given in Table 1. Impact of turbulence | $F_{0}=+\infty$ | $F_{0}=z_{\mathrm{ap}}$
---|---|---
weak | $0.016$ | $0.039$
moderate | $0.080$ | $0.150$
strong | $0.317$ |

We also consider another way to verify statistical independence between the
beam-deflection distance and the beam shape, which does not depend on the
aperture. We characterize the beam shape using the tangential beam-spot width
$W_{r}$ as depicted in Fig. 12. For each sampling event, we find the beam-
centroid position $\mathbf{r}_{0}=(x_{0},y_{0})$ and the corresponding angle
$\chi$ that defines the direction of the vector $\mathbf{r}_{0}$ to the $x$
axis such that $\tan\chi=x_{0}/y_{0}$. In the coordinate frame
$(x_{r},y_{r})$, oriented to the beam centroid,

$\displaystyle x_{r}=x\cos\chi+y\sin\chi,$ (40) $\displaystyle
y_{r}=-x\sin\chi+y\cos\chi,$ (41)

we calculate $S_{x_{r}x_{r}}$ similarly to Eq. (30) and the corresponding beam
width as $W_{r}=\sqrt{S_{x_{r}x_{r}}}$. From the set of sampling events, we
obtain the Pearson correlation coefficient between $W_{r}$ and
$r_{0}=(x_{0}^{2}+y_{0}^{2})^{1/2}$,

$\displaystyle S(r_{0},W_{r})=\frac{\langle\Delta r_{0}\Delta
W_{r}\rangle}{\sqrt{\langle\Delta r_{0}^{2}\rangle\langle\Delta
W_{r}^{2}\rangle}},$ (42)

characterizing statistical correlations between the beam centroid and the beam
shape. The results for different channels are shown in Table 3. They
demonstrate that these statistical correlations can be significant for
channels with moderate and strong impact of turbulence. However, as seen
above, for a wide range of aperture radii this correlation does not
significantly affect the transmittance values.

###  V.4 Statistical distribution of semi-axes, forming elliptic beams

Let us consider the semiaxes of the ellipse, forming the beam shape in the
Gaussian approximation. These quantities can be defined as eigenvalues of the
matrix $\mathbf{S}$ [cf. Eq. (30)],

$\displaystyle
W_{1,2}^{2}=\frac{1}{2}\left[S_{xx}+S_{yy}\pm\sqrt{(S_{xx}-S_{yy})^{2}+4S_{xy}^{2}}\right],$

such that $W_{1}$ is assigned the eigenvalue with the sign ”$+$” or ”$-$” if
$S_{xy}\geq 0$ or $S_{xy}\leq 0$, respectively. As it is assumed in the
elliptic-beam model [85], the quantities

$\displaystyle\Theta_{1,2}=\ln W_{1,2}^{2}/W_{0}^{2}$ (43)

are distributed according to a bivariate Gaussian distribution. This
assumption is verified in this section.

For each sampling event we numerically calculate the matrix $\mathbf{S}$
according to Eq. (30) and the corresponding values of $\Theta_{1}$ and
$\Theta_{2}$. The corresponding coarse-grained scatter plot for the channel
with strong impact of turbulence and $F_{0}=z_{\mathrm{ap}}$ is shown in Fig.
13. It is compared with the covariance ellipse, given by the equation

$\displaystyle\sum\limits_{i,j=1}^{2}\big{(}\Theta_{i}-\langle\Theta_{i}\rangle\big{)}\Sigma_{ij}^{-1}\big{(}\Theta_{j}-\langle\Theta_{j}\rangle\big{)}=4,$
(44)

where $\Sigma_{ij}=\left\langle\Delta\Theta_{i}\Delta\Theta_{j}\right\rangle$
is the covariance matrix for the variables $\Theta_{1}$ and $\Theta_{2}$.
Obviously, the corresponding distribution differs significantly from the
Gaussian distribution, especially due to the pronounced minimum for
$\Theta_{1}=\Theta_{2}$. Another important observation is that the Pearson
correlation coefficient for these quantities is small; the variance ellipse
almost resembles a circle. Consequently, $\Theta_{1}$ and $\Theta_{2}$ can be
considered as uncorrelated in the second order.

Figure 13:  Coarse-grained scatter plot for $\Theta_{1}$ and $\Theta_{2}$ in the elliptic-beam model for the channel with strong impact of turbulence. The dashed line corresponds to the covariance ellipse, described by Eq. (44). Table 4:  Skewness and excess kurtosis for $\Theta_{\mathrm{(c)}}$ and $\Theta_{\mathrm{(s)}}$ in the case of channels with weak, moderate, and strong impact of turbulence, as given in Table 1. Impact of turbulence |  Wave-front radius $F_{0}$ |  Skewness for $\Theta_{(\mathrm{c})}$ and $\Theta_{(\mathrm{s})}$ |  Excess kurtosis for $\Theta_{(\mathrm{c})}$ and $\Theta_{(\mathrm{s})}$
---|---|---|---
weak | $z_{\mathrm{ap}}$ |  $2.5\times 10^{-1}$ $-6.1\times 10^{-4}$ |  $2.3\times 10^{-1}$ $-8.5\times 10^{-1}$
$+\infty$ |  $-9.9\times 10^{-2}$ $-1.0\times 10^{-2}$ |  $-6.4\times 10^{-3}$ $-1.0$
moderate | $z_{\mathrm{ap}}$ |  $2.5\times 10^{-1}$ $9.8\times 10^{-3}$ |  $1.2\times 10^{-1}$ $-8.1\times 10^{-1}$
$+\infty$ |  $1.5\times 10^{-1}$ $-3.7\times 10^{-3}$ |  $-1.1\times 10^{-1}$ $-9.9\times 10^{-1}$
strong | $+\infty$ |  $3.2\times 10^{-1}$ $-9\times 10^{-3}$ |  $2.5\times 10^{-1}$ $-7.7\times 10^{-1}$

In order to characterize non-Gaussianity, we first rotate the coordinate
system $(\Theta_{1},\Theta_{2})$ to another one,
$(\Theta_{\mathrm{(c)}},\Theta_{\mathrm{(s)}})$, as

$\displaystyle\Theta_{\mathrm{(c)}}=\frac{1}{\sqrt{2}}\left(\Theta_{1}+\Theta_{2}\right),$
(45)
$\displaystyle\Theta_{\mathrm{(s)}}=\frac{1}{\sqrt{2}}\left(\Theta_{1}-\Theta_{2}\right).$
(46)

Since the probability distribution function is symmetrical with respect to the
axis $\Theta_{\mathrm{(c)}}$, the quantities $\Theta_{\mathrm{(c)}}$ and
$\Theta_{\mathrm{(s)}}$ should be uncorrelated. We then calculate the skewness
and excess kurtosis for these quantities in order to conclude about non-
Gaussianity of the probability distribution for $\Theta_{1}$ and $\Theta_{2}$.
The results are summarized in Table 4. Significant deviations from the
Gaussian distribution indicate the further need to include this feature in the
elliptic-beam model.

##  VI Application: Transmission of quadrature-squeezed light

Quadrature squeezing is a fundamental example of nonclassical phenomena
exhibited by quantum light when measuring a field quadrature
$\hat{x}=2^{-1/2}(\hat{a}+\hat{a}^{\dagger})$, where $\hat{a}$ is the field
annihilation operator [104, 105, 106]. The quadrature variances can be
expressed as

$\displaystyle\left\langle\Delta\hat{x}^{2}\right\rangle=\frac{1}{2}+\left\langle:\Delta\hat{x}^{2}:\right\rangle,$
(47)

where the term $\frac{1}{2}$ corresponds to the quadrature variance of the
vacuum state and $\left\langle:\Delta\hat{x}^{2}:\right\rangle$ represents the
normal-ordered quadrature variance. A negative value for
$\left\langle:\Delta\hat{x}^{2}:\right\rangle$ indicates that the quadrature
variance is smaller than its vacuum-state value.

A scheme of quadrature measurement with homodyne detection for the light
passing through the turbulent atmosphere has been proposed and implemented in
Refs. [10, 11]. In this case, the local oscillator is transmitted in an
orthogonally polarized but the same spatial mode as the signal. Since the
depolarization effect of the atmosphere is negligible (see, e.g., Ref. [37]),
the transmittance of the local oscillator and the signal are strongly
correlated. This correlation provides an opportunity to control the value of
channel transmittance and postselect the events with high transmittance. Such
a scheme has been implemented in Ref. [16] and theoretically analyzed in Refs.
[55, 85].

As it has been shown in Ref. [55], the normal-ordered quadrature variance at
the receiver site,
$\left\langle:\Delta\hat{x}^{2}:\right\rangle_{\mathrm{out}}$, is related to
the normal-ordered quadrature variance at the transmitter site,
$\left\langle:\Delta\hat{x}^{2}:\right\rangle_{\mathrm{in}}$, as

$\displaystyle\left\langle:\Delta\hat{x}^{2}:\right\rangle_{\mathrm{out}}=\left\langle\eta\right\rangle\left\langle:\Delta\hat{x}^{2}:\right\rangle_{\mathrm{in}}+\left\langle\Delta
T^{2}\right\rangle\left\langle\hat{x}\right\rangle_{\mathrm{in}}^{2},$ (48)

where $T=\sqrt{\eta}$ is the transmission coefficient. It is important to note
that the constant losses on absorption and scattering in the atmosphere as
well as the losses of the optical system should be accounted for in the PDT or
in the numerically simulated data. For the squeezed vacuum state one has
$\left\langle\hat{x}\right\rangle_{\mathrm{in}}=0$ and only the first term in
Eq. (48) plays a role. The mean transmittance for the postselected signal is
given by

$\displaystyle\left\langle\eta\right\rangle=\frac{1}{\overline{\mathcal{F}}(\eta_{\mathrm{min}})}\int\limits_{\eta_{\mathrm{min}}}^{1}d\eta\,\eta\,\mathcal{P}\left(\eta\right),$
(49)

where $\eta_{\mathrm{min}}$ is the minimal value of the postselected
efficiency (the postselection threshold) and
$\overline{\mathcal{F}}(\eta_{\mathrm{min}})=\int_{\eta_{\mathrm{min}}}^{1}d\eta\,\eta\,\mathcal{P}\left(\eta\right)$
is the exceedance. The same quantity can be evaluated from the numerically
simulated data as

$\displaystyle\left\langle\eta\right\rangle=\frac{1}{N_{\mathrm{min}}}\sum\limits_{\\{\eta_{i}\geq\eta_{\mathrm{min}}\\}}\eta_{i},$
(50)

where $\eta_{i}$ are sampled values of the transmittance values and
$N_{\mathrm{min}}$ is the number of events with
$\eta_{i}\geq\eta_{\mathrm{min}}$.

We consider quadrature squeezing for the channels listed in Table 1. A typical
example of the dependence of the squeezing parameter on the postselection
threshold $\eta_{\mathrm{min}}$ is shown in Fig. 14 for the channel with
moderate impact of turbulence and $F_{0}=z_{\mathrm{ap}}$. A property common
to all the examples considered is that the elliptic-beam and beam-wandering
models show a significant discrepancy with the numerically simulated data in
the case of a small postselection threshold $\eta_{\mathrm{min}}$. This is due
to the fact that these models are not bound to the value of the first moment
$\langle\eta\rangle$. Another approach [85, 17] assumes an estimation of the
channel parameters by the method of moments based on the beam-wandering and
elliptic-beam PDTs. This leads to good agreement with experimental
(numerically-simulated) data. However, this method leads to biased values of
the estimated parameters of channels.

Figure 14:  Squeezing parameter at the receiver as a function of the
postselection threshold $\eta_{\mathrm{min}}$, evaluated with the numerically
simulated data ($\mathrm{N}$), the elliptic-beam approximation ($\mathrm{E}$),
the truncated log-normal distribution ($\mathrm{L}$), and the Beta-
distribution model ($B$), shown for the channel with moderate impact of
turbulence (cf. Table 1) and $F_{0}=z_{\mathrm{ap}}$. The initial squeezing at
the transmitter site is $-3$ dB. The constant losses are $0.38$ dB. The
receiver aperture radius is
$R_{\mathrm{ap}}=0.43W_{\mathrm{LT}}=0.019~{}\mathrm{cm}$.  Figure 15:
Squeezing parameter at the receiver as a function of the aperture radius
$R_{\mathrm{ap}}$ shown for different values of the postselection thresholds
$\eta_{\mathrm{min}}$. The channel parameters, the squeezing at the
transmitter, and the constant losses are the same as in Fig. 14.

Another interesting observation from Fig. 14 is a slight discrepancy between
the squeezing parameters obtained from the truncated log-normal distribution
and the numerically-simulated data at low postselection thresholds. This
discrepancy arises due to the influence of the tail cutoff on the first two
moments of the transmittance, leading to an error in Eqs. (23) and (24). The
semianalytical models based on the law of total probability also exhibit a
small discrepancy due to the correlations between beam shape and beam
deflection discussed in Sec. V.3. However, analytical models based on the law
of total probability with the approximation of small beam wandering do not
exhibit this discrepancy since this approximation assumes a strong constraint
on the first two moments of the transmittance.

The Beta-distribution model shows the best agreement with the numerically
simulated data. The corresponding PDT is strictly constrained to the first two
moments of the transmittance, making this model work well for small values of
the postselection threshold. However, this model also gives acceptable results
for larger values of $\eta_{\mathrm{min}}$.

The dependence of the squeezing parameters on the aperture radius is shown in
Fig. 15 for different values of the postselection threshold
$\eta_{\mathrm{min}}$. The behavior of these parameters is generally
predictable. For instance, a small value of the aperture radius may require a
large value of the postselection threshold to achieve an acceptable value of
the squeezing parameter. However, as the aperture radius approaches the long-
term beam-spot radius, the postselection becomes less significant.

##  VII Conclusions

We have presented numerical simulations of the PDT for horizontal atmospheric
quantum channels. Depending on the Rytov parameter values, we considered
scenarios with weak, strong, and moderate impact of turbulence. We performed
numerical simulations using the sparse-spectrum model for the phase-screen
method and compared our results with existing analytical models, including the
truncated log-normal distribution, beam-wandering model, elliptic-beam
approximation, and the model based on the law of total probability.
Additionally, we compared our results with the empirical model presented here
based on the Beta distribution and a corresponding modification of the model
based on the law of total probability. We also tested nontrivial statistical
assumptions on beam parameters at the receiver underlying analytical models of
the PDT.

An important conclusion drawn from this study is that the applicability of
analytical models depends strongly on the aperture radius. At the same time,
the dependence on the turbulence strength has a significantly smaller effect
on this applicability. For instance, both the beam-wandering and elliptic-beam
models demonstrate good agreement with numerically simulated data for aperture
radii of the order of or slightly exceeding the long-term beam-spot radius. If
the aperture radius is smaller or larger than this value, the PDT shape is
similar to beam-wandering and elliptic-beam PDT, but with the left- and right-
hand-side shifted distribution modes, respectively. This implies that the
elliptic-beam model may fit numerically simulated data well within a wider
range of aperture-radius values. However, such a fitting may result in biased
estimations of turbulence parameters.

The empirical Beta-distribution model presented here, as well as the
corresponding model based on the law of total probability, has shown the best
agreement with numerically simulated data over a wide range of aperture radii.
The Beta-distribution model has a simple analytical form and its parameters
depend solely on the first and second moments of the transmittance. It should
be noted that in practical scenarios, determining the range of applicability
of the beam-wandering and elliptic-beam models with precision may be
difficult. Therefore, the use of the Beta-distribution model or the
corresponding model based on the law of total probability may still be
relevant, even if the aperture radius is expected to fall within this range.
Considering this, the Beta-distribution model can be regarded as highly
applicable for theoretical research on atmospheric quantum channels.

We have applied our results to consider the transfer of quadrature squeezing
through atmospheric channels. Proper experimental design allows for monitoring
the actual transmittance value and postselecting events with high
transmittance. The resulting squeezing depends on the mean transmittance for
the postselected data. However, the beam-wandering and elliptic-beam models
usually do not accurately approximate $\langle\eta\rangle$, resulting in worse
agreement with numerically simulated data for small values of the
postselection threshold than other models. Therefore, in scenarios where the
accuracy of both $\langle\eta\rangle$ and $\langle\eta^{2}\rangle$ is crucial,
the Beta-distribution model proves to be a more reliable solution. It
demonstrates better agreement with the numerically simulated data,
particularly for small minimum values of the postselected transmittance, than
the beam-wandering and elliptic-beam models.

This work was supported by the National Research Foundation of Ukraine through
Project No. 2020.02/0111, Nonclassical and hybrid correlations of quantum
systems under realistic conditions. The authors thank D. Vasylyev and M.
Bohmann for insightful discussions.

##  Appendix A Input-output relations

In this appendix we briefly remind reader of the derivation of the input-
output relation (1) and the expression for the transmittance $\eta$ given by
Eq. (6). More details on the derivation of these expressions can be found in
Appendixes A, B, and C of the Supplemental Material in Ref. [61]. As mentioned
in Introduction, this consideration applies only to transmission protocols
with fixed spatial structure of light modes at the transmitter.

Let us consider a quasimonochromatic light mode described by the function

$\displaystyle
V_{\mathrm{in}}\left(\mathbf{r},z;t\right)=\int_{\mathbb{R}}dk\Phi(k)u_{\mathrm{in}}\left(\mathbf{r},z;k\right)e^{ik(z-ct)}.$
(51)

Here $\Phi(k)$ is a normalized spectrally narrow function,
$u_{\mathrm{in}}\left(\mathbf{r},z;k\right)$ is a solution to the paraxial
equation (2) with a fixed $\delta n(\mathbf{r},z)$, and $c$ is the speed of
light. Since $u_{\mathrm{in}}\left(\mathbf{r},z;k\right)$ varies much slower
with $z$ than $e^{ikz}$, the modes with orthonormal time-spectrum functions
can be considered as approximately orthonormal with respect to the $L_{2}$
scalar product in $\mathbb{R}^{3}$.

The light at the transmitter side is prepared in a quantum state of the quasi-
monochromatic mode $V_{\mathrm{in}}\left(\mathbf{r},z;t\right)$, which is
described by the $P$ function

$\displaystyle
P_{\mathrm{in}}(\alpha)=\frac{1}{\pi^{2}}\int_{\mathbb{C}}d^{2}\beta\operatorname{Tr}\left[\hat{\rho}\,e^{(\hat{a}^{\dagger}_{\mathrm{in}}-\alpha^{\ast})\beta}e^{-(\hat{a}_{\mathrm{in}}-\alpha)\beta^{\ast}}\right].$
(52)

Here $\hat{\rho}$ is the density operator and the annihilation operator
$\hat{a}_{\mathrm{in}}$ is expressed in terms of the positive-frequency part
of the electromagnetic-field operator
$\hat{A}_{\mathrm{in}}^{(+)}(\mathbf{r},z;t)$ as

$\displaystyle\hat{a}_{\mathrm{in}}=\int_{\mathbb{R}^{3}}d^{2}{\mathbf{r}}dz\hat{A}_{\mathrm{in}}^{(+)}(\mathbf{r},z;t)V_{\mathrm{in}}^{\ast}\left(\mathbf{r},z;t\right).$
(53)

For brevity, we consider a scalar field, i.e., polarized light, and the
dimensionless field operator.

We also consider another (unnormalized) solution to Eq. (2),

$\displaystyle
u_{\mathrm{out}}\left(\mathbf{r},z;k\right)=\int_{\mathbb{R}^{2}}d^{2}\mathbf{r}^{\prime}T(\mathbf{r},\mathbf{r}^{\prime};z;z_{\mathrm{ap}})u_{\mathrm{in}}\left(\mathbf{r}^{\prime},z;k\right).$
(54)

The integral kernel $T(\mathbf{r},\mathbf{r}^{\prime};z;z_{\mathrm{ap}})$ can
be uniquely defined from the boundary condition at the aperture plane,

$\displaystyle
u_{\mathrm{out}}\left(\mathbf{r},z_{\mathrm{ap}};k\right)=I_{\mathcal{A}}(\mathbf{r})u_{\mathrm{in}}\left(\mathbf{r},z_{\mathrm{ap}};k\right),$
(55)

where $I_{\mathcal{A}}(\mathbf{r})$ is the indicator function equal to unity
inside the aperture opening $\mathcal{A}$ and zero outside. This implies that
after passing through the aperture, the normalized quasimonochromatic mode
reads

$\displaystyle V_{\mathrm{out}}\left(\mathbf{r},z;t\right)$
$\displaystyle=\frac{1}{\sqrt{\eta}}\int_{\mathbb{R}}dk\Phi(k)u_{\mathrm{out}}\left(\mathbf{r},z;k\right)e^{ik(z-ct)},$
(56)

where the transmittance $\eta$ is defined from the normalization condition of
the function $V_{\mathrm{out}}\left(\mathbf{r},z;t\right)$ in
$\mathbb{R}^{3}$. Similarly, the expression

$\displaystyle\hat{A}_{\mathrm{out}}^{(+)}(\mathbf{r},z;t)$ (57)
$\displaystyle=\int_{\mathbb{R}^{2}}d^{2}\mathbf{r}^{\prime}T(\mathbf{r},\mathbf{r}^{\prime};z;z_{\mathrm{ap}})\hat{A}_{\mathrm{in}}^{(+)}(\mathbf{r}^{\prime},z;t)+\hat{C}(\mathbf{r},z;t)$

uniquely defines the positive-frequency field operator
$\hat{A}_{\mathrm{out}}^{(+)}(\mathbf{r},z;t)$, which for $z\geq
z_{\mathrm{ap}}$ describes the electromagnetic field after passing through the
aperture. Here $\hat{C}(\mathbf{r},z;t)$ defines noise modes in the vacuum
state. This operator is defined from the conditions of preserving the
commutation relations.

The detection system analyzes quantum states of the mode
$V_{\mathrm{out}}\left(\mathbf{r},z;t\right)$ with the $P$ function

$\displaystyle
P_{\mathrm{out}}(\alpha)=\frac{1}{\pi^{2}}\int_{\mathbb{C}}d^{2}\beta\operatorname{Tr}\left[\hat{\rho}\,e^{(\hat{a}^{\dagger}_{\mathrm{out}}-\alpha^{\ast})\beta}e^{-(\hat{a}_{\mathrm{out}}-\alpha)\beta^{\ast}}\right],$
(58)

where

$\displaystyle\hat{a}_{\mathrm{out}}=\int_{\mathbb{R}^{3}}d^{2}{\mathbf{r}}dz\hat{A}_{\mathrm{out}}^{(+)}(\mathbf{r},z;t)V_{\mathrm{out}}^{\ast}\left(\mathbf{r},z;t\right).$
(59)

Using Eq. (57) we can show that

$\displaystyle\hat{a}_{\mathrm{out}}=\sqrt{\eta}\hat{a}_{\mathrm{in}}+\sqrt{1-\eta}\hat{c},$
(60)

where $\hat{c}$ is an operator representing a noise mode in a vacuum state,
which is directly related to $\hat{C}(\mathbf{r},z;t)$ . Substituting this
expression into Eq. (58) and using Eq. (52) one gets

$\displaystyle
P_{\mathrm{out}}(\alpha)=\frac{1}{\eta}P_{\mathrm{in}}\left(\frac{\alpha}{\sqrt{\eta}}\right).$
(61)

The detection systems considered here do not distinguish between fluctuating
shapes of the modes $V_{\mathrm{out}}\left(\mathbf{r},z;t\right)$. This
implies that the final expression should be averaged with respect to $\eta$.
As a result, we arrive at the input-output relation (1).

The transmittance $\eta$ can be obtained by normalization of
$V_{\mathrm{out}}\left(\mathbf{r},z;t\right)$ [cf. Eq. (A)]. We employ the
fact that $u_{\mathrm{out}}\left(\mathbf{r},z;k\right)$ varies much slower
with $z$ than $e^{ikz}$. Together with Eq. (55) this gives

$\displaystyle\eta=\int_{\mathbb{R}}dk|\Phi(k)|^{2}\eta(k),$ (62)

where

$\displaystyle\eta(k)=\int_{\mathcal{A}}d^{2}\mathbf{r}|u_{\mathrm{in}}\left(\mathbf{r},z_{\mathrm{ap}};k\right)|^{2}.$
(63)

Since $\Phi(k)$ is a sharp normalized function, we get

$\displaystyle\eta\approx\eta(k_{0}),$ (64)

where $k_{0}$ is the maximum point of $|\Phi(k)|^{2}$. This leads to Eq. (6).
This means that the transmittance does not depend significantly on the time
spectrum $\Phi(k)$ of quasi-monochromatic modes.

##  Appendix B Verification of phase-screen method

In this appendix we present a verification technique of the sparse-spectrum
model for the phase-screen method. For this purpose we consider the phase-
structure function [37, 97]

$\displaystyle D_{\phi}$
$\displaystyle=\left\langle\left[\phi(\mathbf{r})-\phi(\mathbf{r}+\Delta\mathbf{r})\right]^{2}\right\rangle$
$\displaystyle=2\left[\left\langle\phi^{2}(\mathbf{r})\right\rangle-\left\langle\phi(\mathbf{r})\phi(\mathbf{r}+\Delta\mathbf{r})\right\rangle\right]$
(65)

where the correlation functions are given by Eq. (13). We compare these
theoretical values with the sampled values obtained by three different
methods: (i) the standard phase-screen method with a fast Fourier transform
(FFT) [89, 90, 91, 92], (ii) the method of subharmonics [107], and (iii) the
sparse-spectrum model [86, 87, 88].

The results are shown in Fig. 16. It can be clearly seen that the sparse-
spectrum model gives the best agreement between theory and numerical data. For
our simulation we use the following parameters of the atmosphere and the
simulation parameters: $\lambda=808$ nm, $C_{n}^{2}=10^{-14}$ m-2/3,
$\ell_{0}=10^{-3}$ m, $L_{0}=80$ m, $l=z_{m}-z_{m-1}=100$ m [cf. Eq. (12)], a
spatial grid composed of 1024 points for each dimension, a spatial grid step
of 2 mm, and the number of iterations $10^{4}$. For the subharmonic method, we
use six subharmonics. For the sparse-spectrum model, the number of spectral
rings is $N=1024$ and the inner and outer bounds of the spectrum are
$K_{\mathrm{min}}=1/15L_{0}$ and $K_{\mathrm{max}}=2/\ell_{0}$, respectively.

Figure 16:  Phase-structure function depending on the separation distance. The
lines correspond to the theoretical value and to the results of three
different simulations.

##  Appendix C Beam-wandering PDT

As it has been shown in Ref. [61], the beam-wandering PDT is given by the log-
negative Weibull distribution, i.e.,

$\displaystyle\mathcal{P}(\eta)=$
$\displaystyle\frac{R^{2}}{\sigma_{\mathrm{bw}}^{2}\eta\,\vartheta(2/W_{\mathrm{ST}})}\left(\ln\frac{\eta_{0}}{\eta}\right)^{\frac{2}{\vartheta(2/W_{\mathrm{ST}})}-1}$
$\displaystyle\times\exp\left[-\frac{R^{2}}{2\sigma^{2}_{\mathrm{bw}}}\left(\ln\frac{\eta_{0}}{\eta}\right)^{\frac{2}{\vartheta(2/W_{\mathrm{ST}})}}\right]$
(66)

for $\eta\in[0,\eta_{0}]$ and $\mathcal{P}(\eta)=0$ otherwise. Here

$\displaystyle\eta_{0}=1-\exp\left(-2\frac{R_{\mathrm{ap}}^{2}}{W_{\mathrm{ST}}^{2}}\right)$
(67)

is the maximum transmittance,

$\displaystyle R(\zeta)=R_{\mathrm{ap}}$ (68)
$\displaystyle\times\Bigl{[}\ln\Bigl{(}2\frac{\eta_{0}}{1-\exp\bigl{(}-R_{\mathrm{ap}}^{2}\zeta^{2}\bigr{)}}\operatorname{I}_{0}\bigl{(}R_{\mathrm{ap}}^{2}\zeta^{2}\bigr{)}\Bigr{)}\Bigr{]}^{-\frac{1}{\vartheta(\zeta)}}$

is the scale parameter, and

$\displaystyle\vartheta(\zeta)$
$\displaystyle=2R_{\mathrm{ap}}^{2}\zeta^{2}\frac{e^{-R_{\mathrm{ap}}^{2}\zeta^{2}}\operatorname{I}_{1}\bigl{(}R_{\mathrm{ap}}^{2}\zeta^{2}\bigr{)}}{1-\exp\bigl{(}-R_{\mathrm{ap}}^{2}\zeta^{2}\bigr{)}\operatorname{I}_{0}\bigl{(}R_{\mathrm{ap}}^{2}\zeta^{2}\bigr{)}}$
$\displaystyle\times\left[\ln\left(2\frac{\eta_{0}}{1-\exp\bigl{(}-R_{\mathrm{ap}}^{2}\zeta^{2}\bigr{)}\operatorname{I}_{0}\bigl{(}R_{\mathrm{ap}}^{2}\zeta^{2}\bigr{)}}\right)\right]^{-1}$
(69)

is the shape parameter. In these equations, $R_{\mathrm{ap}}$ is the aperture
radius, $W_{\mathrm{ST}}^{2}$ is the short term [cf. Eq. (25)],
$\sigma_{\mathrm{bw}}^{2}$ is the beam-wandering variance [cf. Eq. (27)],
$\zeta=2/W_{\mathrm{ST}}$, and $\operatorname{I}_{n}(x)$ is the modified
Bessel function.

##  Appendix D Analytical approximation for the elliptic-beam model

In order to formulate the analytical approximation used in the elliptic-beam
model [85], we first need an expression for the transmittance of the
elliptical beam (28) through the aperture of the radius $R_{\mathrm{ap}}$, as
a function of the vector
$\mathbf{v}=\begin{pmatrix}x_{0}&y_{0}&\Theta_{1}&\Theta_{2}\end{pmatrix}$ and
the angle $\phi$ between the axis $x$ and the ellipse semiaxis $W_{1}$. Here
$x_{0}$ and $y_{0}$ are the beam-centroid coordinates and $\Theta_{1,2}$ are
related to the length of the ellipse semiaxes, $W_{1,2}$, as given by Eq.
(43). This transmittance reads

$\displaystyle\eta=\widetilde{\eta}_{0}\exp\left\\{-\left[\frac{r_{0}}{R\left(\frac{2}{W_{\rm
eff}(\phi-\chi)}\right)}\right]^{\vartheta\left(\frac{2}{W_{\rm
eff}(\phi-\chi)}\right)}\right\\},$ (70)

where $r_{0}=|\mathbf{r}_{0}|$ is the beam-deflection distance; the scale
parameter $R(\zeta)$ and the shape parameter $\vartheta(\zeta)$ are given by
Eqs. (68) and (C), respectively;

$\displaystyle\widetilde{\eta}_{0}$
$\displaystyle=1{-}\operatorname{I}_{0}\left(R_{\mathrm{ap}}^{2}\frac{W_{1}^{2}-W_{2}^{2}}{W_{1}^{2}W_{2}^{2}}\right)e^{-R_{\mathrm{ap}}^{2}\frac{W_{1}^{2}+W_{2}^{2}}{W_{1}^{2}W_{2}^{2}}}$
$\displaystyle-2\Bigl{[}1{-}e^{-\frac{R_{\mathrm{ap}}^{2}}{2}\Bigl{(}\frac{1}{W_{1}}{-}\frac{1}{W_{2}}\Bigr{)}^{2}}\Bigr{]}$
(71)
$\displaystyle\times\exp\left[-\left\\{\frac{R_{\mathrm{ap}}\frac{(W_{1}+W_{2})^{2}}{|W_{1}^{2}-W_{2}^{2}|}}{R\left(\frac{1}{W_{1}}-\frac{1}{W_{2}}\right)}\right\\}^{\vartheta\left(\frac{1}{W_{1}}-\frac{1}{W_{2}}\right)}\right],$

is the maximum transmittance for the elliptic beam;

$\displaystyle W_{\rm
eff}^{2}(\phi-\chi)=4R_{\mathrm{ap}}^{2}\Bigl{[}\mathcal{W}\Bigl{(}\frac{4R_{\mathrm{ap}}^{2}}{W_{1}W_{2}}e^{2R_{\mathrm{ap}}^{2}\left(\frac{1}{W_{1}^{2}}+\frac{1}{W_{2}^{2}}\right)}\Bigr{.}\Bigr{.}$
(72)
$\displaystyle\qquad\times\Bigl{.}\Bigl{.}e^{R_{\mathrm{ap}}^{2}\left(\frac{1}{W_{1}^{2}}-\frac{1}{W_{2}^{2}}\right)\cos(2\phi-2\chi)}\Bigr{)}\Bigr{]}^{-1};$

$\chi$ is the angle that defines the direction of the vector $\mathbf{r}_{0}$
to the $x$ axis (cf. Fig. 12); $\operatorname{I}_{n}(x)$ is the modified
Bessel function; and $\mathcal{W}$ is the Lambert function.

The elliptic-beam PDT is given by

$\displaystyle\mathcal{P}\left(\eta\right){=}\frac{2}{\pi}\int_{\mathbb{R}^{4}}d^{4}\mathbf{v}\int\limits_{0}^{\pi/2}d\phi\,\rho_{G}(\mathbf{v}|\boldsymbol{\mu},\Sigma)\delta\left[\eta{-}\eta\left(\mathbf{v},\phi\right)\right].$
(73)

Here $\eta\left(\mathbf{v},\phi\right)$ is the transmittance defined by Eq.
(70) as a function of random parameters and
$\rho_{G}(\mathbf{v}|\boldsymbol{\mu},\Sigma)$ is the Gaussian probability
density function of the vector $\mathbf{v}$ with the mean $\boldsymbol{\mu}$
and the covariance matrix $\Sigma$. This matrix is clearly divided into two
diagonal blocks since $\langle\Delta\mathbf{r}_{0}\Delta\Theta_{1,2}\rangle=0$
due to symmetry. The first block, related to the beam-centroid coordinates
$x_{0}$ and $y_{0}$, reads $\sigma_{\mathrm{bw}}^{2}\mathbb{I}_{2}$, where
$\mathbb{I}_{2}$ is the $2\times 2$ unity matrix. The second block consists of
the elements

$\displaystyle\langle\Delta\Theta_{i}\Delta\Theta_{j}\rangle=\ln\left[1+\frac{\langle\Delta
W_{i}^{2}\Delta W_{j}^{2}\rangle}{\langle W_{i}^{2}\rangle\langle
W_{j}^{2}\rangle}\right].$ (74)

Two first components of the vector $\boldsymbol{\mu}$ are zero in the case
when $\langle\mathbf{r}_{0}\rangle=0$. The two other components are given by

$\displaystyle\langle\Theta_{i}\rangle=\ln\left[\frac{\langle
W_{i}^{2}\rangle}{W_{0}^{2}}\left(1+\frac{\langle(\Delta
W_{i}^{2})^{2}\rangle}{\langle W_{i}^{2}\rangle^{2}}\right)^{-1/2}\right].$
(75)

The moments of $W_{i}^{2}$ in Eqs. (74) and (75) are expressed in terms of
field correlation functions
$\Gamma_{2}\\!\left(\mathbf{r};z_{\mathrm{ap}}\right)$ and
$\Gamma_{4}\\!\left(\mathbf{r}_{1},\mathbf{r}_{2};z_{\mathrm{ap}}\right)$ as

$\displaystyle\langle$ $\displaystyle
W_{i}^{2}\rangle{=}4\left[\int_{\mathbb{R}^{2}}d^{2}\mathbf{r}\,x^{2}\Gamma_{2}\\!\left(\mathbf{r};z_{\mathrm{ap}}\right){-}\langle
x_{0}^{2}\rangle\right],$ (76) $\displaystyle\langle
W_{i}^{2}W_{j}^{2}\rangle=8\Big{[}{-}8\,\delta_{ij}\langle
x_{0}^{2}\rangle^{2}{-}\langle x_{0}^{2}\rangle\langle W_{i}^{2}\rangle$ (77)
$\displaystyle+\int_{\mathbb{R}^{4}}d^{4}\mathbf{r}\,\left[x_{1}^{2}x_{2}^{2}\left(4\delta_{ij}{-}1\right)-x_{1}^{2}y_{2}^{2}\left(4\delta_{ij}{-}3\right)\right]$
$\displaystyle\times\,\Gamma_{4}\\!\left(\mathbf{r}_{1},\mathbf{r}_{2};z_{\mathrm{ap}}\right)\Big{]}.$

These moments can be calculated analytically, for instance, using the phase
approximation of the Huygens-Kirchhoff method, as demonstrated in Ref. [85].
However, in this paper we calculated these moments from the numerically
simulated data.

##  Appendix E Approximation of weak beam-wandering for the model based on the
law of total probability

In both variants of the model based on the law of total probability, i.e.,
those based on the truncated log-normal distribution (cf. Sec. II.3.4) and
those based on the Beta distribution (cf. Sec. III), one needs to evaluate the
conditional moments of the transmittance, $\langle\eta\rangle_{r_{0}}$ and
$\langle\eta^{2}\rangle_{r_{0}}$. In the approximation of weak beam wandering
[70] these quantities are evaluated as

$\displaystyle\langle\eta\rangle_{r_{0}}=\eta_{0}\exp\left[-\left(\frac{r_{0}}{R(2/W_{\mathrm{ST}})}\right)^{\vartheta(2/W_{\mathrm{ST}})}\right],$
(78)
$\displaystyle\langle\eta^{2}\rangle_{r_{0}}=h_{0}^{2}\exp\left[-2\left(\frac{r_{0}}{R(2/W_{\mathrm{ST}})}\right)^{\vartheta(2/W_{\mathrm{ST}})}\right],$
(79)

where the scale parameter $R(\zeta)$ and the shape parameter
$\vartheta(\zeta)$ are given by Eqs. (68) and (C), respectively. The
efficiency $\eta_{0}$ can be obtained from Eq. (67) if the field correlation
function $\Gamma_{2}(\mathbf{r};z_{\mathrm{ap}})$ is Gaussian. Otherwise, it
is obtained as

$\displaystyle\eta_{0}=\frac{\langle\eta\rangle}{\displaystyle{\int\limits_{0}^{\infty}du\,u\,e^{-\frac{u^{2}}{2}}e^{-\left(\frac{\sigma_{\mathrm{bw}}}{R(2/W_{\mathrm{ST}})}u\right)^{\vartheta(2/W_{\mathrm{ST}})}}}}.$
(80)

The parameter $h_{0}$ is evaluated as

$\displaystyle
h_{0}^{2}=\frac{\langle\eta^{2}\rangle}{\displaystyle{\int\limits_{0}^{\infty}du\,u\,e^{-\frac{u^{2}}{2}}e^{-2\left(\frac{\sigma_{\mathrm{bw}}}{R(2/W_{\mathrm{ST}})}u\right)^{\vartheta(2/W_{\mathrm{ST}})}}}}.$
(81)

In addition, this approximation ensures that two important statistical moments
of transmittance, namely, $\langle\eta\rangle$ and $\langle\eta^{2}\rangle$,
retain their original definitions given by Eqs. (9) and (10), respectively.

## References

  * Gisin  _et al._ [2002] N. Gisin, G. Ribordy, W. Tittel, and H. Zbinden, Quantum cryptography, Rev. Mod. Phys. 74, 145 (2002).
  * Xu  _et al._ [2020] F. Xu, X. Ma, Q. Zhang, H.-K. Lo, and J.-W. Pan, Secure quantum key distribution with realistic devices, Rev. Mod. Phys. 92, 025002 (2020).
  * Pirandola  _et al._ [2020] S. Pirandola, U. L. Andersen, L. Banchi, M. Berta, B. D., R. Colbeck, D. Englund, T. Gehring, C. Lupo, C. Ottaviani, J. L. Pereira, M. Razavi, J. Shamsul Shaari, M. Tomamichel, V. C. Usenko, G. Vallone, P. Villoresi, and P. Wallden, Advances in quantum cryptography, Adv. Opt. Photon. 12, 1012 (2020).
  * Renner and Wolf [2023] R. Renner and R. Wolf, Quantum advantage in cryptography, AIAA Journal 61, 1895 (2023).
  * Gottesman and Chuang [2001] D. Gottesman and I. Chuang, Quantum digital signatures (2001), arXiv:arXiv:quant-ph/0105032 [quant-ph] .
  * Bennett  _et al._ [1993] C. H. Bennett, G. Brassard, C. Crépeau, R. Jozsa, A. Peres, and W. K. Wootters, Teleporting an unknown quantum state via dual classical and einstein-podolsky-rosen channels, Phys. Rev. Lett. 70, 1895 (1993).
  * Braunstein and Kimble [1998] S. L. Braunstein and H. J. Kimble, Teleportation of continuous quantum variables, Phys. Rev. Lett. 80, 869 (1998).
  * Żukowski  _et al._ [1993] M. Żukowski, A. Zeilinger, M. A. Horne, and A. K. Ekert, “event-ready-detectors” bell experiment via entanglement swapping, Phys. Rev. Lett. 71, 4287 (1993).
  * Ursin  _et al._ [2007] R. Ursin _et al._ , Entanglement-based quantum communication over 144 km, Nat. Phys. 3, 481 (2007).
  * Elser  _et al._ [2009] D. Elser, T. Bartley, B. Heim, C. Wittmann, D. Sych, and G. Leuchs, Feasibility of free space quantum key distribution with coherent polarization states, New J. Phys. 11, 045014 (2009).
  * Heim  _et al._ [2010] B. Heim, D. Elser, T. Bartley, M. Sabuncu, C. Wittmann, D. Sych, C. Marquardt, and G. Leuchs, Atmospheric channel characteristics for quantum communication with continuous polarization variables, Appl. Phys. B 98, 635 (2010).
  * Fedrizzi  _et al._ [2009] A. Fedrizzi, R. Ursin, T. Herbst, M. Nespoli, R. Prevedel, T. Scheidl, F. Tiefenbacher, T. Jennewein, and A. Zeilinger, High-fidelity transmission of entanglement over a high-loss free-space channel, Nat. Phys. 5, 389 (2009).
  * Capraro  _et al._ [2012] I. Capraro, A. Tomaello, A. Dall’Arche, F. Gerlin, R. Ursin, G. Vallone, and P. Villoresi, Impact of turbulence in long range quantum and classical communications, Phys. Rev. Lett. 109, 200502 (2012).
  * Yin  _et al._ [2012] J. Yin _et al._ , Quantum teleportation and entanglement distribution over 100-kilometre free-space channels, Nature (London) 488, 185 (2012).
  * Ma  _et al._ [2012] X.-S. Ma _et al._ , Quantum teleportation over 143 kilometres using active feed-forward, Nature (London) 489, 269 (2012).
  * Peuntinger  _et al._ [2014] C. Peuntinger, B. Heim, C. R. Müller, C. Gabriel, C. Marquardt, and G. Leuchs, Distribution of squeezed states through an atmospheric channel, Phys. Rev. Lett. 113, 060502 (2014).
  * Vasylyev  _et al._ [2017] D. Vasylyev, A. A. Semenov, W. Vogel, K. Günthner, A. Thurn, O. Bayraktar, and C. Marquardt, Free-space quantum links under diverse weather conditions, Phys. Rev. A 96, 043856 (2017).
  * Jin  _et al._ [2019] J. Jin, J.-P. Bourgoin, R. Tannous, S. Agne, C. J. Pugh, K. B. Kuntz, B. L. Higgins, and T. Jennewein, Genuine time-bin-encoded quantum key distribution over a turbulent depolarizing free-space channel, Opt. Express 27, 37214 (2019).
  * Bulla  _et al._ [2023] L. Bulla, M. Pivoluska, K. Hjorth, O. Kohout, J. Lang, S. Ecker, S. P. Neumann, J. Bittermann, R. Kindler, M. Huber, M. Bohmann, and R. Ursin, Nonlocal temporal interferometry for highly resilient free-space quantum communication, Phys. Rev. X 13, 021001 (2023).
  * Nauerth  _et al._ [2013] S. Nauerth, F. Moll, M. Rau, C. Fuchs, J. Horwath, S. Frick, and H. Weinfurter, Air-to-ground quantum communication, Nat. Photonics 7, 382 (2013), letter.
  * Wang  _et al._ [2013] J.-Y. Wang _et al._ , Direct and full-scale experimental verifications towards ground-satellite quantum key distribution, Nat. Photonics 7, 387 (2013).
  * Bourgoin  _et al._ [2015] J.-P. Bourgoin, B. L. Higgins, N. Gigov, C. Holloway, C. J. Pugh, S. Kaiser, M. Cranmer, and T. Jennewein, Free-space quantum key distribution to a moving receiver, Opt. Express 23, 33437 (2015).
  * Vallone  _et al._ [2015] G. Vallone, D. Bacco, D. Dequal, S. Gaiarin, V. Luceri, G. Bianco, and P. Villoresi, Experimental satellite quantum communications, Phys. Rev. Lett. 115, 040502 (2015).
  * Dequal  _et al._ [2016] D. Dequal, G. Vallone, D. Bacco, S. Gaiarin, V. Luceri, G. Bianco, and P. Villoresi, Experimental single-photon exchange along a space link of 7000 km, Phys. Rev. A 93, 010301(R) (2016).
  * Vallone  _et al._ [2016] G. Vallone, D. Dequal, M. Tomasin, F. Vedovato, M. Schiavon, V. Luceri, G. Bianco, and P. Villoresi, Interference at the single photon level along satellite-ground channels, Phys. Rev. Lett. 116, 253601 (2016).
  * Carrasco-Casado  _et al._ [2016] A. Carrasco-Casado, H. Kunimori, H. Takenaka, T. Kubo-Oka, M. Akioka, T. Fuse, Y. Koyama, D. Kolev, Y. Munemasa, and M. Toyoshima, LEO-to-ground polarization measurements aiming for space qkd using small optical transponder (SOTA), Opt. Express 24, 12254 (2016).
  * Takenaka  _et al._ [2017] H. Takenaka, A. Carrasco-Casado, M. Fujiwara, M. Kitamura, M. Sasaki, and M. Toyoshima, Satellite-to-ground quantum-limited communication using a 50-kg-class microsatellite, Nat. Photonics 11, 502 (2017).
  * Liao  _et al._ [2017] S.-K. Liao _et al._ , Satellite-to-ground quantum key distribution, Nature (London) 549, 43 (2017).
  * Yin  _et al._ [2017a] J. Yin _et al._ , Satellite-based entanglement distribution over 1200 kilometers, Science 356, 1140 (2017a).
  * Günthner  _et al._ [2017] K. Günthner _et al._ , Quantum-limited measurements of optical signals from a geostationary satellite, Optica 4, 611 (2017).
  * Ren  _et al._ [2017] J.-G. Ren _et al._ , Ground-to-satellite quantum teleportation, Nature (London) 549, 70 (2017).
  * Yin  _et al._ [2017b] J. Yin, Y. Cao, Y.-H. Li, J.-G. Ren, S.-K. Liao, L. Zhang, _et al._ , Satellite-to-ground entanglement-based quantum key distribution, Phys. Rev. Lett. 119, 200501 (2017b).
  * Liao  _et al._ [2018] S.-K. Liao, W.-Q. Cai, J. Handsteiner, B. Liu, J. Yin, L. Zhang, _et al._ , Satellite-relayed intercontinental quantum network, Phys. Rev. Lett. 120, 030501 (2018).
  * Vasylyev  _et al._ [2019] D. Vasylyev, W. Vogel, and F. Moll, Satellite-mediated quantum atmospheric links, Phys. Rev. A 99, 053830 (2019).
  * Ecker  _et al._ [2021] S. Ecker, B. Liu, J. Handsteiner, M. Fink, D. Rauch, F. Steinlechner, T. Scheidl, A. Zeilinger, and R. Ursin, Strategies for achieving high key rates in satellite-based QKD, npj Quantum Information 7, 5 (2021).
  * Ecker  _et al._ [2023] S. Ecker, J. Pseiner, J. Piris, and M. Bohmann, Advances in entanglement-based QKD for space applications, in  _International Conference on Space Optics — ICSO 2022_, Vol. 12777, edited by K. Minoglou, N. Karafolas, and B. Cugny, International Society for Optics and Photonics (SPIE, 2023) p. 1277727.
  * Tatarskii [1971] V. Tatarskii,  _The Effect of the Turbulent Atmosphere on Wave Propagation_ (Israel Program for Scientific Translations, Jerusalem, 1971).
  * Fante [1975] R. Fante, Electromagnetic beam propagation in turbulent media, Proc. IEEE 63, 1669 (1975).
  * Fante [1980] R. Fante, Electromagnetic beam propagation in turbulent media: An update, Proc. IEEE 68, 1424 (1980).
  * D’Ambrosio  _et al._ [2012] V. D’Ambrosio, E. Nagali, S. P. Walborn, L. Aolita, S. Slussarenko, L. Marrucci, and F. Sciarrino, Complete experimental toolbox for alignment-free quantum communication, Nature Commun. 3, 961 (2012).
  * Vallone  _et al._ [2014] G. Vallone, V. D’Ambrosio, A. Sponselli, S. Slussarenko, L. Marrucci, F. Sciarrino, and P. Villoresi, Free-space quantum key distribution by rotation-invariant twisted photons, Phys. Rev. Lett. 113, 060503 (2014).
  * Krenn  _et al._ [2016] M. Krenn, J. Handsteiner, M. Fink, R. Fickler, R. Ursin, M. Malik, and A. Zeilinger, Twisted light transmission over 143 km, PNAS 113, 13648 (2016).
  * Paterson [2005] C. Paterson, Atmospheric turbulence and orbital angular momentum of single photons for optical communication, Phys. Rev. Lett. 94, 153901 (2005).
  * Pors  _et al._ [2011] B.-J. Pors, C. H. Monken, E. R. Eliel, and J. P. Woerdman, Transport of orbital-angular-momentum entanglement through a turbulent atmosphere, Opt. Express 19, 6671 (2011).
  * Roux [2011] F. S. Roux, Infinitesimal-propagation equation for decoherence of an orbital-angular-momentum-entangled biphoton state in atmospheric turbulence, Phys. Rev. A 83, 053822 (2011).
  * Roux  _et al._ [2015] F. S. Roux, T. Wellens, and V. N. Shatokhin, Entanglement evolution of twisted photons in strong atmospheric turbulence, Phys. Rev. A 92, 012326 (2015).
  * Leonhard  _et al._ [2015] N. D. Leonhard, V. N. Shatokhin, and A. Buchleitner, Universal entanglement decay of photonic-orbital-angular-momentum qubit states in atmospheric turbulence, Phys. Rev. A 91, 012345 (2015).
  * Lavery [2018] M. P. J. Lavery, Vortex instability in turbulent free-space propagation, New J. Phys. 20, 043023 (2018).
  * Sorelli  _et al._ [2019] G. Sorelli, N. Leonhard, V. N. Shatokhin, C. Reinlein, and A. Buchleitner, Entanglement protection of high-dimensional states by adaptive optics, New J. Phys. 21, 023003 (2019).
  * Peřina  _et al._ [1973] J. Peřina, V. Peřinová, M. C. Teich, and P. Diament, Two descriptions for the photocounting detection of radiation passed through a random medium: A comparison for the turbulent atmosphere, Phys. Rev. A 7, 1732 (1973).
  * Milonni  _et al._ [2004] P. W. Milonni, J. H. Carter, C. G. Peterson, and R. J. Hughes, Effects of propagation through atmospheric turbulence on photon statistics, J. Opt. B: Quantum Semiclass. Opt. 6, S742 (2004).
  * Semenov and Vogel [2010] A. A. Semenov and W. Vogel, Entanglement transfer through the turbulent atmosphere, Phys. Rev. A 81, 023835 (2010).
  * Gumberidze  _et al._ [2016] M. O. Gumberidze, A. A. Semenov, D. Vasylyev, and W. Vogel, Bell nonlocality in the turbulent atmosphere, Phys. Rev. A 94, 053801 (2016).
  * Usenko  _et al._ [2012] V. C. Usenko, B. Heim, C. Peuntinger, C. Wittmann, C. Marquardt, G. Leuchs, and R. Filip, Entanglement of Gaussian states and the applicability to quantum key distribution over fading channels, New J. Phys. 14, 093048 (2012).
  * Semenov  _et al._ [2012] A. A. Semenov, F. Töppel, D. Y. Vasylyev, H. V. Gomonay, and W. Vogel, Homodyne detection for atmosphere channels, Phys. Rev. A 85, 013826 (2012).
  * Heim  _et al._ [2014] B. Heim, C. Peuntinger, N. Killoran, I. Khan, C. Wittmann, C. Marquardt, and G. Leuchs, Atmospheric continuous-variable quantum communication, New J. Phys. 16, 113018 (2014).
  * Shapiro [2011] J. H. Shapiro, Scintillation has minimal impact on far-field bennett-brassard 1984 protocol quantum key distribution, Phys. Rev. A 84, 032340 (2011).
  * Bonato  _et al._ [2009] C. Bonato, A. Tomaello, V. D. Deppo, G. Naletto, and P. Villoresi, Feasibility of satellite quantum key distribution, New J. Phys. 11, 045017 (2009).
  * Lanning  _et al._ [2021] R. N. Lanning, M. A. Harris, D. W. Oesch, M. D. Oliker, and M. T. Gruneisen, Quantum communication over atmospheric channels: A framework for optimizing wavelength and filtering, Phys. Rev. Applied 16, 044027 (2021).
  * Semenov and Vogel [2009] A. A. Semenov and W. Vogel, Quantum light in the turbulent atmosphere, Phys. Rev. A 80, 021802(R) (2009).
  * Vasylyev  _et al._ [2012] D. Y. Vasylyev, A. A. Semenov, and W. Vogel, Toward global quantum communication: Beam wandering preserves nonclassicality, Phys. Rev. Lett. 108, 220501 (2012).
  * Glauber [1963] R. J. Glauber, Coherent and incoherent states of the radiation field, Phys. Rev. 131, 2766 (1963).
  * Sudarshan [1963] E. C. G. Sudarshan, Equivalence of semiclassical and quantum mechanical descriptions of statistical light beams, Phys. Rev. Lett. 10, 277 (1963).
  * Semenov  _et al._ [2018] A. A. Semenov, M. Bohmann, D. Vasylyev, and W. Vogel, Nonclassicality and Bell nonlocality in atmospheric links, in  _Quantum Communications and Quantum Imaging XVI_, Vol. 10771, edited by R. E. Meyers, Y. Shih, and K. S. Deacon, International Society for Optics and Photonics (SPIE, 2018) pp. 155 – 170\.
  * Bohmann  _et al._ [2016] M. Bohmann, A. A. Semenov, J. Sperling, and W. Vogel, Gaussian entanglement in the turbulent atmosphere, Phys. Rev. A 94, 010302(R) (2016).
  * Hosseinidehaj and Malaney [2015a] N. Hosseinidehaj and R. Malaney, Gaussian entanglement distribution via satellite, Phys. Rev. A 91, 022304 (2015a).
  * Bohmann  _et al._ [2017] M. Bohmann, J. Sperling, A. A. Semenov, and W. Vogel, Higher-order nonclassical effects in fluctuating-loss channels, Phys. Rev. A 95, 012324 (2017).
  * Hosseinidehaj and Malaney [2015b] N. Hosseinidehaj and R. Malaney, Entanglement generation via non-gaussian transfer over atmospheric fading channels, Phys. Rev. A 92, 062336 (2015b).
  * Wang  _et al._ [2018a] W. Wang, F. Xu, and H.-K. Lo, Prefixed-threshold real-time selection method in free-space quantum key distribution, Phys. Rev. A 97, 032337 (2018a).
  * Vasylyev  _et al._ [2018] D. Vasylyev, W. Vogel, and A. A. Semenov, Theory of atmospheric quantum channels based on the law of total probability, Phys. Rev. A 97, 063852 (2018).
  * Guo  _et al._ [2017] Y. Guo, C. Xie, Q. Liao, W. Zhao, G. Zeng, and D. Huang, Entanglement-distillation attack on continuous-variable quantum key distribution in a turbulent atmospheric channel, Phys. Rev. A 96, 022320 (2017).
  * Papanastasiou  _et al._ [2018] P. Papanastasiou, C. Weedbrook, and S. Pirandola, Continuous-variable quantum key distribution in uniform fast-fading channels, Phys. Rev. A 97, 032311 (2018).
  * Derkach  _et al._ [2020] I. Derkach, V. C. Usenko, and R. Filip, Squeezing-enhanced quantum key distribution over atmospheric channels, New J. Phys. 22, 053006 (2020).
  * Hosseinidehaj  _et al._ [2019] N. Hosseinidehaj, Z. Babar, R. Malaney, S. X. Ng, and L. Hanzo, Satellite-based continuous-variable quantum communications: State-of-the-art and a predictive outlook, IEEE Commun. Surv. Tutor. 21, 881 (2019).
  * Wang  _et al._ [2018b] S. Wang, P. Huang, T. Wang, and G. Zeng, Atmospheric effects on continuous-variable quantum key distribution, New J. Phys. 20, 083037 (2018b).
  * Chai  _et al._ [2019] G. Chai, Z. Cao, W. Liu, S. Wang, P. Huang, and G. Zeng, Parameter estimation of atmospheric continuous-variable quantum key distribution, Phys. Rev. A 99, 032326 (2019).
  * Derkach and Usenko [2021] I. Derkach and V. C. Usenko, Applicability of squeezed- and coherent-state continuous-variable quantum key distribution over satellite links, Entropy 23, 10.3390/e23010055 (2021).
  * Pirandola [2021a] S. Pirandola, Composable security for continuous variable quantum key distribution: Trust levels and practical key rates in wired and wireless networks, Phys. Rev. Research 3, 043014 (2021a).
  * Hosseinidehaj  _et al._ [2021] N. Hosseinidehaj, N. Walk, and T. C. Ralph, Composable finite-size effects in free-space continuous-variable quantum-key-distribution systems, Phys. Rev. A 103, 012605 (2021).
  * Pirandola [2021b] S. Pirandola, Satellite quantum communications: Fundamental bounds and practical security, Phys. Rev. Research 3, 023130 (2021b).
  * Pirandola [2021c] S. Pirandola, Limits and security of free-space quantum communications, Phys. Rev. Research 3, 013279 (2021c).
  * Hofmann  _et al._ [2019] K. Hofmann, A. A. Semenov, W. Vogel, and M. Bohmann, Quantum teleportation through atmospheric channels, Phys. Scr. 94, 125104 (2019).
  * Sheng-Li  _et al._ [2017] Sheng-Li, C.-H. Jin, J.-H. Shi, J.-S. Guo, X.-B. Zou, and G.-C. Guo, Continuous variable quantum teleportation in beam-wandering modeled atmosphere channel, Chinese Phys. Lett. 34, 040302 (2017).
  * Villaseñor  _et al._ [2021] E. Villaseñor, M. He, Z. Wang, R. Malaney, and M. Z. Win, Enhanced uplink quantum communication with satellites via downlink channels, IEEE Trans. Quant. Eng. 2, 1 (2021).
  * Vasylyev  _et al._ [2016] D. Vasylyev, A. A. Semenov, and W. Vogel, Atmospheric quantum channels with weak and strong turbulence, Phys. Rev. Lett. 117, 090501 (2016).
  * Charnotskii [2013a] M. Charnotskii, Sparse spectrum model for a turbulent phase, J. Opt. Soc. Am. A 30, 479 (2013a).
  * Charnotskii [2013b] M. Charnotskii, Statistics of the sparse spectrum turbulent phase, J. Opt. Soc. Am. A 30, 2455 (2013b).
  * Charnotskii [2020] M. Charnotskii, Comparison of four techniques for turbulent phase screens simulation, J. Opt. Soc. Am. A 37, 738 (2020).
  * Fleck  _et al._ [1976] J. A. Fleck, J. R. Morris, and M. D. Feit, Time-dependent propagation of high energy laser beams through the atmosphere, Appl. Phys. 10, 129 (1976).
  * Frehlich [2000] R. Frehlich, Simulation of laser propagation in a turbulent atmosphere, Appl. Opt. 39, 393 (2000).
  * Lukin and Fortes [2002] V. P. Lukin and B. V. Fortes,  _Adaptive Beaming and Imaging in the Turbulent Atmosphere_ (SPIE, 2002).
  * Schmidt [2010] J. Schmidt,  _Numerical simulations of optical wave propagation with examples in MATLAB_ (SPIE, Bellingham, 2010).
  * Lushnikov and Vladimirova [2018] P. M. Lushnikov and N. Vladimirova, Toward defeating diffraction and randomness for laser beam propagation in turbulent atmosphere, JETP Letters 108, 571 (2018).
  * Klug  _et al._ [2023] A. Klug, C. Peters, and A. Forbes, Robust structured light in atmospheric turbulence, Adv. Phot. 5, 016006 (2023).
  * Bachmann  _et al._ [2023] D. Bachmann, M. Isoard, V. Shatokhin, G. Sorelli, N. Treps, and A. Buchleitner, Highly transmitting modes of light in dynamic atmospheric turbulence, Phys. Rev. Lett. 130, 073801 (2023).
  * [96] See Ancillary files for the nummerically-simulated data, an interactive tool for their visualization, and the corresponding Python 3 codes.
  * Andrews and Phillips [2005] L. Andrews and R. Phillips,  _Laser Beam Propagation Through Random Media_ , Online access with subscription: SPIE Digital Library (Society of Photo Optical, 2005).
  * Johnson  _et al._ [1995] N. L. Johnson, S. Kotz, and N. Balakrishnan,  _Continuous Univariate Distributions_ , Vol. 2 (Wiley, 1995).
  * Scarani  _et al._ [2009] V. Scarani, H. Bechmann-Pasquinucci, N. J. Cerf, M. Dušek, N. Lütkenhaus, and M. Peev, The security of practical quantum key distribution, Rev. Mod. Phys. 81, 1301 (2009).
  * Berger and Zhou [2014] V. W. Berger and Y. Zhou, Kolmogorov–Smirnov Test: Overview, in  _Wiley StatsRef: Statistics Reference Online_ (John Wiley & Sons, Ltd, 2014).
  * Martin and Flatté [1988] J. M. Martin and S. M. Flatté, Intensity images and statistics from numerical simulation of wave propagation in 3-d random media, Appl. Opt. 27, 2111 (1988).
  * Kaushal  _et al._ [2017] H. Kaushal, V. K. Jain, and S. Kar,  _Free Space Optical Communication_ (Springer India, New Delhi, 2017).
  * Joanes and Gill [1998] D. N. Joanes and C. A. Gill, Comparing measures of sample skewness and kurtosis, J. R. Stat. Soc.: Series D 47, 183 (1998).
  * Walls [1983] D. F. Walls, Squeezed states of light, Nature (London) 306, 141 (1983).
  * Slusher  _et al._ [1985] R. E. Slusher, L. W. Hollberg, B. Yurke, J. C. Mertz, and J. F. Valley, Observation of squeezed states generated by four-wave mixing in an optical cavity, Phys. Rev. Lett. 55, 2409 (1985).
  * Wu  _et al._ [1986] L.-A. Wu, H. J. Kimble, J. L. Hall, and H. Wu, Generation of squeezed states by parametric down conversion, Phys. Rev. Lett. 57, 2520 (1986).
  * Lane  _et al._ [1992] R. . G. Lane, A. Glindemann, and J. C. Dainty, Simulation of a Kolmogorov phase screen, Waves in Random Media 2, 209 (1992).