# Fourth-order moment of the light field in atmosphere for moderate and strong
turbulence

Roman Baskov  Corresponding address<EMAIL_ADDRESS>Bogolyubov
Institute for Theoretical Physics of the National Academy of Sciences of
Ukraine,
Metrolohichna Street 14-b, Kyiv 03143, Ukraine  Institute of Physics of the
National Academy of Sciences of Ukraine,
prospekt Nauky 46, Kyiv 03028, Ukraine

###### Abstract

Collisionless Boltzmann equation is used to describe the intensity
correlations in partially saturated and fully saturated regimes in terms of
photon distribution function in the phase space. Explicit expression for
fourth moment of the light fields is obtained for the case of moderate and
strong turbulence. Such expression consists of two terms accounting for two
regions in the phase space that independently contribute to the correlation
function. It is shown that present solution agrees with previous results for
fully saturated regime. Additionally it embodies the effect of partially
saturated radiation where the correlations of photon trajectories are
important and the magnitude of the scintillation index is well above the
unity. Fourth moment is used to study the fluctuations of transmittance which
consider the effect of finite detector aperture.

##  I Introduction

Propagation of the light in atmospheric channel is an essential part of many
new-age applied areas such as quantum key distribution Capraro  _et al._
(2012); Usenko  _et al._ (2012), satellite-ground communication Hosseinidehaj
and Malaney (2015); Aspelmeyer  _et al._ (2003), quantum teleportation Ma  _et
al._ (2012); Ren  _et al._ (2017); Hofmann  _et al._ (2019), propagation of
entangled and squeezed states Ursin  _et al._ (2007); Yin  _et al._ (2017);
Peuntinger  _et al._ (2014); Vasylyev  _et al._ (2016). Altogether these areas
contribute to the the development of the next-generation quantum and classical
communication systems including quantum internet Meyers (2015), quantum
protocols Cao  _et al._ (2020); Kish  _et al._ (2020), etc. However, spatio-
temporal properties of the light in atmosphere are modified drastically in a
way of propagation limiting current applications.

Fluctuations of refractive index in Earth’s atmosphere introduce random
distortions to the phase of waves. Since range of sizes of optical
inhomogeneities is very wide, from millimeters (inner scale of turbulence,
$l_{0}$) to hundreds of meters (outer scale of turbulence, $L_{0}$), laser
beam is exposed to the bunch of various negative effects: beam spreading, beam
wandering, fragmentation, beam jitter, intensity fluctuations, etc Fante
(1975, 1980); Tatarskii (1971); Andrews and Phillips (2005). All of them
affect statistical and spatio-temporal properties of the light radiation
causing additional to absorption and scattering losses in atmosphere and
impairing the performance of free-space communication systems.

Although intensity fluctuations plays critical role in the transmission of
optical signal, their description remains one of the most challenging problem
in free-space optics. For atmospheric channel the initially coherent laser
radiation acquires some properties of Gaussian statistics DeWolf (1968). This
leads to saturation effect of the optical wave Gracheva and Gurvich (1965).
Such gradual change of statistical properties complicates theoretical analysis
of correlation properties (fourth-order moment) of propagating radiation.

Some general approach involving equations of evolution for fourth-order
moments was proposed in early studies Fante (1975, 1980). However, the
applicability of these equations is limited due to their complexity. Besides
that correlations of intensity were studied in terms of scintillation index
and covariance function Andrews and Phillips (2005); Andrews  _et al._ (1997);
Dashen (1979); Banakh and Mironov (1979). Nevertheless, there is still no
rigorous theory of fourth moment for optical fields in atmosphere for moderate
and strong turbulence (partially saturated and fully saturated regimes).

Here we use method of photon distribution function (PDF) in the phase space
Berman and Chumak (2006) to describe laser radiation in atmosphere. PDF is
defined as photon density in coordinate-momentum space (phase space).
Intensity of the light and correlation properties of radiation are derived
from PDF moments. Method of PDF was successfully applied to the problem of
light propagation in lossy channels Berman and Chumak (2009); Berman  _et al._
(2007); Chumak and Stolyarov (2013, 2014). Particularly it proved to be
effective for study of intensity fluctuations (scintillation) in the range of
moderate turbulence regime Baskov and Chumak (2018); Chumak and Baskov (2016)
and fourth order moment problem in asymptotic case of large propagation
distances, $z\rightarrow\infty$ Baskov and Chumak (2020). In current paper,
for the description of partially saturated and saturated regimes where only
substantial change of photon momenta due to atmospheric turbulence takes place
we use Boltzmann kinetic equation for PDF Berman and Chumak (2006). In this
case the effect of turbulence is enclosed in random force originating from the
gradient of refractive index.

The paper is devoted to derivation of the expression for fourth moment
applicable for the case of moderate-to-strong and strong turbulence regimes.
Such function plays leading role in many applied researches including
intensity correlation Vellekoop  _et al._ (2010); Newman and Webb (2014),
enhanced focusing Popoff  _et al._ (2014); Vellekoop and Mosk (2007),
different imaging problems Katz  _et al._ (2012); Hardy and Shapiro (2013);
Zhang  _et al._ (2010); Wang  _et al._ (2010); Shi  _et al._ (2013),
reconstruction of the probability distribution of transmittance Semenov and
Vogel (2009) and its applications to communication protocols Xue  _et al._
(2020). In paper we utilize the effect of multiple collisions with turbulent
inhomogeneities which leads to the Gaussian statistics of radiation. In this
case only particular volume in the phase space contribute to the correlation
function.

The remainder of this paper is organized as follows. In Sec. II, we provide
review of photon distribution function approach applied to laser beam
propagation in atmosphere. In Sec. III, explicit expression for fourth moment
of the light fields and its analysis are presented. Section IV devoted to
estimation of transmittance fluctuations and its dependence on the size of
detector aperture. In Appendix A, we give detailed derivation of the
expressions for fourth moment terms.

##  II Preliminaries

### Photon distribution function.

The photon distribution function resembling the idea of distribution functions
in physics of solids Tarasenko  _et al._ (1992) is given by Berman and Chumak
(2006); Rarenko  _et al._ (1992); *ujp2012

$\hat{f}({\bf r},{\bf q},t)=\frac{1}{V}\sum_{\bf k}e^{-i{\bf k\cdot
r}}b^{\dagger}_{{\bf q}+{\bf k}/2}b_{{\bf q}-{\bf k}/2},$ (1)

where $b^{\dagger}_{\bf q}$ and $b_{\bf q}$ are the quantum amplitudes of
bosonic photon field with the wave vector ${\bf q}$; $V\equiv
L_{x}L_{y}L_{z}\equiv SL_{z}$ is the normalizing volume. All operators are
given in the Heisenberg representation. The laser beam propagates in the $z$
direction. It is assumed that ${\bf k}_{\perp}\,,{\bf q}_{\perp}\ll q_{0}$
where $q_{0}$ is the wave vector corresponding to the central frequency
$\omega_{0}$ of radiation $\omega_{0}=cq_{0}$, $c$ is the speed of light in a
vacuum. Such assumption justifies the paraxial approximation. The initial
polarization of light left out of consideration in this case as for a wide
range of propagation distances it remains almost constant (see Ref. Strohbehn
and Clifford (1967)).

The Hamiltonian of photons in a medium with a fluctuating refractive index
could be derived from representation of energy in inhomogeneous media Landau
and Lifshitz (1960)

$H=\sum_{\bf k}\hbar\omega_{\bf k}b^{\dagger}_{\bf k}b_{\bf k}-\sum_{{\bf
k},{\bf k^{\prime}}}\hbar\omega_{\bf k}n_{\bf k^{\prime}}b^{\dagger}_{\bf
k}b_{{\bf k}+{\bf k^{\prime}}}$ (2)

where $\hbar\omega_{\bf k}\equiv\hbar ck$ is the photon energy, and $n_{\bf
k}$ is the Fourier transform of the refractive index fluctuations $\delta
n({\bf r})$. The Fourier transform is defined by

$n_{\bf k}=\frac{1}{V}\int dVe^{i{\bf kr}}\delta n({\bf r}).$ (3)

Usually, $\delta n$ is assumed to be a Gaussian random variable with known
covariance $\langle\delta n({\bf r})\delta n({\bf r}^{\prime})\rangle$. The
covariance is defined by its Fourier transform, $\psi({\bf g})$, with respect
to the difference ${\bf r}-{\bf r}^{\prime}$. In a statistically homogeneous
atmosphere it can be written as

$\langle\delta n({\bf r}-{\bf r}^{\prime})\delta n(0)\rangle=\int d{\bf
g}e^{-i{\bf g}(\bf{r-r^{\prime}})}\psi({\bf g}).$ (4)

The evolution equation for PDF is derived from Heisenberg’s equation. Although
in more general case its evolution is gathered by Boltzmann-Langevin equation
Baskov and Chumak (2018) which takes to account whole range of possible
changes in photon momentum due to collisions with atmosphere inhomogeneities,
for reasonably long distances (see Chumak and Baskov (2016)) description with
collisionless Boltzmann equation

$\partial_{t}\hat{f}({\bf r},{\bf q},t)+{\bf c_{q}}\cdot\partial_{\bf
r}\hat{f}({\bf r},{\bf q},t)+{\bf F}({\bf r})\cdot\partial_{\bf q}\hat{f}({\bf
r},{\bf q},t)=0$ (5)

is justified. In this case the effect of atmospheric turbulence is enclosed in
random smooth force ${\bf F}({\bf r})=\omega_{0}\partial_{\bf r}n({\bf r})$.

The general solution of (5) is obtained by characteristics method

$\hat{f}({\bf r},{\bf q},t)=\phi\Bigg{\\{}{\bf
r}-\int_{0}^{t}dt^{\prime}\frac{\partial{\bf r}(t^{\prime})}{\partial
t^{\prime}};{\bf q}-\int_{0}^{t}dt^{\prime}\frac{\partial{\bf
q}(t^{\prime})}{\partial t^{\prime}}\Bigg{\\}},$ (6)

where evolution of PDF is described in terms of the classical trajectories of
photons

$\displaystyle\frac{\partial{\bf r}(t^{\prime})}{\partial t^{\prime}}={\bf
c}[{\bf q}(t^{\prime})],$ (7) $\displaystyle\frac{\partial{\bf
q}(t^{\prime})}{\partial t^{\prime}}={\bf F}[{\bf r}(t^{\prime})],$ (8)

the function $\phi({\bf r},{\bf q})$ is the “initial” value of $\hat{f}({\bf
r},{\bf q},t)$ in the aperture plane of the source, i.e.

$\phi({\bf r},{\bf q})=\frac{1}{V}\sum_{\bf k}e^{-i{\bf kr}}(b^{+}_{{\bf
q}+\frac{\bf k}{2}}b_{{\bf q}-\frac{\bf k}{2}})|_{t=0}\equiv\sum_{\bf
k}e^{-i{\bf kr}}\phi({\bf k},{\bf q}).$ (9)

Under paraxial approximation atmosphere mostly affects divergence of the beam,
${\bf q}_{\perp}$, and has negligible influence on the longitudinal components
($z$-axis). Then, Eq. (6) can be written as

$\displaystyle\hat{f}({\bf r},{\bf q},t){=}\phi\Bigg{\\{}{\bf r}{-}{\bf
c_{q}}t{+}\frac{c}{q_{0}}\int\limits_{0}^{t}dt^{\prime}t^{\prime}{\bf
F}_{\bot}[{\bf r}(t^{\prime})];\hskip 28.45274pt$ (10) $\displaystyle{\bf
q}{-}\int\limits_{0}^{t}dt^{\prime}{\bf F}_{\bot}[{\bf
r}(t^{\prime})]\Bigg{\\}}.$

Assuming the initial configuration for laser radiation is known, the first and
second moments of $\hat{f}$, which describe beam intensity and its
correlations, could be obtained by means of the iterative procedure for powers
of ${\bf F}_{\perp}({\bf r})$ Berman and Chumak (2006); Chumak and Baskov
(2016).

##  III Intensity correlations

The intensity (density of photons in the spatial domain) is derived from
operator $\hat{f}({\bf r},{\bf q},t)$ by summation over all values of $\bf q$

$\hat{I}({\bf r},t)=\sum_{\bf q}\hat{f}({\bf r},{\bf
q},t)=\frac{1}{V}\sum_{\bf q,k}e^{-i{\bf k\cdot r}}b^{\dagger}_{{\bf q}+{\bf
k}/2}b_{{\bf q}-{\bf k}/2}.$ (11)

Consequently, the intensity correlations, fourth-order moment for the field
operators, is defined by

$\displaystyle\Gamma_{4}({\bf r},{\bf r}^{\prime})$ $\displaystyle\equiv$
$\displaystyle\langle\hat{I}({\bf r},t)\hat{I}({\bf r^{\prime}},t)\rangle$
$\displaystyle{=}$ $\displaystyle\frac{1}{V^{2}}\sum_{{\bf q,k},\atop{\bf
q^{\prime},k^{\prime}}}e^{-i({\bf k}\cdot{\bf r}{+}{\bf k^{\prime}}\cdot{\bf
r^{\prime}})}\langle b^{\dagger}_{{\bf q}{+}\frac{\bf k}{2}}b_{{\bf
q}{-}\frac{\bf k}{2}}b^{\dagger}_{{\bf q^{\prime}}{+}\frac{\bf
k^{\prime}}{2}}b_{{\bf q^{\prime}}{-}\frac{\bf k^{\prime}}{2}}\rangle.$

The averaging $\langle\dots\rangle$ includes the quantum-mechanical averaging
of operators $\hat{I}$ and averaging over different configurations of
atmospheric turbulence. Both averaging can be performed independently.

It was shown in the recent research Baskov and Chumak (2020) that due to
saturation effect for fluctuations in asymptotic case of large distances,
$z\rightarrow\infty$, fourth moment may be expressed via second moments

$\displaystyle\langle b^{\dagger}_{{\bf q}{+}\frac{\bf k}{2}}b_{{\bf
q}{-}\frac{\bf k}{2}}b^{\dagger}_{{\bf q^{\prime}}{+}\frac{\bf
k^{\prime}}{2}}b_{{\bf q^{\prime}}{-}\frac{\bf
k^{\prime}}{2}}\rangle{\approx}\hskip 113.81102pt$ (13) $\displaystyle\langle
b_{{\bf q}{+}\frac{\bf k}{2}}^{\dagger}b_{{\bf q}{-}\frac{\bf
k}{2}}\rangle\langle b_{{\bf q}^{\prime}{+}\frac{{\bf
k}^{\prime}}{2}}^{\dagger}b_{{\bf q}^{\prime}{-}\frac{{\bf
k}^{\prime}}{2}}\rangle{+}\langle b_{{\bf q}+\frac{\bf k}{2}}^{\dagger}b_{{\bf
q}^{\prime}-\frac{{\bf k}^{\prime}}{2}}\rangle\langle b_{{\bf
q}^{\prime}+\frac{{\bf k}^{\prime}}{2}}^{\dagger}b_{{\bf q}{-}\frac{\bf
k}{2}}\rangle$ $\displaystyle=n_{\bf q}n_{{\bf q}^{\prime}}\delta_{{\bf
k},0}\delta_{{\bf k}^{\prime},0}+n_{{\bf q}+\frac{\bf k}{2}}n_{{\bf
q}-\frac{\bf k}{2}}\delta_{{\bf q},{\bf q}^{\prime}}\delta_{{\bf k},-{\bf
k}^{\prime}},$

where $n_{\bf q}\equiv\langle b^{\dagger}_{\bf q}b_{\bf q}\rangle$. Also it is
assumed that initial laser radiation is in a multiphoton coherent state, so
the shot-noise term is omitted. Expression (13) is legitimate if the
amplitudes $b^{\dagger}$ and $b$ obey Gaussian statistics. In other words for
$t\rightarrow\infty$ we assume that each primary coherent electromagnetic wave
experiences multiple scatterings by randomly distributed turbulent eddies
Dashen (1979); DeWolf (1968) and the radiation becomes fully saturated.

At large but finite $z$, in partially saturated regime, one should consider
“nondiagonal” terms in four-wave correlations (13). In this case there are two
regions where pair correlations of the field operators should be taken into
account Berman and Chumak (2006): $(i)$ $k,k^{\prime}\leq R_{b}^{-1}$, $(ii)$
$|{\bf q}-{\bf q}^{\prime}+({\bf k}+{\bf k}^{\prime})/2|,|{\bf q}-{\bf
q}^{\prime}+({\bf k}+{\bf k}^{\prime})/2|\leq R_{b}^{-1}$, where
$R_{b}^{2}\equiv\langle{\bf r}^{2}\rangle_{T}$ (see Baskov and Chumak (2020))
is the ”turbulent” part of beam radius,
$R_{b}^{2}={8z^{3}c\alpha}/{(3r_{0}^{2}\omega_{0}^{2})}$,
$\alpha=0.5{\pi\omega_{0}^{2}}{c}^{-1}\int d{\bf g}g^{2}\psi(g)$. “Turbulent”
term should be dominant to initial radius of beam and diffraction term,
$R_{b}^{2}>r_{0}^{2},4z^{2}q_{0}^{-2}r_{0}^{-2}$, to distinguish curtain level
of saturation of the fluctuations, i.e. partially saturated and fully
saturated regime. Exploiting the approach from Berman and Chumak (2006);
Chumak and Baskov (2016) fourth moments for these two regions are obtained for
Gaussian beams (see Appendix A for details), $\Gamma_{4}({\bf r},{\bf
r}^{\prime})=\Gamma_{4}^{(i)}({\bf r},{\bf r}^{\prime})+\Gamma_{4}^{(ii)}({\bf
r},{\bf r}^{\prime})$:

$\displaystyle\Gamma_{4}^{(i)}({\bf r},{\bf r}^{\prime})=2\pi
C\int\limits_{0}^{\infty}d\tilde{q}\tilde{q}\left[F_{1}^{\tilde{q},\rho_{\parallel}}F_{2}^{\tilde{q},\rho_{\parallel}}\big{(}F_{3}^{\tilde{q},\rho_{\parallel}}H_{1}^{\tilde{q},\rho_{\parallel}}-{G_{1}^{\tilde{q},\rho_{\parallel}}}^{2}\big{)}\big{(}F_{4}^{\tilde{q},\rho_{\parallel}}H_{2}^{\tilde{q},\rho_{\parallel}}-{G_{2}^{\tilde{q},\rho_{\parallel}}}^{2}\big{)}\right]^{-\frac{1}{2}}\times\hskip
85.35826pt$ (14)
$\displaystyle\exp\left\\{-\frac{\left(\tilde{q}-\rho_{\parallel}\frac{G_{1}^{\tilde{q},\rho_{\parallel}}}{2F_{3}^{\tilde{q},\rho_{\parallel}}}\right)^{2}}{(H_{1}^{\tilde{q},\rho_{\parallel}}-{G_{1}^{\tilde{q},\rho_{\parallel}}}^{2}/F_{3}^{\tilde{q},\rho_{\parallel}})}\right\\}\exp\left\\{-\frac{\rho_{\perp}^{2}H_{2}^{\tilde{q},\rho_{\parallel}}}{(H_{2}^{\tilde{q},\rho_{\parallel}}F_{4}^{\tilde{q},\rho_{\parallel}}-{G_{2}^{\tilde{q},\rho_{\parallel}}}^{2})}\right\\}\exp\left\\{-\left(\frac{\rho_{\parallel}^{\prime
2}}{F_{1}^{\tilde{q},\rho_{\parallel}}}+\frac{\rho_{\perp}^{\prime
2}}{F_{2}^{\tilde{q},\rho_{\parallel}}}+\frac{\rho_{\parallel}^{2}}{4F_{3}^{\tilde{q},\rho_{\parallel}}}\right)\right\\},$
$\displaystyle\Gamma_{4}^{(ii)}({\bf r},{\bf r}^{\prime})=2\pi
C\int\limits_{0}^{\infty}d\tilde{q}\tilde{q}\left[F_{1}^{\tilde{q}}F_{2}^{\tilde{q}}\big{(}F_{3}^{\tilde{q}}H_{1}^{\tilde{q}}-{G_{1}^{\tilde{q}}}^{\,2}\big{)}\big{(}F_{4}^{\tilde{q}}H_{2}^{\tilde{q}}-{G_{2}^{\tilde{q}}}^{\,2}\big{)}\right]^{-\frac{1}{2}}\times\hskip
184.9429pt$ (15)
$\displaystyle\exp\left\\{-\frac{{\tilde{q}}^{2}}{H_{1}^{\tilde{q}}-{G_{1}^{\tilde{q}}}^{\,2}/F_{3}^{\tilde{q}}}\right\\}\exp\left\\{-\left(\frac{\rho_{\parallel}^{\prime
2}}{F_{1}^{\tilde{q}}}+\frac{\rho_{\perp}^{\prime
2}}{F_{2}^{\tilde{q}}}\right)\right\\}\exp\left\\{-i2{\tilde{q}}\rho_{\parallel}\frac{q_{0}}{z}\right\\},$

where $\boldsymbol{\rho}={\bf r}-{\bf r}^{\prime}$,
$\boldsymbol{\rho}^{\prime}=({\bf r}+{\bf r}^{\prime})/2$ are two-dimensional
vectors transverse to propagation direction $z$; $F$, $G$, $H$ are functions
of $\rho_{\parallel}$ and ${\tilde{q}}$ for $(i)$ and functions of
${\tilde{q}}$ for $(ii)$ contributions to $\Gamma_{4}$; $C$ is constant
derived from total flux. It is worth to emphasize that although intensity
correlations are evaluated for two points ${\bf r}$ and ${\bf r}^{\prime}$,
the expression for fourth moment is expressed via the difference ${\bf r}-{\bf
r}^{\prime}$ and the center of mass position $({\bf r}+{\bf r}^{\prime})/2$.
The dependence on $\boldsymbol{\rho}^{\prime}$ has a simple Gaussian form. At
the same time dependence on $\boldsymbol{\rho}$ is quite intricate and
accounts for the correlations of different trajectories. Such a situation
could be favorable for the calculation of integral quantities that consider
spatial distribution of the radiation in detector aperture plane.

Fourth moment $\Gamma_{4}({\bf r},{\bf r^{\prime}})$ characterize
spatiotemporal properties of the laser beam in $(x,y)$-plane in atmosphere.
Functions $F$, $G$, $H$ (see Appendix A) incorporate the effect of correlation
for different photon trajectories with $\\{{\bf r},{\bf q}\\}$ and $\\{{\bf
r}^{\prime},{\bf q}^{\prime}\\}$ on intensity correlations. As it was shown in
Berman and Chumak (2006); Chumak and Baskov (2016) such correlations of
trajectories are responsible for intensity fluctuations in the range of
moderate and strong turbulence.

For the case of asymptotically large distances, $z\rightarrow\infty$, cross-
correlation term vanishes because ${\Delta{r}}\rightarrow\infty$ (see (A)) due
to randomization of the particle displacements from the straight lines, so
functions $F$, $G$, and $H$ do not depend on ${\tilde{q}}$ and $\rho$.
Therefore, the values of functions are expressed via $R_{b}$:
$F\approx\frac{1}{2}R_{b}^{2}$, $G\approx\frac{3}{4}R_{b}^{2}$,
$H\approx\frac{3}{2}R_{b}^{2}$. (Free-space terms are omitted since turbulence
is assumed to give dominant contribution.) Also since cross-correlation is
vanished the integration over directions for $\tilde{\bf q}$ is preserved

$\displaystyle\Gamma_{4}^{(i)}({\bf r},{\bf r}^{\prime})=C\iint d\tilde{\bf
q}\left[F\big{(}FH-{G}^{2}\big{)}\right]^{-{1}}\times\hskip 28.45274pt$ (16)
$\displaystyle\exp\left\\{-\frac{\left(\tilde{\bf q}-({\bf r}-{\bf
r^{\prime}})\frac{G}{2F}\right)^{2}}{(H-{G}^{2}/F)}\right\\}\exp\left\\{-\left(\frac{{{\bf
r}}^{2}}{2F}+\frac{{{\bf r^{\prime}}}^{2}}{2F}\right)\right\\},$
$\displaystyle\Gamma_{4}^{(ii)}({\bf r},{\bf r}^{\prime})=\hskip 170.71652pt$
(17) $\displaystyle C\iint d\tilde{\bf
q}\left[F\big{(}FH-{G}^{2}\big{)}\right]^{-{1}}\exp\left\\{-\frac{\tilde{\bf
q}^{2}}{(H-{G}^{2}/F)}\right\\}\times$
$\displaystyle\exp\left\\{-\left(\frac{({\bf r}+{\bf
r^{\prime}})^{2}}{4F}\right)\right\\}\exp\left\\{-i2\tilde{\bf q}({\bf r}-{\bf
r^{\prime}})\frac{q_{0}}{z}\right\\},$

In this case it is easy to perform integration analytically. The contribution
$\Gamma_{4}^{(i)}({\bf r},{\bf r}^{\prime})$ can be expressed via average
intensity, $\langle\hat{I}({\bf
r})\rangle\propto\exp\left\\{-\left(\frac{{{\bf
r}}^{2}}{R_{b}^{2}}\right)\right\\}\big{/}R_{b}^{2}$ (see Baskov and Chumak
(2020)) and $\Gamma_{4}^{(ii)}({\bf r},{\bf r}^{\prime})$ has a simple
Gaussian form, so

$\displaystyle\Gamma_{4}({\bf r},{\bf r}^{\prime})=\langle\hat{I}({\bf
r},t)\rangle\langle\hat{I}({\bf r^{\prime}},t)\rangle+$ (18) $\displaystyle
C\pi\left(\frac{1}{F}\right)^{2}\exp\left[{-}\frac{({\bf r}{+}{\bf
r^{\prime}})^{2}}{4F}{-}\frac{({\bf r}{-}{\bf
r^{\prime}})^{2}q_{0}^{2}(H-G^{2}/F)}{z^{2}}\right],$

which is exactly the result of Baskov and Chumak (2020), taking to account
that $2F\approx R_{b}^{2}$, $q_{0}^{2}(H-G^{2}/F)/z^{2}\approx\langle{\bf
q}^{2}\rangle_{T}/8$, $\langle{\bf q}^{2}\rangle_{T}=4\alpha t$, and the
relation between constants.

### Applicability of approximation

The applicability of expressions (14) and (15) is inherent from main
approximations of the approach, i.e., collisionless Boltzmann equation and the
concept of photon trajectories. As it was pointed out in Berman and Chumak
(2006) the momentum of photons should be much bigger than characteristic wave
vectors of turbulence. Since upper limit of the spectrum is defined by inner
scale of turbulence, $l_{0}$, the relation $\langle{\bf
q}^{2}\rangle_{T}l_{0}^{2}$ should be large enough. On the other hand,
throughout the paper we account for the effect of correlation of photon
trajectories, so such concept should be justified. To consider photons as
particles, whose density in the $({\bf r},{\bf q})$ domain is defined by the
distribution function $\hat{f}({\bf r},{\bf q},t)$, the uncertainty of the
momentum ${\bf q}$ should be small Chumak and Baskov (2016). The value of the
uncertainty can be estimated from the definition of the distribution function
(1) as ${\bf k}/2$. For large distances and Tatarskii spectrum such
uncertainty is estimated by the relation

$\langle{\bf q}^{2}\rangle_{T}R_{b}^{2}\approx 15\cdot
q_{0}^{2}l_{0}^{-2/3}C_{n}^{4}z^{4},$ (19)

which also should be large enough.

##  IV Fluctuations of transmitted radiation

For many practical cases, e.g., development classical and quantum
communication, the fluctuation of transmittance in Earth’s atmosphere is a key
parameter that defines the properties of atmospheric channel Vasylyev  _et
al._ (2019); Barrios and Dios (2012); Vetelino  _et al._ (2007). The magnitude
of fluctuations is estimated via variance

$\sigma_{\eta}^{2}=\frac{\langle\hat{\eta}^{2}\rangle-\langle\hat{\eta}\rangle^{2}}{\langle\hat{\eta}\rangle^{2}},$
(20)

where transmittance of the optical channel is defined as

$\hat{\eta}=(4C\pi^{3})^{-\frac{1}{2}}\int\limits_{\mathcal{A}}d{\bf
r}\hat{I}({\bf r},t)$ (21)

and accounts for the finite size of detector aperture. The normalizing
condition $\langle\hat{\eta}\rangle=1$ for $\mathcal{A}$ which is much larger
than the beam cross section is used.

Figure 1: Effect of “nondiagonal” terms in $\Gamma_{4}$ on aperture averaged
scintillations. In each pair of lines lower curve depicts the results from
Baskov and Chumak (2020), where it is assumed that intensity fluctuations are
fully saturated; upper curve represents the results for (14) and (15). Dash-
dotted lines: $z=20\,km$, $C_{n}^{2}=2.5\times 10^{-14}\,m^{-2/3}$, Rytov
variance $\sigma_{R}^{2}=62$; dashed lines: $z=17\,km$, $C_{n}^{2}=5.8\times
10^{-15}\,m^{-2/3}$, $\sigma_{R}^{2}=60$; solid lines: $z=100\,km$,
$C_{n}^{2}=2.5\times 10^{-16}\,m^{-2/3}$, $\sigma_{R}^{2}=66$. Common
parameters of the beam and channel for all curves: $r_{0}=0.01\,m$,
$l_{0}/2\pi=10^{-3}\,m$, $q_{0}=10^{7}\,m^{-1}$.

Two moments for $\eta$ are defined as Fried (1967); Vasylyev  _et al._ (2018)

$\displaystyle\langle\hat{\eta}\rangle$ $\displaystyle=$
$\displaystyle(4C\pi^{3})^{-\frac{1}{2}}\int\limits_{\mathcal{A}}d{\bf
r}\Gamma_{2}({\bf r}),$ (22) $\displaystyle\langle\hat{\eta}^{2}\rangle$
$\displaystyle=$ $\displaystyle(4C\pi^{3})^{-1}\int\limits_{\mathcal{A}}d{\bf
r}\int\limits_{\mathcal{A}}d{\bf r}^{\prime}\Gamma_{4}({\bf r},{\bf
r}^{\prime}),$ (23)

where $\Gamma_{2}({\bf r})\equiv\langle\hat{I}({\bf r},t)\rangle$. To obtain
the fluctuations of transmitted radiation $\sigma_{\eta}^{2}$ we calculate
numerically (five-fold integration)
$\langle\hat{\eta}^{2}\rangle=\langle\hat{\eta}^{2}\rangle^{(i)}+\langle\hat{\eta}^{2}\rangle^{(ii)}$
for circular aperture with radius $R$ using the expressions (14) and (15).

Figure 2:  Contributions of two regions $(i)$ (solid lines) and $(ii)$ (dashed
lines) for partially saturated (upper pair of curves) and saturated regimes
(lower pair of curves) for $\langle{\bf q}^{2}\rangle_{T}^{1/2}=600\,m^{-1}$
from Fig. 1.

First of all, it is informative to compare approximations where partially
saturated and fully saturated regimes were assumed (Fig. 1) for atmospheric
channels considered in Baskov and Chumak (2020). Such comparison allows both
to estimate the accuracy of asymptotic approximation, $z\rightarrow\infty$,
and to emphasize the differentiation of correlation properties for two
regimes. For partially saturated radiation scintillations are slightly larger
because of the additional contribution of “nondiagonal” terms in fourth moment
which incorporate larger region of the phase space and account for residual
effect of the initial statistics of radiation. Remarkably, the values of
aperture averaged scintillations differ from unity for point-like aperture.
That is a natural outcome for partially saturated regime where radiation still
preserve some properties of initial statistics and does not fully acquire
Gaussian statistics. In contrast to approximation of fully saturated radiation
there is clear dependence of the values of scintillations,
$\sigma_{\eta}^{2}$, for small aperture sizes on Rytov parameter. In
consistent with other studies, there is such size of aperture where detector
could not be considered as point-like and steep reduction of the fluctuations
is observed. These values strongly dependent on the transverse momentum of
photons $\langle{\bf q}^{2}\rangle_{T}$. In addition, unlikely to fully
saturated approximation the beam spreading $R_{b}^{2}$ also plays significant
role to the behavior of aperture-averaged scintillations via both (14) and
(15). Figure 2 shows that present result more adequately account for the
values of $\Gamma_{4}^{(i)}$ and corresponding correlation length. Since
Baskov and Chumak (2020) considers fully saturated regime, $\Gamma_{4}^{(i)}$
does not contribute to the fluctuations there. In contrast term (14) have
sizable effect till the sizes of aperture are less than beam radius.

Also, there are two basic properties for scintillations averaging that are
preserved in current approximation: for point-like detectors, ${\bf r}={\bf
r^{\prime}}=0$, (14) and (15) contribute equally which repeats the result of
previous works Berman and Chumak (2006); Baskov and Chumak (2020); for
aperture sizes reasonably larger than beam radius the fluctuations of
transmittance tend to zero.

Figure 3: Aperture-averaged scintillation index vs. radius of detector
aperture. Parameters for the both graphs: $r_{0}=0.01\,m$,
$l_{0}/2\pi=10^{-3}\,m$, $q_{0}=1.29\times 10^{7}\,m^{-1}$, $z=3\,km$.

For atmospheric channels under partially saturated regime (Fig. 3) asymptotic
result (18) significantly differs from one that considers “nondiagonal” terms
in the expression for $\Gamma_{4}$. First of all, $\Gamma_{4}^{(i)}$ term
contributes substantially to the values of transmittance fluctuations.
Particularly it is responsible for long “tails” of the curves at large
detector apertures. This effect is reminiscent to the leveling effect
mentioned in Churnside (1991), where there are two characteristic scales,
$\rho_{0}\propto\langle{\bf q}^{2}\rangle_{T}^{1/2}$ and
$\rho_{0}z/q_{0}\propto R_{b}$, which defines correlation properties of the
light radiation. Particularly one may see from Fig. 3 that correlation length
in term $\Gamma_{4}^{(ii)}$ being proportional to $\langle{\bf
q}^{2}\rangle_{T}^{1/2}$ is much smaller than correlation length of
$\Gamma_{4}^{(i)}$ which comes from the values of beam radius.

Figure 4: Aperture-averaged scintillation index vs. Rytov parameter.
Parameters of the channel: $r_{0}=0.01\,m$, $l_{0}/2\pi=10^{-3}\,m$,
$q_{0}=10^{7}\,m^{-1}$. Figure 5: Aperture-averaged scintillation index vs.
Rytov parameter, $q_{0}=1.29\times 10^{7}\,m^{-1}$, other parameters are the
same as in Fig. 3.

Figures 4 and 5 depict aperture effect on scintillations in relation to values
of Rytov parameter. One might see that both magnitude of fluctuations and
steepness of its decrease with larger sizes of the detector aperture strongly
depends on parameters of the atmospheric channel. Generally, for bigger values
of Rytov parameter the effect of saturation of the fluctuations responsible
for decrease of $\sigma_{\eta}^{2}$. That is, magnitude of intensity
fluctuations decrease for atmospheric channels with stronger turbulence
regime.

##  V Summary and Conclusions

The fourth moment of light fields in atmosphere was derived under
collisionless Boltzmann equation approximation where smooth random force
represents the effect of turbulence on the laser beam. Such approximation is
justified for the case of moderate and strong turbulence where fluctuations
are saturated for particular degree. Randomization of the electromagnetic
waves allows to distinguish two separate regions in the phase space that
contribute to noise level in optical system and its correlation properties.
For the case of partially saturated regime intensity correlations are not
factorized to product of two-wave correlations and additional “nondiagonal”
terms should be considered. For Gaussian beams the effect of photon
trajectories correlations is enclosed in a set of functions that define
distribution of photon density in the phase space. Values of such functions
strongly depend on properties of the atmospheric channel.

The fourth moment is applied for the problem of transmittance fluctuations. It
is shown that current calculations describe more wide range of parameters of
the atmospheric channels. There are two characteristic scales whose relation
defines correlation length for detected radiation. Found results could be
useful for practical purposes.

## Acknowledgments

The author thanks D. Vasylyev, A. Semenov, and O. O. Chumak for useful
discussions and comments. The author acknowledges support from National
Research Foundation of Ukraine through the Project 2020.02/0111.

##  Appendix A Calculation of $\Gamma_{4}$

In the main text we consider Gaussian beams and Tatarskii spectrum for the
fluctuations of index-of-refraction. Therefore, for Gaussian beam the initial
configuration of radiation defines Baskov and Chumak (2020)

$\langle\phi({\bf k},{\bf q})\rangle_{qm}=\frac{2\pi
r_{0}^{2}}{VL_{x}L_{y}}\langle
b^{\dagger}b\rangle_{qm}e^{-(q_{\perp}^{2}+k_{\perp}^{2}/4)r_{0}^{2}/2},$ (24)

where the symbol $\langle...\rangle_{qm}$ indicates a quantum-mechanical
averaging of operators in the angle brackets. Coefficient $C^{\prime}$ can be
obtained if the total photon flux is known. Tatarskii spetrum is defined as

$\psi({\bf
g})=0.033C_{n}^{2}{\exp\left\\{-(gl_{0}^{\prime})^{2}\right\\}}{g^{-11/3}}.$
(25)

In the similar manner to the works Berman and Chumak (2006); Chumak and Baskov
(2016) we can write down the general expression for fourth moment in terms of
random fluctuating force ${\bf F}_{\perp}$

$\displaystyle\langle I({\bf r},t)I({\bf
r^{\prime}},t)\rangle{=}\Bigg{(}\frac{2\pi
r_{0}^{2}}{VL_{x}L_{y}}\Bigg{)}^{2}\langle
b^{\dagger}b^{\dagger}bb\rangle_{qm}\sum_{\begin{subarray}{c}{\bf q},{\bf
k},\\\ {\bf q^{\prime}},{\bf k^{\prime}}\end{subarray}}\bigg{\langle}e^{-i{\bf
k}[{\bf r}-{\bf c_{q}}t+\frac{c}{q_{0}}\int_{0}^{t}dt^{\prime}t^{\prime}{\bf
F}_{\perp}({\bf r}({\bf q},t^{\prime}))]-i{\bf k^{\prime}}[{\bf
r^{\prime}}-{\bf c_{\bf
q^{\prime}}}t+\frac{c}{q_{0}}\int_{0}^{t}dt^{\prime}t^{\prime}{\bf
F}_{\perp}({\bf r}({\bf q^{\prime}},t^{\prime}))]}\times$
$\displaystyle\exp{\left(-\left(Q^{2}+Q^{\prime
2}+\frac{k^{2}+{k^{\prime}}^{2}}{4}\right)\frac{r_{0}^{2}}{2}\right)}\bigg{\rangle},$
(26)

where ${\bf Q}={\bf q}-\int_{0}^{t}dt^{\prime}{\bf F}_{\perp}({\bf r}({\bf
q},t))$ and ${\bf Q^{\prime}}={\bf q^{\prime}}-\int_{0}^{t}dt^{\prime}{\bf
F}_{\perp}({\bf r^{\prime}}({\bf q^{\prime}},t))$ are the solutions of
evolution equations for momenta. It is reasonable to express all turbulent
part in the exponent in linear form for ${\bf F}_{\perp}$. So the factor
$e^{-(Q_{\perp}^{2}+{Q_{\perp}^{\prime}}^{2})\frac{r_{0}^{2}}{2}}$ in (A) is
expressed in the integral form as

$e^{-(Q_{\perp}^{2}+{Q^{\prime}_{\perp}}^{2})\frac{r_{0}^{2}}{2}}=\int\frac{d{\bf
p}d{\bf p^{\prime}}}{(2\pi r_{0}^{2})^{2}}e^{i{\bf p}\cdot{\bf
Q}_{\perp}+i{\bf p^{\prime}}\cdot{\bf
Q}_{\perp}^{\prime}-(p^{2}+{p^{\prime}}^{2})/{2r_{0}^{2}}}.$ (27)

In this case (A) can be expressed as

$\displaystyle\langle I({\bf r},t)I({\bf r^{\prime}},t)\rangle{=}\frac{\langle
b^{\dagger}b^{\dagger}bb\rangle_{qm}}{(VL_{x}L_{y})^{2}}\int{d{\bf p}d{\bf
p^{\prime}}}\sum_{\begin{subarray}{c}{\bf q},{\bf k},\\\ {\bf q^{\prime}},{\bf
k^{\prime}}\end{subarray}}e^{-i{\bf k}[{\bf r}-{\bf c_{q}}t]-i{\bf
k^{\prime}}[{\bf r^{\prime}}-{\bf c_{\bf
q^{\prime}}}t]}e^{-{(k^{2}+{k^{\prime}}^{2})r_{0}^{2}}/{8}}e^{i{\bf
p}\cdot{\bf Q}_{\perp}+i{\bf p^{\prime}}\cdot{\bf
Q}_{\perp}^{\prime}-(p^{2}+{p^{\prime}}^{2})/{2r_{0}^{2}}}\langle M\rangle,$
(28)

where the factor

$M{=}\exp\left({-}i\int_{0}^{t}dt^{\prime}\\{({\bf p}{+}{\bf
k}t^{\prime}\frac{c}{q_{0}}){\bf F}[{\bf r}({\bf q},t^{\prime})]{+}({\bf
p}^{\prime}{+}{\bf k}^{\prime}t^{\prime}\frac{c}{q_{0}}){\bf F}[{\bf r}({\bf
q}^{\prime},t^{\prime})]\\}\right)$ (29)

includes all fluctuating part in fourth moment. (For the sake of brevity from
herein we omit $\perp$ notation for perpendicular to propagation direction.)
Using similar to Berman and Chumak (2006); Chumak and Baskov (2016) approach,
the average of $M$ is expressed as

$\displaystyle\langle
M\rangle=\exp\Bigg{\\{}-0.033C_{n}^{2}\pi^{2}q_{0}^{2}\int_{0}^{z}dx\int d{\bf
g}g^{-11/3}e^{-g^{2}l_{0}^{\prime 2}}\Big{[}(({\bf p}+{\bf k}x/q_{0})\cdot{\bf
g})^{2}+(({\bf p}^{\prime}+{\bf k}^{\prime}x/q_{0})\cdot{\bf g})^{2}$
$\displaystyle+2({\bf p}+{\bf k}x/q_{0})\cdot{\bf g}({\bf p}^{\prime}+{\bf
k}^{\prime}x/q_{0})\cdot{\bf g}\exp\left\\{i{\bf g}\cdot\boldsymbol{\Delta
r}-\frac{R_{b}^{2}}{480{l_{0}^{\prime}}^{2}}{\Delta
r}^{2}(1-x/z)^{3}g^{2}\left[1+\frac{R_{b}^{2}(1-x/z)^{3}}{672{l_{0}^{\prime}}^{2}}\right]\right\\}\Big{]}\Bigg{\\}},$
(30)

where $g=|{\bf g}|$ and we also used Markov approximation index-of-refraction
spectrum is delta-correlated in the direction of propagation, which was
rigorously justified in Chumak and Baskov (2016). After integration over
direction of ${\bf g}$

$\displaystyle\langle
M\rangle=\exp\Bigg{\\{}-0.033C_{n}^{2}\pi^{2}q_{0}^{2}\int_{0}^{z}dx\int_{0}^{\infty}dgg^{-2/3}e^{-g^{2}l_{0}^{\prime
2}}\Big{[}({\bf p}+{\bf k}x/q_{0})^{2}+({\bf p}^{\prime}+{\bf
k}^{\prime}x/q_{0})^{2}\hskip 56.9055pt$
$\displaystyle+\big{(}2(p_{\parallel}+k_{\parallel}x/q_{0})(p_{\parallel}^{\prime}+k_{\parallel}^{\prime}x/q_{0})(J_{0}(g\Delta
r)-J_{2}(g\Delta
r))+2(p_{\perp}+k_{\perp}x/q_{0})(p_{\perp}^{\prime}+k_{\perp}^{\prime}x/q_{0})(J_{0}(g\Delta
r)+J_{2}(g\Delta r))\big{)}$ (31)
$\displaystyle\times\exp\left\\{-\frac{R_{b}^{2}}{480{l_{0}^{\prime}}^{2}}{\Delta
r}^{2}(1-x/z)^{3}g^{2}\left[1+\frac{R_{b}^{2}(1-x/z)^{3}}{672{l_{0}^{\prime}}^{2}}\right]\right\\}\Big{]}\Bigg{\\}},$

where $J_{0}$ and $J_{2}$ are zeroth and second order Bessel functions, vector
$\boldsymbol{\Delta r}=({\bf r}-{\bf r^{\prime}})-({\bf q}-{\bf
q}^{\prime})(z-x)/q_{0}$, $\Delta r=|\boldsymbol{\Delta r}|$. The indices
$\\{_{\parallel}\\}$ and $\\{_{\perp}\\}$ indicate the parallel and
perpendicular to $\boldsymbol{\Delta r}$ components of the corresponding 2D
vectors. First two terms in square brackets represent correlation of waves
with same pairs $\\{{\bf r},{\bf q}\\}$ and $\\{{\bf r}^{\prime},{\bf
q}^{\prime}\\}$. The terms that include Bessel functions consider cross-
correlation of photon trajectories ($\\{{\bf r},{\bf q}\\}$ and $\\{{\bf
r}^{\prime},{\bf q}^{\prime}\\}$ are different) which were extensively
reviewed in Chumak and Baskov (2016). In (A) and (A) we omit the dependence on
angle between $\boldsymbol{\Delta r}$ and ${\bf g}$ (see Eq. (40) in Chumak
and Baskov (2016)) in the last exponent as its contribution is small. After
substitution of (A) to (28) and change of variables, most of integrations
could be performed analytically. As a result expression for fourth moment can
be written as

$\displaystyle\Gamma_{4}^{(i)}({\bf r},{\bf r}^{\prime})=2\pi
C\int\limits_{0}^{\infty}d\tilde{q}\tilde{q}\left[F_{1}^{\tilde{q},\rho_{\parallel}}F_{2}^{\tilde{q},\rho_{\parallel}}\big{(}F_{3}^{\tilde{q},\rho_{\parallel}}H_{1}^{\tilde{q},\rho_{\parallel}}-{G_{1}^{\tilde{q},\rho_{\parallel}}}^{2}\big{)}\big{(}F_{4}^{\tilde{q},\rho_{\parallel}}H_{2}^{\tilde{q},\rho_{\parallel}}-{G_{2}^{\tilde{q},\rho_{\parallel}}}^{2}\big{)}\right]^{-\frac{1}{2}}\times\hskip
85.35826pt$ (32)
$\displaystyle\exp\left\\{-\frac{\left(\tilde{q}-\rho_{\parallel}\frac{G_{1}^{\tilde{q},\rho_{\parallel}}}{2F_{3}^{\tilde{q},\rho_{\parallel}}}\right)^{2}}{(H_{1}^{\tilde{q},\rho_{\parallel}}-{G_{1}^{\tilde{q},\rho_{\parallel}}}^{2}/F_{3}^{\tilde{q},\rho_{\parallel}})}\right\\}\exp\left\\{-\frac{\rho_{\perp}^{2}H_{2}^{\tilde{q},\rho_{\parallel}}}{(H_{2}^{\tilde{q},\rho_{\parallel}}F_{4}^{\tilde{q},\rho_{\parallel}}-{G_{2}^{\tilde{q},\rho_{\parallel}}}^{2})}\right\\}\exp\left\\{-\left(\frac{\rho_{\parallel}^{\prime
2}}{F_{1}^{\tilde{q},\rho_{\parallel}}}+\frac{\rho_{\perp}^{\prime
2}}{F_{2}^{\tilde{q},\rho_{\parallel}}}+\frac{\rho_{\parallel}^{2}}{4F_{3}^{\tilde{q},\rho_{\parallel}}}\right)\right\\},$

functions $F$, $G$ and $H$ are defined as

$\displaystyle
F_{1,2,3,4}=\frac{r_{0}^{2}}{4}+\frac{z^{2}}{q_{0}^{2}r_{0}^{2}}+\varphi_{1,2,3,4}$
(33) $\displaystyle G_{1,2}=\frac{z^{2}}{q_{0}^{2}r_{0}^{2}}+\gamma_{1,2}$
(34) $\displaystyle H_{1,2}=\frac{z^{2}}{q_{0}^{2}r_{0}^{2}}+\chi_{1,2}$ (35)

where

$\displaystyle\varphi_{1}{=}\alpha\int\limits_{0}^{1}d\tau\tau^{2}[1{+}(1{+}\beta{\Delta\tilde{r}}^{2}\tau^{3})^{-\frac{1}{6}}[{{}_{1}F_{1}}(\frac{1}{6},1;\frac{-{\Delta\tilde{r}}^{2}}{4{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})}){+}\frac{{\Delta\tilde{r}}^{2}}{48{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})}{{}_{1}F_{1}}(\frac{7}{6},3;\frac{-{\Delta\tilde{r}}^{2}}{{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})})]]$
(36)
$\displaystyle\varphi_{2}{=}\alpha\int\limits_{0}^{1}d\tau\tau^{2}[1{+}(1{+}\beta{\Delta\tilde{r}}^{2}\tau^{3})^{-\frac{1}{6}}[{{}_{1}F_{1}}(\frac{1}{6},1;\frac{-{\Delta\tilde{r}}^{2}}{4{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})}){-}\frac{{\Delta\tilde{r}}^{2}}{48{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})}{{}_{1}F_{1}}(\frac{7}{6},3;\frac{-{\Delta\tilde{r}}^{2}}{{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})})]]$
(37)
$\displaystyle\varphi_{3}{=}\alpha\int\limits_{0}^{1}d\tau\tau^{2}[1{-}(1{+}\beta{\Delta\tilde{r}}^{2}\tau^{3})^{-\frac{1}{6}}[{{}_{1}F_{1}}(\frac{1}{6},1;\frac{-{\Delta\tilde{r}}^{2}}{4{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})}){-}\frac{{\Delta\tilde{r}}^{2}}{48{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})}{{}_{1}F_{1}}(\frac{7}{6},3;\frac{-{\Delta\tilde{r}}^{2}}{{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})})]]$
(38)
$\displaystyle\varphi_{4}{=}\alpha\int\limits_{0}^{1}d\tau\tau^{2}[1{-}(1{+}\beta{\Delta\tilde{r}}^{2}\tau^{3})^{-\frac{1}{6}}[{{}_{1}F_{1}}(\frac{1}{6},1;\frac{-{\Delta\tilde{r}}^{2}}{4{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})}){+}\frac{{\Delta\tilde{r}}^{2}}{48{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})}{{}_{1}F_{1}}(\frac{7}{6},3;\frac{-{\Delta\tilde{r}}^{2}}{{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})})]]$
(39)
$\displaystyle\gamma_{1}{=}\alpha\int\limits_{0}^{1}d\tau\tau[1{-}(1{+}\beta{\Delta\tilde{r}}^{2}\tau^{3})^{-\frac{1}{6}}[{{}_{1}F_{1}}(\frac{1}{6},1;\frac{-{\Delta\tilde{r}}^{2}}{4{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})}){-}\frac{{\Delta\tilde{r}}^{2}}{48{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})}{{}_{1}F_{1}}(\frac{7}{6},3;\frac{-{\Delta\tilde{r}}^{2}}{{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})})]]$
(40)
$\displaystyle\gamma_{2}{=}\alpha\int\limits_{0}^{1}d\tau\tau[1{-}(1{+}\beta{\Delta\tilde{r}}^{2}\tau^{3})^{-\frac{1}{6}}[{{}_{1}F_{1}}(\frac{1}{6},1;\frac{-{\Delta\tilde{r}}^{2}}{4{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})}){+}\frac{{\Delta\tilde{r}}^{2}}{48{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})}{{}_{1}F_{1}}(\frac{7}{6},3;\frac{-{\Delta\tilde{r}}^{2}}{{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})})]]$
(41)
$\displaystyle\chi_{1}{=}\alpha\int\limits_{0}^{1}d\tau[1{-}(1{+}\beta{\Delta\tilde{r}}^{2}\tau^{3})^{-\frac{1}{6}}[{{}_{1}F_{1}}(\frac{1}{6},1;\frac{-{\Delta\tilde{r}}^{2}}{4{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})}){-}\frac{{\Delta\tilde{r}}^{2}}{48{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})}{{}_{1}F_{1}}(\frac{7}{6},3;\frac{-{\Delta\tilde{r}}^{2}}{{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})})]]$
(42)
$\displaystyle\chi_{2}{=}\alpha\int\limits_{0}^{1}d\tau[1{-}(1{+}\beta{\Delta\tilde{r}}^{2}\tau^{3})^{-\frac{1}{6}}[{{}_{1}F_{1}}(\frac{1}{6},1;\frac{-{\Delta\tilde{r}}^{2}}{4{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})}){+}\frac{{\Delta\tilde{r}}^{2}}{48{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})}{{}_{1}F_{1}}(\frac{7}{6},3;\frac{-{\Delta\tilde{r}}^{2}}{{l_{0}^{\prime}}^{2}(1+\beta{\Delta\tilde{r}}^{2}\tau^{2})})]],$
(43)

and parameters are defined as $\alpha=3R_{b}^{2}$,
$\beta=\frac{R_{b}^{2}}{480{l_{0}^{\prime}}^{2}}\left[1+\frac{R_{b}^{2}\tau^{3}}{672{l_{0}^{\prime}}^{2}}\right]$,
${\Delta\tilde{r}}=\rho_{\parallel}-2\tilde{q}\tau$ and
${}_{1}F_{1}(a,b;z)=\sum\limits_{n=0}^{\infty}\frac{a^{(n)}z^{n}}{b^{(n)}n!}$
is a confluent hypergeometric function (Kummer’s function), $a^{(n)}$,
$b^{(n)}$ are the Pochhammer symbols. Here notations $\boldsymbol{\rho}={\bf
r}-{\bf r}^{\prime}$, $\boldsymbol{\rho}^{\prime}=({\bf r}+{\bf
r}^{\prime})/2$ are used since ${\bf r}$ and ${\bf r}^{\prime}$ enters
expression for $\Gamma_{4}$ only in such combinations (see more details in
main text). Functions (36)-(43) depend on $\rho_{\parallel}$ and $\tilde{q}$.

For the region of the phase space defined with conditions (ii) one should
perform change of indices similar to one in Berman and Chumak (2009); Baskov
and Chumak (2020)

${\bf k}\rightarrow{\bf q}-{\bf q^{\prime}}+\frac{{\bf k}+{\bf
k^{\prime}}}{2},\quad{\bf k^{\prime}}\rightarrow{\bf q^{\prime}}-{\bf
q}+\frac{{\bf k}+{\bf k^{\prime}}}{2}$ ${\bf
q}\rightarrow\frac{1}{2}\left({\bf q}+{\bf q^{\prime}}+\frac{{\bf k}-{\bf
k^{\prime}}}{2}\right),\quad{\bf q^{\prime}}\rightarrow\frac{1}{2}\left({\bf
q}+{\bf q^{\prime}}-\frac{{\bf k}-{\bf k^{\prime}}}{2}\right).$

Then $\Gamma_{4}^{(ii)}$ is expressed as

$\displaystyle\Gamma_{4}^{(ii)}({\bf r},{\bf r^{\prime}})=\sum_{{\bf q},{\bf
q^{\prime}}}e^{i({\bf q^{\prime}}-{\bf q})\cdot({\bf r}-{\bf
r^{\prime}})}\langle f(\frac{{\bf r}+{\bf r^{\prime}}}{2},{\bf q})f(\frac{{\bf
r}+{\bf r^{\prime}}}{2},{\bf q^{\prime}})\rangle,$ (44)

which is reminiscent to the asymptotic expression for fluctuations of
intensity Baskov and Chumak (2020). However, for smaller distances, accounting
for partial saturation of fluctuations, it is not factorized to the product of
first moments for PDF. Calculation of (A) is similar to one for term (28).
However since PDFs in (A) depend on the same coordinate $({\bf r}+{\bf
r^{\prime}})/{2}$ both auto-correlations and cross-correlation terms in (A) do
not include space coordinates. After integrations

$\displaystyle\Gamma_{4}^{(ii)}({\bf r},{\bf r}^{\prime})=2\pi
C\int\limits_{0}^{\infty}d\tilde{q}\tilde{q}\left[F_{1}^{\tilde{q}}F_{2}^{\tilde{q}}\big{(}F_{3}^{\tilde{q}}H_{1}^{\tilde{q}}-{G_{1}^{\tilde{q}}}^{\,2}\big{)}\big{(}F_{4}^{\tilde{q}}H_{2}^{\tilde{q}}-{G_{2}^{\tilde{q}}}^{\,2}\big{)}\right]^{-\frac{1}{2}}\times\hskip
184.9429pt$ (45)
$\displaystyle\exp\left\\{-\frac{{\tilde{q}}^{2}}{H_{1}^{\tilde{q}}-{G_{1}^{\tilde{q}}}^{\,2}/F_{3}^{\tilde{q}}}\right\\}\exp\left\\{-\left(\frac{\rho_{\parallel}^{\prime
2}}{F_{1}^{\tilde{q}}}+\frac{\rho_{\perp}^{\prime
2}}{F_{2}^{\tilde{q}}}\right)\right\\}\exp\left\\{-i2{\tilde{q}}\rho_{\parallel}\frac{q_{0}}{z}\right\\},$

where functions $F$,$G$ and $H$ does not depend on ${\bf r}$, ${\bf
r^{\prime}}$ and one should put $2\tilde{q}\tau$ instead of
${\Delta\tilde{r}}$ in corresponding expressions.

## References

  * Capraro  _et al._ (2012) I. Capraro, A. Tomaello, A. Dall’Arche, F. Gerlin, R. Ursin, G. Vallone, and P. Villoresi, Phys. Rev. Lett. 109, 200502 (2012).
  * Usenko  _et al._ (2012) V. C. Usenko, B. Heim, C. Peuntinger, C. Wittmann, C. Marquardt, G. Leuchs, and R. Filip, New Journal of Physics 14, 093048 (2012).
  * Hosseinidehaj and Malaney (2015) N. Hosseinidehaj and R. Malaney, Phys. Rev. A 91, 022304 (2015).
  * Aspelmeyer  _et al._ (2003) M. Aspelmeyer, T. Jennewein, M. Pfennigbauer, W. Leeb, and A. Zeilinger, IEEE Journal of Selected Topics in Quantum Electronics 9, 1541 (2003).
  * Ma  _et al._ (2012) X.-S. Ma _et al._ , Nature 489, 269 (2012).
  * Ren  _et al._ (2017) J.-G. Ren _et al._ , Nature 549, 70 (2017).
  * Hofmann  _et al._ (2019) K. Hofmann, A. A. Semenov, W. Vogel, and M. Bohmann, Physica Scripta 94, 125104 (2019).
  * Ursin  _et al._ (2007) R. Ursin _et al._ , Nature Physics 3, 481 (2007).
  * Yin  _et al._ (2017) J. Yin _et al._ , Science 356, 1140 (2017).
  * Peuntinger  _et al._ (2014) C. Peuntinger, B. Heim, C. R. Müller, C. Gabriel, C. Marquardt, and G. Leuchs, Phys. Rev. Lett. 113, 060502 (2014).
  * Vasylyev  _et al._ (2016) D. Vasylyev, A. A. Semenov, and W. Vogel, Phys. Rev. Lett. 117, 090501 (2016).
  * Meyers (2015) R. E. Meyers, Free-space and atmospheric quantum communications, in  _Advanced Free Space Optics (FSO): A Systems Approach_ (Springer New York, New York, NY, 2015) pp. 343–387.
  * Cao  _et al._ (2020) Y. Cao _et al._ , Phys. Rev. Lett. 125, 260503 (2020).
  * Kish  _et al._ (2020) S. P. Kish, E. Villasenor, R. Malaney, K. A. Mudge, and K. J. Grant, Quantum Engineering 2, e50 (2020).
  * Fante (1975) R. Fante, Proceedings of the IEEE 63, 1669 (1975).
  * Fante (1980) R. Fante, Proceedings of the IEEE 68, 1424 (1980).
  * Tatarskii (1971) V. I. Tatarskii,  _The effects of the turbulent atmosphere on wave propagation_ (Israel Program for Scientific Translations ; Reproduced by National Technical Information Service, U.S. Dept. of Commerce, Jerusalem; Springfield, VA., 1971).
  * 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).
  * DeWolf (1968) D. A. DeWolf, J. Opt. Soc. Am. 58, 461 (1968).
  * Gracheva and Gurvich (1965) M. E. Gracheva and A. S. Gurvich, Soviet Radiophysics 8, 511 (1965).
  * Andrews  _et al._ (1997) L. C. Andrews, R. L. Phillips, and A. R. Weeks, J. Opt. Soc. Am. A 14, 1938 (1997).
  * Dashen (1979) R. Dashen, Journal of Mathematical Physics 20, 894 (1979).
  * Banakh and Mironov (1979) V. A. Banakh and V. L. Mironov, Opt. Lett. 4, 259 (1979).
  * Berman and Chumak (2006) G. P. Berman and A. A. Chumak, Phys. Rev. A 74, 013805 (2006).
  * Berman and Chumak (2009) G. P. Berman and A. A. Chumak, Phys. Rev. A 79, 063848 (2009).
  * Berman  _et al._ (2007) G. P. Berman, A. A. Chumak, and V. N. Gorshkov, Phys. Rev. E 76, 056606 (2007).
  * Chumak and Stolyarov (2013) O. O. Chumak and E. V. Stolyarov, Phys. Rev. A 88, 013855 (2013).
  * Chumak and Stolyarov (2014) O. O. Chumak and E. V. Stolyarov, Phys. Rev. A 90, 063832 (2014).
  * Baskov and Chumak (2018) R. A. Baskov and O. O. Chumak, Phys. Rev. A 97, 043817 (2018).
  * Chumak and Baskov (2016) O. O. Chumak and R. A. Baskov, Phys. Rev. A 93, 033821 (2016).
  * Baskov and Chumak (2020) R. A. Baskov and O. O. Chumak, Journal of Optics 22, 105603 (2020).
  * Vellekoop  _et al._ (2010) I. M. Vellekoop, A. Lagendijk, and A. P. Mosk, Nature Photonics 4, 320 (2010).
  * Newman and Webb (2014) J. A. Newman and K. J. Webb, Phys. Rev. Lett. 113, 263903 (2014).
  * Popoff  _et al._ (2014) S. M. Popoff, A. Goetschy, S. F. Liew, A. D. Stone, and H. Cao, Phys. Rev. Lett. 112, 133903 (2014).
  * Vellekoop and Mosk (2007) I. M. Vellekoop and A. P. Mosk, Opt. Lett. 32, 2309 (2007).
  * Katz  _et al._ (2012) O. Katz, E. Small, and Y. Silberberg, Nature Photonics 6, 549 (2012).
  * Hardy and Shapiro (2013) N. D. Hardy and J. H. Shapiro, Phys. Rev. A 87, 023820 (2013).
  * Zhang  _et al._ (2010) P. Zhang, W. Gong, X. Shen, and S. Han, Phys. Rev. A 82, 033817 (2010).
  * Wang  _et al._ (2010) F. Wang, Y. Cai, and O. Korotkova, in  _Atmospheric and Oceanic Propagation of Electromagnetic Waves IV_, Vol. 7588, edited by O. Korotkova, International Society for Optics and Photonics (SPIE, 2010) pp. 137 – 144\.
  * Shi  _et al._ (2013) D. Shi, C. Fan, P. Zhang, H. Shen, J. Zhang, C. Qiao, and Y. Wang, Opt. Express 21, 2050 (2013).
  * Semenov and Vogel (2009) A. A. Semenov and W. Vogel, Phys. Rev. A 80, 021802 (2009).
  * Xue  _et al._ (2020) Y. Xue, L. Shi, W. Chen, Z. Yin, G.-j. Fan-Yuan, H. Fu, Q. Lu, and J. Wei, Phys. Rev. A 102, 062602 (2020).
  * Tarasenko  _et al._ (1992) A. A. Tarasenko, P. M. Tomchuk, and A. A. Chumak,  _Fluctuations in the bulk and on the surface of solids_ (Naukova Dumka, Kyiv, 1992) (in Russian).
  * Rarenko  _et al._ (1992) A. I. Rarenko, A. A. Tarasenko, , and A. A. Chumak, Ukr. J. Phys. 37, 1577 (1992).
  * Chumak and Sushkova (2012) O. Chumak and N. Sushkova, Ukr. J. Phys. 57, 30 (2012).
  * Strohbehn and Clifford (1967) J. Strohbehn and S. Clifford, IEEE Transactions on Antennas and Propagation 15, 416 (1967).
  * Landau and Lifshitz (1960) L. D. Landau and E. M. Lifshitz,  _Electrodynamics of continuous media vol. 8._ (Pergamon Press, London, 1960).
  * Vasylyev  _et al._ (2019) D. Vasylyev, W. Vogel, and F. Moll, Phys. Rev. A 99, 053830 (2019).
  * Barrios and Dios (2012) R. Barrios and F. Dios, Opt. Express 20, 13055 (2012).
  * Vetelino  _et al._ (2007) F. S. Vetelino, C. Young, L. Andrews, and J. Recolons, Appl. Opt. 46, 2099 (2007).
  * Fried (1967) D. L. Fried, J. Opt. Soc. Am. 57, 169 (1967).
  * Vasylyev  _et al._ (2018) D. Vasylyev, W. Vogel, and A. A. Semenov, Phys. Rev. A 97, 063852 (2018).
  * Churnside (1991) J. H. Churnside, Appl. Opt. 30, 1982 (1991).