# Partially Coherent Electromagnetic Beam Propagation in Random Media

Josselin Garnier Centre de Mathématiques Appliquées, Ecole Polytechnique,
91128 Palaiseau Cedex, France<EMAIL_ADDRESS>Knut Sølna
Department of Mathematics, University of California, Irvine CA 92697
<EMAIL_ADDRESS>
###### Abstract

A theory for the characterization of the fourth moment of electromagnetic wave
beams is presented in the case when the source is partially coherent. A
Gaussian-Schell model is used for the partially coherent random source. The
white-noise paraxial regime is considered, which holds when the wavelength is
much smaller than the correlation radius of the source, the beam radius of the
source, and the correlation length of the medium, which are themselves much
smaller than the propagation distance. The complex wave amplitude field can
then be described by the Itô-Schrödinger equation. This equation gives closed
evolution equations for the wave field moments at all orders and here the
fourth order equations are considered. The general fourth moment equations are
solved explicitly in the scintillation regime (when the correlation radius of
the source is of the same order as the correlation radius of the medium, but
the beam radius is much larger) and the result gives a characterization of the
intensity covariance function. The form of the intensity covariance function
derives from the solution of the transport equation for the Wigner
distribution associated with the second wave moment. The fourth moment results
for polarized waves is used in an application to imaging of partially coherent
sources.

###### Contents

  1. 1 Introduction
  2. 2 Electromagnetic Waves in the White-Noise Paraxial Regime
  3. 3 The Intensity Covariance for Polarized Waves
  4. 4 Estimation of Partially Coherent Electromagnetic Sources
  5. 5 Conclusions
  6. A The Scintillation Regime for the Electromagnetic Waves
  7. B The Mean Polarized-Wigner Transform for a Partially Coherent Beam
  8. C The Mutual Coherence Function in the Strongly Scattering Regime
  9. D The Intensity Covariance Function

##  1 Introduction

We consider beam wave propagation in complex media in the situation when we
model both the source and the medium as being random. Modeling with a random
or complex medium in the context of wave propagation is natural in many
situations. In the early foundational work [32, 33] and also in [1, 26] a main
motivation was propagation through the turbulent atmosphere, but there are
many other important applications as well. In cases where one considers
propagation through the fluctuating ocean, the earth’s crust or through
biological tissue it is also natural and convenient in many cases to model in
terms of a random medium [12, 26, 28]. In these cases the medium may be too
complex to describe pointwise, but one can hope to be able to describe or
model the statistics of the medium fluctuations. The challenge is then to
capture the complicated coupling between the medium fluctuations and the wave
field and to understand how a particular model for the random medium
statistics affects the statistics of the wave field which becomes a random
field due to random scattering. In this paper we consider the case of beam
waves or paraxial waves in a high frequency long range propagation scenario,
when the wavelength is much smaller than the correlation radius of the source,
the beam radius of the source, and the correlation length of the medium, which
are themselves much smaller than the propagation distance. In this case we can
approximate the wave field in terms of the solution of the Itô-Schrödinger
equation. This equation was analyzed in [11] and derived from the Helmholtz
equation in [16]. Despite the long history of the theory of waves in random
media a rigorous and explicit description of the fourth moment was only
obtained in [20] in the scintillation regime (when the correlation radius of
the source is of the same order as the correlation radius of the medium, but
the beam radius is much larger). Here we generalize this fourth moment theory
to the case of polarized waves. This is a deep and quite surprising result in
the theory of waves in random media. It states that even though the polarized
wave has only partially lost its coherence due to scattering it behaves from
the point of view of the fourth moment as if it was a Gaussian field. In some
sense this quasi-Gaussian property explains some of the success or robustness
of the theory of waves in random media since the second-order
characterization, which in general is relatively easy to obtain, also explains
the behavior of fourth-order wave functionals.

An important aspect of our modeling in this paper is that we also consider the
source as being random. Modeling with a random source may be motivated by the
complexity of the source as when one considers emission from a star. A second
motivation for understanding and analyzing beam wave propagation from a random
source is that such sources have been promoted as being desirable for
scintillation reduction when beaming through a complex medium, see [24, 3, 2,
29, 35, 36] and references therein. Scintillation here refers to the situation
that the transmitted beam intensity fluctuates rapidly due to scattering over
the propagation path. One intuition is that by using a complex source one gets
a better mixing over wave ray paths and a scintillation reduction. Our
objective here is to develop an analytic framework where in particular such
questions can be analyzed. We consider the Gaussian-Schell model for the
source when the source coherence statistics is defined in terms of Gaussian
envelopes. This gives rather simple and convenient forms for the wave
statistics, but the theory can easily be modified to more general forms for
the source statistics. The situation with a partially coherent source, but a
homogeneous medium was considered in [14], while we here consider the case
when also the medium is random. Note that we focus on the fourth moment, while
in [13] moments of all orders were considered under some simplifying
assumptions. The assumption that allows us to get explicit expressions for the
fourth moment is to assume the scintillation regime, when the correlation
radius of the source is of the same order as the correlation radius of the
medium, but the beam radius is much larger. Fluctuations of scalar wave
intensity was considered in [27, 25] for a beam type propagation when the
random medium fluctuations are Gaussian, while we assume Gaussianity for the
source, but not for the medium. In [5] a characterization of intensity
fluctuation and how it depends on the regularity of deterministic initial data
is presented. Here we consider the case with random initial data with smooth
Gaussian statistics.

We remark that there are a number of approaches to model high-order moments of
the wave field that are based on perturbative approaches. Indeed, the
derivation of such approximations is based on the premise that the waves is
only perturbed or affected by the scattering to lower order [1, 4, 10]. In
this paper we describe an analytic framework that gives a rigorous scaling
limit identification of the fourth-order moment in the saturated regime when
the incoherent or scattered part of the wave field is of order one and is as
large as or larger than the coherent component. In fact this description also
captures the situation when the wave is fully incoherent and the coherent part
of the wave energy is essentially fully scattered.

We describe next a main result in the paper. We consider the situation when
the time-harmonic electric field in the source plane $z=0$ correspond to a
partially coherent beam and has the form

$\vec{\itbf E}(z=0,{{\itbf x}})=\sum_{j=1}^{2}f_{j}({{\itbf x}})\hat{\itbf
e}_{j},$ (1)

where $\hat{\itbf e}_{1}$ and $\hat{\itbf e}_{2}$ are two orthogonal unit
vectors in the transverse plane, $z$ is the beam propagation direction and
${{\itbf x}}$ the lateral spatial coordinates. The functions $f_{1}$, $f_{2}$
are zero-mean Gaussian processes with covariance

$\mathbb{E}\big{[}f_{j}\big{(}{{\itbf r}}+\frac{{{\itbf
q}}}{2}\big{)}\overline{f_{l}}\big{(}{{\itbf r}}-\frac{{{\itbf
q}}}{2}\big{)}\big{]}=\left\\{\begin{array}[]{ll}A_{j}^{2}\exp\Big{(}-\frac{|{{\itbf
r}}|^{2}}{r_{o}^{2}}-\frac{|{{\itbf q}}|^{2}}{4\rho_{o}^{2}}\Big{)}&\mbox{ if
}j=l,\\\ A_{j}A_{l}\chi\exp\Big{(}-\frac{|{{\itbf
r}}|^{2}}{r_{o}^{2}}-\frac{|{{\itbf q}}|^{2}}{4\rho_{1}^{2}}\Big{)}&\mbox{ if
}j\neq l,\end{array}\right.$ (2)

where the polarization degree $\chi\in[-1,1]$, $r_{o}$ is the beam radius of
the source, $\rho_{o}$ is the correlation radius of the source and we have
$\rho_{1}\leq\rho_{o}$. We refer to [9] for further background on modeling
with polarized waves. Consider first the second field moment in the form of
the mutual coherence function in the white-noise paraxial regime (with the
wavelength smaller than the correlation radii of the medium and the beam which
are moreover on the scale of the the beam radius, which in turn is smaller
than the propagation distance, see Appendix A). As discussed in Section 2 the
mutual coherence function of the wave field is in this regime given by

$\mu_{2,jl}(z,{{\itbf r}},{{\itbf q}})=\mathbb{E}\big{[}u_{j}\big{(}z,{{\itbf
r}}+\frac{{{\itbf q}}}{2}\big{)}\overline{u_{l}}\big{(}z,{{\itbf
r}}-\frac{{{\itbf
q}}}{2}\big{)}\big{]}=\left\\{\begin{array}[]{ll}A_{j}^{2}{\cal
H}_{\rho_{o}}(z,{{\itbf r}},{{\itbf q}})&\mbox{ if }j=l,\\\
A_{j}A_{l}\chi{\cal H}_{\rho_{1}}(z,{{\itbf r}},{{\itbf q}})&\mbox{ if }j\neq
l,\end{array}\right.$ (3)

where the fundamental second-order lateral scattering function is defined by

$\displaystyle{\cal H}_{\rho_{o}}(z,{{\itbf r}},{{\itbf q}})=$
$\displaystyle\frac{r_{o}^{2}}{4\pi}\int_{\mathbb{R}^{2}}\exp\Big{(}i{\boldsymbol{\zeta}}\cdot{{\itbf
r}}-\frac{r_{o}^{2}|{\boldsymbol{\zeta}}|^{2}}{4}-\frac{|{\boldsymbol{\zeta}}\frac{z}{k_{o}}-{{\itbf
q}}|^{2}}{4\rho_{o}^{2}}\Big{)}$
$\displaystyle\times\exp\Big{(}\frac{k_{o}^{2}}{4}\int_{0}^{z}C\big{(}\frac{{\boldsymbol{\zeta}}z^{\prime}}{k_{o}}-{{\itbf
q}}\big{)}-C({\bf 0})dz^{\prime}\Big{)}d{\boldsymbol{\zeta}},$ (4)

with $k_{o}$ the central wavenumber and $C$ the covariance function of the
medium fluctuations when integrated in the $z$-dimension (see Eq. (18)). The
intensity is defined by

$I(z,{{\itbf r}})=\sum_{j=1}^{2}|u_{j}(z,{{\itbf r}})|^{2},$ (5)

and the mean intensity is then

$\mathbb{E}\big{[}I({z},{{{\itbf r}}})\big{]}=(A_{1}^{2}+A_{2}^{2}){\cal
H}_{\rho_{o}}(z,{{\itbf r}},{\bf 0}).$ (6)

Consider next the fourth field moment in form of the covariance of the
intensity. As discussed in Section 3, in the scintillation regime where the
wavelength is smaller than the correlation radii of the medium and the beam,
which are smaller than the beam radius, which is smaller than the propagation
distance, the intensity covariance has the form

$\displaystyle{\rm Cov}\big{(}I\big{(}{z},{{{\itbf r}}}+\frac{{{\itbf
q}}}{2}\big{)},I\big{(}{z},{{{\itbf r}}}-\frac{{{\itbf q}}}{2}\big{)}\big{)}$
$\displaystyle=$ $\displaystyle\big{(}A_{1}^{4}+A_{2}^{4}\big{)}|{\cal
H}_{\rho_{o}}(z,{{\itbf r}},{{\itbf q}})|^{2}$ (7)
$\displaystyle\hbox{}+2A_{1}^{2}A_{2}^{2}\chi^{2}|{\cal
H}_{\rho_{1}}(z,{{\itbf r}},{{\itbf q}})|^{2}.$

This representation proves the quasi-Gaussian property, that the fourth moment
(intensity covariance) derives from the second moment (the mutual coherence
function), as is the case in general for Gaussian random fields. From this
description we can identify the intensity decoherence scale and spreading
scale and we discuss this explicitly in the strong scattering case in Appendix
C. More generally the quasi-Gaussian property provides the basis for analysis
of a number of wave propagation challenges and we discuss one such application
in Section 4. We remark that fourth moments of the wave fields that are more
general than the intensity covariance can be obtained via Proposition D.1 in
Appendix D and that the quasi-Gaussian property holds for such general
moments.

The outline of the paper is as follows. In Section 2 we describe modeling of
the partially coherent source in the polarized case and we present the Itô-
Schrödinger equation describing wave propagation in the white-noise paraxial
regime. In Section 3 we present the fourth moment result for the polarized
waves. The main result that allows us to do this is that the Itô-Schrödinger
equations for the polarization modes are dynamically uncoupled, however,
statistically coupled. In Section 4 we discuss an application of the theory
that we have developed to source imaging (using the intensity covariance
function instead of the mutual coherence function as in [6]). The details of
the derivation and also a discussion of the second-order moment or mutual
coherence function can be found in the appendices.

##  2 Electromagnetic Waves in the White-Noise Paraxial Regime

We consider propagation of a partially coherent electromagnetic beam waves
through a three-dimensional random medium. Maxwell’s equations for the three-
dimensional electric field $\vec{\itbf E}$ and the three-dimensional magnetic
field strength $\vec{\itbf H}$ are in the time-harmonic case (with frequency
$\omega_{o}$):

$\displaystyle\nabla\times\vec{\itbf E}=-i\omega_{o}\mu(z,{{\itbf
x}})\vec{\itbf H},$ (8) $\displaystyle\nabla\cdot({\epsilon}(z,{{\itbf
x}})\vec{\itbf E})=\rho(z,{{\itbf x}}),$ (9)
$\displaystyle\nabla\times\vec{\itbf H}=\vec{\itbf J}^{(s)}(z,{{\itbf
x}})+i\omega_{o}{\epsilon}(z,{{\itbf x}})\vec{\itbf E},$ (10)
$\displaystyle\nabla\cdot(\mu(z,{{\itbf x}})\vec{\itbf H})=0.$ (11)

The term $\vec{\itbf J}^{(s)}$ is a current source term, ${\epsilon}$ is the
dielectric permittivity of the medium, and $\mu$ is the magnetic permeability
of the medium. Note that the equation of continuity of charge
$i\omega_{o}\rho+\nabla\cdot\vec{\itbf J}^{(s)}=0$ is automatically satisfied.

We assume that

  * •

The medium is randomly heterogeneous:

$\displaystyle{\epsilon}(z,{{\itbf x}})$
$\displaystyle={\epsilon}_{o}\big{[}1+m_{\epsilon}(z,{{\itbf x}})\big{]},$
(12) $\displaystyle\mu(z,{{\itbf x}})$
$\displaystyle=\mu_{o}\big{[}1+m_{\mu}(z,{{\itbf x}})\big{]}.$ (13)

The random processes $m_{\epsilon}(z,{{\itbf x}})$ and $m_{\mu}(z,{{\itbf
x}})$ are bounded, stationary, and zero-mean and they satisfy ergodic (mixing)
conditions in $z$.

  * •

We consider a partially coherent source $\vec{\itbf f}({{\itbf x}})$, which is
a field with Gaussian statistics and mean zero that is localized in the plane
$z=0$. We address the case of a Gauss-Schell model for the source. The source
is then $\vec{\itbf J}^{(s)}(z,{{\itbf
x}})=2\mu_{o}^{-1/2}{\epsilon}_{0}^{1/2}\vec{\itbf f}({{\itbf x}})\delta(z)$,
where $f_{3}=0$ and $f_{1}$, $f_{2}$ are zero-mean Gaussian processes with
covariance

$\mathbb{E}\big{[}f_{j}\big{(}{{\itbf r}}+\frac{{{\itbf
q}}}{2}\big{)}\overline{f_{l}}\big{(}{{\itbf r}}-\frac{{{\itbf
q}}}{2}\big{)}\big{]}=\left\\{\begin{array}[]{ll}A_{j}^{2}\exp\Big{(}-\frac{|{{\itbf
r}}|^{2}}{r_{o}^{2}}-\frac{|{{\itbf q}}|^{2}}{4\rho_{o}^{2}}\Big{)}&\mbox{ if
}j=l,\\\ A_{j}A_{l}\chi\exp\Big{(}-\frac{|{{\itbf
r}}|^{2}}{r_{o}^{2}}-\frac{|{{\itbf q}}|^{2}}{4\rho_{1}^{2}}\Big{)}&\mbox{ if
}j\neq l.\end{array}\right.$ (14)

All parameters are real and positive, with $\chi\in[-1,1]$. The parameters
have to satisfy several constraints to ensure that we deal with a well-defined
covariance function:

$\rho_{o}\leq
r_{o},\quad\quad\chi^{2}\rho_{o}^{-2}+(1-\chi^{2})r_{o}^{-2}\leq\rho_{1}^{-2}\leq\rho_{o}^{-2}.$

Here $r_{o}$ is the radius of the beam and $\rho_{o}$ is the correlation
radius. The special case $\rho_{o}=r_{o},\chi=0$ corresponds to a coherent
source with components $f_{j}({{\itbf x}})=f_{o,j}\exp(-|{{\itbf
x}}|^{2}/(2r_{o}^{2}))$ where $f_{o,j},j=1,2$ are two independent complex-
valued zero-mean Gaussian random variables with variance $A_{j}^{2}$. In the
general case with $\rho_{o}<r_{o}$ the field is partially coherent and the
field components have the form of correlated speckle patterns with speckle
spots with typical radius $\rho_{o}$ and with an overall intensity envelope
that is $\exp(-|{{\itbf r}}|^{2}/r_{o}^{2})$.

In the white-noise paraxial regime (which holds when the wavelength is much
smaller than the correlation radius of the source, the correlation radius of
the medium, and the beam radius, which are themselves much smaller than the
propagation distance) the electric field modulo a range-dependent phase (see
Appendix A) has the form

$\vec{\itbf E}(z,{{\itbf x}})=\sum_{j=1}^{2}u_{j}(z,{{\itbf x}})\hat{\itbf
e}_{j},$ (15)

where $\hat{\itbf e}_{1}$ and $\hat{\itbf e}_{2}$ are the unit vectors in the
transverse plane pointing in the $x$ and $y$ directions and the complex
amplitude fields $u_{j}$ are the solution of the following statistically
coupled Itô-Schrödinger equations [16, 18, 19]:

$d{u}_{j}(z,{{\itbf x}})=\frac{i}{2k_{o}}\Delta_{{{\itbf x}}}{u}_{j}(z,{{\itbf
x}})dz+\frac{ik_{o}}{2}{u}_{j}(z,{{\itbf x}})\circ d{B}(z,{{\itbf x}}),\quad
j=1,2,$ (16)

with the initial condition in the plane $z=0$:

${u}_{j}(z=0,{{\itbf x}})=f_{j}({{\itbf x}}).$

Here the random process $B(z,{{\itbf x}})$ is a real-valued Brownian field
with a covariance that derives from two-point statistics in the model for the
medium fluctuations in (12-13)

$\mathbb{E}[{B}(z,{{\itbf x}}){B}(z^{\prime},{{\itbf
x}}^{\prime})]={\min\\{z,z^{\prime}\\}}{C}({{\itbf x}}-{{\itbf x}}^{\prime}),$
(17)

where

$\displaystyle C({{\itbf
x}})=\int_{\mathbb{R}}\mathbb{E}[(m_{\epsilon}+m_{\mu})({{\itbf
x}}^{\prime}+{{\itbf x}},z^{\prime}+z)(m_{\epsilon}+m_{\mu})({{\itbf
x}}^{\prime},z^{\prime})]dz.$ (18)

Note therefore that the evolution equations for the lateral components of the
electromagnetic field is driven by the same Brownian field. We remark that
$C({\bf 0})$ can be interpreted as the product of the variance of the
fluctuations of the random medium times its longitudinal correlation radius:

$C({\bf 0})=\sigma^{2}\ell_{z},$

for $\sigma$ the standard deviation of the medium fluctuations:
$\sigma^{2}=\mathbb{E}[(m_{\epsilon}+m_{\mu})^{2}({{\itbf
x}}^{\prime},z^{\prime})]$.

The derivation of (16) from the random three-dimensional scalar wave equation
is presented in [16]. Its derivation for the Maxwell’s equations (8-11) is
presented in [18, 19]. In Eq. (16) the symbol $\circ$ stands for the
Stratonovich stochastic integral [16]. The first- and second-order moments of
the wave field solution of (16) have been studied in [16, 18] and recover
results derived in [26]. In view of the centered partial coherent source we
find that the first-order moment of the wave field is zero. The governing
equations for the higher order moments can be identified via Itô calculus for
Hilbert space valued processes. We find in particular that the second-order
moment of the wave field (mutual coherence function) defined by

$\mu_{2,jl}(z,{{\itbf r}},{{\itbf q}})=\mathbb{E}\Big{[}u_{j}\big{(}z,{{\itbf
r}}+\frac{{{\itbf q}}}{2}\big{)}\overline{u_{l}\big{(}z,{{\itbf
r}}-\frac{{{\itbf q}}}{2}\big{)}}\Big{]}$ (19)

satisfies [21]

$\frac{\partial\mu_{2,jl}}{\partial z}=\frac{i}{k_{o}}\nabla_{{{\itbf
r}}}\cdot\nabla_{{{\itbf
q}}}\mu_{2,jl}+\frac{k_{o}^{2}}{4}{U}_{2}\big{(}{{\itbf q}}\big{)}\mu_{2,jl},$
(20)

with the potential $U_{2}({{\itbf q}})=C({{\itbf q}})-C({\bf 0})$ and the
initial condition determined by (14). As shown in Appendix B, it then follows
that the mutual coherence function is given by

$\mu_{2,jl}(z,{{\itbf r}},{{\itbf
q}})=\left\\{\begin{array}[]{ll}A_{j}^{2}{\cal H}_{\rho_{o}}(z,{{\itbf
r}},{{\itbf q}})&\mbox{ if }j=l,\\\ A_{j}A_{l}\chi{\cal
H}_{\rho_{1}}(z,{{\itbf r}},{{\itbf q}})&\mbox{ if }j\neq
l,\end{array}\right.$ (21)

where ${\cal H}_{\rho_{o}}(z,{{\itbf r}},{{\itbf q}})$ is defined by (4).

We next show how our main quantity of interest - the intensity covariance -
can be expressed in terms of the mutual coherence function.

##  3 The Intensity Covariance for Polarized Waves

The intensity is defined by (5) and by (3), its mean is given by (6). The
second-order moment of the intensity is

$\mathbb{E}\big{[}I\big{(}z,{{\itbf x}}_{1}\big{)}I\big{(}z,{{\itbf
x}}_{2}\big{)}\big{]}=\sum_{j,l=1}^{2}\mu_{4,jl}(z,{{\itbf x}}_{1},{{\itbf
x}}_{2},{{\itbf x}}_{1},{{\itbf x}}_{2}),$ (22)

where the $\mu_{4,jl}$’s are defined by

$\mu_{4,jl}(z,{{\itbf x}}_{1},{{\itbf x}}_{2},{{\itbf y}}_{1},{{\itbf
y}}_{2})=\mathbb{E}\big{[}u_{j}(z,{{\itbf x}}_{1})\overline{u_{j}(z,{{\itbf
y}}_{1})}u_{l}(z,{{\itbf x}}_{2})\overline{u_{l}(z,{{\itbf y}}_{2})}\big{]}.$

These moments satisfy the equation (48) and have the initial conditions:

$\displaystyle\mu_{4,jj}(z=0,{{\itbf x}}_{1},{{\itbf x}}_{2},{{\itbf
y}}_{1},{{\itbf y}}_{2})$ $\displaystyle=A_{j}^{4}\exp\Big{(}-\frac{|{{\itbf
x}}_{1}+{{\itbf y}}_{1}|^{2}}{4r_{o}^{2}}-\frac{|{{\itbf x}}_{1}-{{\itbf
y}}_{1}|^{2}}{4\rho_{o}^{2}}-\frac{|{{\itbf x}}_{2}+{{\itbf
y}}_{2}|^{2}}{4r_{o}^{2}}-\frac{|{{\itbf x}}_{2}-{{\itbf
y}}_{2}|^{2}}{4\rho_{o}^{2}}\Big{)}$
$\displaystyle\quad+A_{j}^{4}\exp\Big{(}-\frac{|{{\itbf x}}_{1}+{{\itbf
y}}_{2}|^{2}}{4r_{o}^{2}}-\frac{|{{\itbf x}}_{1}-{{\itbf
y}}_{2}|^{2}}{4\rho_{o}^{2}}-\frac{|{{\itbf x}}_{2}+{{\itbf
y}}_{1}|^{2}}{4r_{o}^{2}}-\frac{|{{\itbf x}}_{2}-{{\itbf
y}}_{1}|^{2}}{4\rho_{o}^{2}}\Big{)},$

for $j=1,2$ and

$\displaystyle\mu_{4,jl}(z=0,{{\itbf x}}_{1},{{\itbf x}}_{2},{{\itbf
y}}_{1},{{\itbf y}}_{2})$
$\displaystyle=A_{j}^{2}A_{l}^{2}\exp\Big{(}-\frac{|{{\itbf x}}_{1}+{{\itbf
y}}_{1}|^{2}}{4r_{o}^{2}}-\frac{|{{\itbf x}}_{1}-{{\itbf
y}}_{1}|^{2}}{4\rho_{o}^{2}}-\frac{|{{\itbf x}}_{2}+{{\itbf
y}}_{2}|^{2}}{4r_{o}^{2}}-\frac{|{{\itbf x}}_{2}-{{\itbf
y}}_{2}|^{2}}{4\rho_{o}^{2}}\Big{)}$
$\displaystyle\quad+A_{j}^{2}A_{l}^{2}\chi^{2}\exp\Big{(}-\frac{|{{\itbf
x}}_{1}+{{\itbf y}}_{2}|^{2}}{4r_{o}^{2}}-\frac{|{{\itbf x}}_{1}-{{\itbf
y}}_{2}|^{2}}{4\rho_{1}^{2}}-\frac{|{{\itbf x}}_{2}+{{\itbf
y}}_{1}|^{2}}{4r_{o}^{2}}-\frac{|{{\itbf x}}_{2}-{{\itbf
y}}_{1}|^{2}}{4\rho_{1}^{2}}\Big{)},$

for $j\neq l$.

Consequently, as shown in Appendix D, in the scintillation regime (which holds
in the white-noise paraxial regime when additionally the correlation radius of
the source is of the same order as the correlation radius of the medium, but
the beam radius is much larger), the second-order moment of the intensity has
the form

$\displaystyle\mathbb{E}\big{[}I\big{(}z,{{\itbf r}}+\frac{{{\itbf
q}}}{2}\big{)}I\big{(}z,{{\itbf r}}-\frac{{{\itbf q}}}{2}\big{)}\big{]}=$
$\displaystyle\big{(}A_{1}^{4}+A_{2}^{4}\big{)}\big{[}{\cal
I}_{\rho_{o}}(2{{\itbf r}},{\bf 0})+{\cal J}_{\rho_{o}}(2{{\itbf r}},{{\itbf
q}})\big{]}$ $\displaystyle+2A_{1}^{2}A_{2}^{2}\big{[}{\cal
I}_{\rho_{o}}(2{{\itbf r}},{\bf 0})+\chi^{2}{\cal J}_{\rho_{1}}(2{{\itbf
r}},{{\itbf q}})\big{]},$ (23)

where ${\cal I}_{\rho_{o}}$ and ${\cal J}_{\rho_{o}}$ are defined by (67) and
(68). This can also be written as

$\displaystyle\mathbb{E}\big{[}I\big{(}z,{{\itbf r}}+\frac{{{\itbf
q}}}{2}\big{)}I\big{(}z,{{\itbf r}}-\frac{{{\itbf q}}}{2}\big{)}\big{]}=$
$\displaystyle\big{(}A_{1}^{4}+A_{2}^{4}\big{)}\big{[}{\cal
H}_{\rho_{o}}(z,{{\itbf r}},{\bf 0})^{2}+|{\cal H}_{\rho_{o}}(z,{{\itbf
r}},{{\itbf q}})|^{2}\big{]}$ $\displaystyle+2A_{1}^{2}A_{2}^{2}\big{[}{\cal
H}_{\rho_{o}}(z,{{\itbf r}},{\bf 0})^{2}+\chi^{2}|{\cal
H}_{\rho_{1}}(z,{{\itbf r}},{{\itbf q}})|^{2}\big{]},$ (24)

or equivalently (7). This shows that the field satisfies the Gaussian rule for
the fourth-order moment in the scintillation regime:

$\displaystyle\mathbb{E}\big{[}I\big{(}z,{{\itbf r}}+\frac{{{\itbf
q}}}{2}\big{)}I\big{(}z,{{\itbf r}}-\frac{{{\itbf q}}}{2}\big{)}\big{]}$
$\displaystyle=\sum_{j,l=1}^{2}\mathbb{E}\big{[}u_{j}\big{(}z,{{\itbf
r}}+\frac{{{\itbf q}}}{2}\big{)}\overline{u_{j}}\big{(}z,{{\itbf
r}}+\frac{{{\itbf q}}}{2}\big{)}\big{]}\mathbb{E}\big{[}u_{l}\big{(}z,{{\itbf
r}}-\frac{{{\itbf q}}}{2}\big{)}\overline{u_{l}}\big{(}z,{{\itbf
r}}-\frac{{{\itbf q}}}{2}\big{)}\big{]}$
$\displaystyle\quad+\sum_{j,l=1}^{2}\mathbb{E}\big{[}u_{j}\big{(}z,{{\itbf
r}}+\frac{{{\itbf q}}}{2}\big{)}\overline{u_{l}}\big{(}z,{{\itbf
r}}-\frac{{{\itbf q}}}{2}\big{)}\big{]}\mathbb{E}\big{[}u_{l}\big{(}z,{{\itbf
r}}-\frac{{{\itbf q}}}{2}\big{)}\overline{u_{j}}\big{(}z,{{\itbf
r}}+\frac{{{\itbf q}}}{2}\big{)}\big{]}.$ (25)

The quasi-Gaussianity property and intensity covariance expression that we
have just identified are useful in various applications. For instance in
applications to imaging based on wave field coherence the signal-to-noise
ratio will in general depend on fourth-order wave field moments [7, 15]. These
moments can then be related to the mutual coherence function for the wave
field via the theory presented in this paper. We discuss next an application
to the inverse source problem for partially coherent beam sources in complex
media using information that can be extracted from the mean intensity
covariance.

##  4 Estimation of Partially Coherent Electromagnetic Sources

We consider the problem of characterizing a partially coherent source. The
goal is to determine the parameters of a Gauss-Schell source from the
measurements of the intensity after a propagation distance $z$. Note that here
we base the estimate on measurements of intensity only while in for instance
in [6] the estimate is based on coherence of the measurements of the field
itself. In [31] the authors use a cross phase in the source for stable
transmission of the coherence pattern in the source field. This then allows
for transmission of information. In order to relate the field coherence to
coherence in the intensity the authors use a Gaussian approximation for the
field. It follows from the analysis in our paper that such an approximations
is valid in the scintillation regime. We remark that in the case with active
multifrequency imaging one can obtain partial phase information from an
appropriate illumination stategy and use of a polarization identity [30]. Here
we consider passive single-frequency intensity based imaging and no phase
information is available.

In the source estimation context considered here we assume the scintillation
regime, moreover, we use the strongly scattering approximation studied in
Appendix C (which holds when $k_{o}^{2}C({\bf 0})z\gg 1$ and $C$ is smooth).
We then find

$\mathbb{E}\big{[}I(z,{{\itbf
r}})\big{]}=(A_{1}^{2}+A_{2}^{2})\frac{r_{o}^{2}}{R^{2}(z;\rho_{o})}\exp\Big{(}-\frac{|{{\itbf
r}}|^{2}}{R^{2}(z;\rho_{o})}\Big{)},$ (26)

where the beam radius is

$R^{2}(z;\rho_{o})=r_{o}^{2}\Big{(}1+\frac{\bar{\gamma}_{2}z^{3}}{3r_{o}^{2}}+\frac{z^{2}}{k_{o}^{2}r_{o}^{2}\rho_{o}^{2}}\Big{)},$
(27)

with $\bar{\gamma}_{2}=-\Delta_{{\itbf x}}C({\bf 0})/4$ being a parameter that
governs the strength of random lateral scattering. In this expression the
original beam radius is $r_{o}$, the third term in the right-hand side gives
the spreading due to diffraction ${\cal O}(z)$ in the homogeneous medium case.
The second term gives the anomalous spreading ${\cal O}(z^{3/2})$ due to
random scattering. This expression for the beam spreading was also derived in
[23] using the Huygens-Fresnel principle.

The covariance intensity function is

$\displaystyle{\rm Cov}\big{(}I\big{(}z,{{\itbf r}}+\frac{{{\itbf
q}}}{2}\big{)},I\big{(}z,{{\itbf r}}-\frac{{{\itbf q}}}{2}\big{)}\big{)}=$
$\displaystyle\big{(}A_{1}^{4}+A_{2}^{4}\big{)}\frac{r_{o}^{4}}{R(z;\rho_{o})^{4}}\exp\Big{(}-\frac{2|{{\itbf
r}}|^{2}}{R^{2}(z;\rho_{o})}-\frac{|{{\itbf
q}}|^{2}}{2\rho^{2}(z;\rho_{o})}\Big{)}$
$\displaystyle+2A_{1}^{2}A_{2}^{2}\chi^{2}\frac{r_{o}^{4}}{R(z;\rho_{1})^{4}}\exp\Big{(}-\frac{2|{{\itbf
r}}|^{2}}{R^{2}(z;\rho_{1})}-\frac{|{{\itbf
q}}|^{2}}{2\rho^{2}(z;\rho_{1})}\Big{)},$ (28)

where the correlation radius of the beam is

$\rho^{2}(z;\rho_{o})=\rho_{o}^{2}\frac{1+\frac{\bar{\gamma}_{2}z^{3}}{3r_{o}^{2}}+\frac{z^{2}}{k_{o}^{2}\rho_{o}^{2}r_{o}^{2}}}{1+k_{o}^{2}\rho_{o}^{2}\bar{\gamma}_{2}z\big{(}1+\frac{\bar{\gamma}_{2}z^{3}}{12r_{o}^{2}}+\frac{z^{2}}{3k_{o}^{2}\rho_{o}^{2}r_{o}^{2}}\big{)}}.$
(29)

In the homogeneous medium case we have
$\rho(z;\rho_{o})/R(z;\rho_{o})=\rho_{o}/r_{o}$ and the lateral correlation
radius increases with the propagation distance while in the random case
$\rho(z;\rho_{o})={\cal O}(z^{-1/2})$ and the lateral correlation radius
decreases with the propagation distance due to random scattering.

The structure of the intensity covariance in lateral and range coordinates can
be used for imaging of the partially coherent source. Here we consider
measurements at one range $z$ only. Note that the expressions in (26) and (28)
are computed based on averaging with respect to the statistics of both the
source and of the medium. We assume that these means can be identified with a
high signal-to-noise ratio. This is the case if both the partially coherent
source and the medium fluctuate in time and we average the measurements over a
time interval that is long compared to the turnover times of the medium and of
the source. In cases when the averaging is not efficient or the medium is
stationary (time-independent) it may be necessary with some form of filtering
to enhance statistical stability [8].

In the statistically stable case with long time averaging at the detector the
observation of the mean intensity and the intensity covariance function makes
it possible to extract the beam radii $R(z;\rho_{o})$, $R(z;\rho_{1})$,
correlation radii $\rho(z;\rho_{o})$, $\rho(z;\rho_{1})$ and the intensity
amplitude $(A_{1}^{2}+A_{2}^{2})r_{o}^{2}$. If $\rho_{o}\neq\rho_{1}$, then we
can also extract $(A_{1}^{4}+A_{2}^{4})r_{o}^{4}$ and
$2A_{1}^{2}A_{2}^{2}\chi^{2}r_{o}^{4}$. If $\rho_{o}=\rho_{1}$, then we can
only extract the sum
$(A_{1}^{4}+A_{2}^{4}+2A_{1}^{2}A_{2}^{2}\chi^{2})r_{o}^{4}$.

Given the values of $z$ and $\bar{\gamma}_{2}$, we can then estimate the beam
radius $r_{o}$ of the source, the correlation radii $\rho_{o}$ and $\rho_{1}$,
the total intensity $A_{1}^{2}+A_{2}^{2}$. If $\rho_{o}\neq\rho_{1}$ we can
also extract the polarization degree $\chi$ and the ellipticity
$e^{2}=(A_{2}^{2}-A_{1}^{2})^{2}/(A_{1}^{2}+A_{2}^{2})^{2}$. Indeed, if we
introduce the three following quantities that can be extracted from data
$Y_{1}=(a_{1}^{2}+a_{2}^{2})$ and $Y_{2}=a_{1}^{4}+a_{2}^{4}$,
$Y_{3}=2a_{1}^{2}a_{2}^{2}\chi^{2}$ with $a_{j}=A_{j}r_{o}$, we have

$e^{2}=\frac{(a_{1}^{2}-a_{2}^{2})^{2}}{(a_{1}^{2}+a_{2}^{2})^{2}}=2\frac{Y_{2}}{Y_{1}^{2}}-1,\quad\quad\chi^{2}=\frac{Y_{3}}{Y_{1}^{2}-Y_{2}}.$
(30)

Then, with an estimate of $r_{o}$ we can identify $A_{1}$ and $A_{2}$.

The estimation is possible and reliable provided the propagation distance is
not too large, i.e., $z$ should not be much larger than

$\displaystyle\max\big{(}(r_{o}^{2}/\bar{\gamma}_{2})^{1/3},1/(\bar{\gamma}_{2}k_{o}^{2}\rho_{o}^{2}),1/\bar{\gamma}_{2}\big{)},$

because the statistics of the beam (beam radius and correlation radius) then
becomes essentially independent of the initial values $r_{o}$ and $\rho_{o}$
as shown in Appendix C. We remark also that when $\bar{\gamma}_{2}z\gg 1$ the
paraxial approximation is not valid anymore, as shown in Appendix C, so it is
not possible to use this approach. Note that that in this regime one can then
use a least-squares misfit approach [6] to estimate $r_{o},\rho_{o}$ and
$\rho_{1}$ from $R(z;\rho_{o})$, $R(z;\rho_{1})$, $\rho(z;\rho_{o})$ and
$\rho(z;\rho_{1})$.

##  5 Conclusions

We have considered time-harmonic electromagnetic wave propagation from
partially coherent sources in random media. In many applications of waves it
is of interest to evaluate the fourth moment of the wave field. Here we
present a theory that allows us to describe such moments and we focus on the
specific fourth moments corresponding to the intensity covariance. We present
here such a description for polarized waves. The results follow from the Itô-
Schrödinger equation for the wave field valid in the white-noise paraxial
regime. An important aspect of these Itô-Schrödinger equations is that the
equations describing the evolutions of the transversely polarized modes are
driven by the same Brownian motion, however, such that they are dynamically
uncoupled. The explicit expressions for the fourth moments are derived in a
subsequent scaling regime that we denote the scintillation regime. An
important aspect of the fourth moment analysis is the proof of the quasi-
Gaussian property which means that the fourth moments can be obtained from the
second wave field moments as if the field had Gaussian statistics. We note
that this property holds true even if the wave field is partially coherent. We
also describe an application to the inverse source problem using information
extracted from the observed intensity covariance. We moreover give explicit
expressions for the decorrelation and spreading scales deriving from the
mutual coherence function for probing through strong clutter, these scales
characterize the statistical structure of the wave field in view of the quasi-
Gaussian property (up to fourth order).

##  Appendix A The Scintillation Regime for the Electromagnetic Waves

We discuss here the white-noise paraxial scaling regime that leads to the Itô-
Schrödinger equation in (16). We refer to [18] for the full derivation. The
electromagnetic wave equations have the form (8-11) with the dielectric
permittivity $\epsilon$ and the magnetic permeability $\mu$ of the medium
modeled by (12-13). We denote $\vec{\itbf E}=({E}_{j})_{j=1,2,3}$ and
$\vec{\itbf H}=({H}_{j})_{j=1,2,3}$. The four-dimensional vector
$({E}_{1},{H}_{2},{E}_{2},{H}_{1})$ then satisfies a closed system as shown in
[18]. Let $c_{o}=\mu_{o}^{-1/2}\epsilon_{o}^{-1/2}$ and
$\zeta_{o}=\mu_{o}^{1/2}\epsilon_{o}^{-1/2}$ be the homogeneous propagation
speed and impedance, then we introduce the decomposition

$\displaystyle{E}_{1}(z,{{\itbf x}})$ $\displaystyle=$
$\displaystyle\zeta_{o}^{\frac{1}{2}}\big{(}{a}_{1}(z,{{\itbf
x}})e^{i\frac{\omega z}{c_{o}}}+{b}_{1}(z,{{\itbf x}})e^{-i\frac{\omega
z}{c_{o}}}\big{)},$ $\displaystyle{H}_{2}(z,{{\itbf x}})$ $\displaystyle=$
$\displaystyle\zeta_{o}^{-\frac{1}{2}}\big{(}{a}_{1}(z,{{\itbf
x}})e^{i\frac{\omega z}{c_{o}}}-{b}_{1}(z,{{\itbf x}})e^{-i\frac{\omega
z}{c_{o}}}\big{)},$ $\displaystyle{E}_{2}(z,{{\itbf x}})$ $\displaystyle=$
$\displaystyle\zeta_{o}^{\frac{1}{2}}\big{(}{a}_{2}(z,{{\itbf
x}})e^{i\frac{\omega z}{c_{o}}}+{b}_{2}(z,{{\itbf x}})e^{-i\frac{\omega
z}{c_{o}}}\big{)},$ $\displaystyle{H}_{1}(z,{{\itbf x}})$ $\displaystyle=$
$\displaystyle\zeta_{o}^{-\frac{1}{2}}\big{(}-{a}_{2}(z,{{\itbf
x}})e^{i\frac{\omega z}{c_{o}}}+{b}_{2}(z,{{\itbf x}})e^{-i\frac{\omega
z}{c_{o}}}\big{)}.$

Here, $a_{j},b_{j},j=1,2$ are coefficients of locally forward and backward (in
$z$) propagating plane waves. In the case of a homogeneous medium with
$m_{\epsilon}\equiv 0,m_{\mu}\equiv 0$ this gives an exact decomposition into
forward and backward plane waves with constant coefficients, while with random
medium fluctuations the coefficients $a_{j}$ and $b_{j}$ satisfy coupled
equations.

Let $\sigma$ be the standard deviation of the fluctuations of the medium.
Moreover, assume that the random fluctuations in the index of refraction are
isotropic and denote by $\ell_{\rm c}$ the correlation length of the
fluctuations, $\lambda_{o}$ the carrier wavelength (equal to $2\pi/k_{o}$,
$k_{o}=\omega_{o}/c_{o}$), $L$ the typical propagation distance, $\rho_{o}$
the correlation radius of the source, and $r_{o}$ the radius of the initial
transverse beam-source. In this framework the variance $C({\bf 0})$ of the
Brownian field in the Itô-Schrödinger equation (16) is of order
$\sigma^{2}\ell_{\rm c}$ and the transverse scale of variation of the
covariance function $C({{\itbf x}})$ in (18) is of order $\ell_{\rm c}$.

We next discuss the scintillation regime in more detail. First, we consider
the primary scaling (white-noise paraxial regime) that leads to the Itô-
Schrödinger equation (16), when the propagation distance is much larger than
the correlation length of the medium, the correlation radius of the source and
the beam radius, which are themselves much larger than the wavelength,
moreover, the medium fluctuations are small. Explicitly, we assume the primary
scaling when

$\displaystyle\frac{\rho_{o}}{\ell_{\rm c}}\sim
1\,,\quad\quad\frac{r_{o}}{\ell_{\rm c}}\sim 1\,,\quad\quad\frac{L}{\ell_{\rm
c}}\sim\alpha^{-1}\,,\quad\quad\frac{\lambda_{o}}{\ell_{\rm
c}}\sim\alpha\,,\quad\quad\sigma^{2}\sim\alpha^{3}\,,$

where $\alpha$ is a small dimensionless parameter. We introduce dimensionless
coordinates by:

$\displaystyle{{\itbf x}}={\ell_{\rm c}}{{{\itbf x}}^{\prime}},\quad\quad
z=Lz^{\prime},\quad\quad k_{o}=\frac{k_{o}^{\prime}}{\ell_{\rm c}\alpha},$
$\displaystyle m_{\varepsilon}(Lz^{\prime},\ell_{\rm c}{{\itbf
x}}^{\prime})=\alpha^{3/2}m_{\varepsilon}^{\prime}\left(z^{\prime},{{\itbf
x}}^{\prime}\right),\quad\quad m_{\mu}(Lz^{\prime},\ell_{\rm c}{{\itbf
x}}^{\prime})=\alpha^{3/2}m_{\mu}^{\prime}\left(z^{\prime},{{\itbf
x}}^{\prime}\right).$

We look for the behavior of the coefficient $u_{j}(z^{\prime},{{\itbf
x}}^{\prime})=a_{j}\big{(}z^{\prime}L,{{\itbf x}}^{\prime}\ell_{\rm c}\big{)}$
for long propagation distances of the order of $\alpha^{-1}$. We obtain a
Schrödinger-type equation in which the potential fluctuates in $z^{\prime}$ on
the scale $\alpha$ and is of amplitude $\alpha^{-1/2}$. This diffusion
approximation scaling gives the Brownian field and the model (16). As follows
from our analysis in [18], the backward propagating wave components $b_{j}$,
$j=1,2$ are small compared to the forward propagating wave components $a_{j}$,
$j=1,2$ in this forward beam propagation regime, We remark also that the local
propagation speed is

$c=\frac{1}{\sqrt{\mu^{\prime}{\epsilon}^{\prime}}}=c_{o}\big{[}1-\alpha^{3/2}\frac{m^{\prime}_{\mu}+m^{\prime}_{\epsilon}}{2}+O(\alpha^{3})\big{]},$

and the local impedance

$\zeta=\sqrt{\mu^{\prime}/{\epsilon}^{\prime}}=\zeta_{o}\big{[}1+\alpha^{3/2}\frac{m^{\prime}_{\mu}-m^{\prime}_{\epsilon}}{2}+O(\alpha^{3})\big{]}.$

In view of (18) it then follows that the effective Brownian field is
determined by the fluctuations of the local propagation speed, but not by the
fluctuations of the local impedance.

In Appendix D we address the subsequent scaling regime in which the
correlation length of the medium $\ell_{\rm c}$ and the correlation radius
$\rho_{o}$ of the source are much smaller than the beam radius $r_{o}$ of the
source, moreover, the medium fluctuations are weak and the beam propagates
deep into the medium. We then get the modified scaling picture

$\frac{\rho_{o}}{\ell_{\rm c}}\sim 1\,,\quad\quad\frac{r_{o}}{\ell_{\rm
c}}\sim{\varepsilon}^{-1}\,,\quad\quad\frac{L}{\ell_{\rm
c}}\sim\alpha^{-1}{\varepsilon}^{-1}\,,\quad\quad\frac{\lambda_{o}}{\ell_{\rm
c}}\sim\alpha\,,\quad\quad\sigma^{2}\sim\alpha^{3}{\varepsilon}\,,$ (31)

and we assume $\alpha\ll{\varepsilon}\ll 1$. This means that the paraxial
white-noise limit $\alpha\to 0$ is taken first, and we find

$2ik_{o}d{u}_{j}^{\varepsilon}+\Delta_{{{\itbf
x}}}{u}_{j}^{\varepsilon}\,dz+k_{o}^{2}{u}_{j}^{\varepsilon}\circ
d{B}^{\varepsilon}(z,{{\itbf x}})=0,$

where the Brownian field ${B}^{\varepsilon}$ has covariance $C^{\varepsilon}$.
Then the limit ${\varepsilon}\to 0$ is applied, corresponding to the
scintillation regime. In the scintillation regime (31) the effective strength
$k_{o}^{2}C^{\varepsilon}({\bf 0})L$ of the Brownian field is of order one
since $\sigma^{2}\ell_{\rm c}L/\lambda_{o}^{2}\sim 1$. We also have that
$L\lambda_{o}/r_{o}^{2}$ is of order ${\varepsilon}$. That is, the typical
propagation distance is smaller than the Rayleigh length associated to a
coherent beam with radius $r_{o}$. Here the Rayleigh length corresponds to the
distance when the transverse radius of a coherent beam with radius $r_{o}$ has
roughly doubled by diffraction in the homogeneous medium case and it is given
by $r_{o}^{2}/\lambda_{o}$. The typical propagation distance is, however, of
the same order as the Rayleigh length associated to a partially coherent beam
with beam radius $r_{o}$ and correlation radius $\rho_{o}$, which is given by
$r_{o}\rho_{o}/\lambda_{o}$ [14]. The scintillation regime is, therefore, a
regime where diffractive and random effects are both effective and their
combination results in non-trivial effects. In this regime we are also able to
derive explicit expressions for the fourth moment, see Appendix D.

##  Appendix B The Mean Polarized-Wigner Transform for a Partially Coherent
Beam

We discuss here the Wigner transform which is convenient in order to describe
second order field moments. Let $j,l\in\\{1,2\\}$. The mean Wigner transform
is defined by

${W}_{\rm m}(z,{{\itbf
r}},{\boldsymbol{\xi}}):=\int_{\mathbb{R}^{2}}\exp\big{(}-i{\boldsymbol{\xi}}\cdot{{\itbf
q}}\big{)}\mathbb{E}\left[{u_{j}}\big{(}z,{{\itbf r}}+\frac{{{\itbf
q}}}{2}\big{)}\overline{u}_{l}\big{(}z,{{\itbf r}}-\frac{{{\itbf
q}}}{2}\big{)}\right]d{{\itbf q}}.$ (32)

In view of Eq. (20) it satisfies the closed system

$\frac{\partial{W}_{\rm m}}{\partial
z}+\frac{1}{k_{o}}{{\boldsymbol{\xi}}}\cdot\nabla_{{{\itbf r}}}W_{\rm
m}=\frac{k_{o}^{2}}{4(2\pi)^{2}}\int_{\mathbb{R}^{2}}\hat{C}({{\itbf
k}})\Big{[}W_{\rm m}({\boldsymbol{\xi}}-{{\itbf k}})-W_{\rm
m}({\boldsymbol{\xi}})\Big{]}d{{\itbf k}},$ (33)

starting from ${W}_{\rm m}(z=0,{{\itbf r}},{\boldsymbol{\xi}})=W_{{\rm
m}0}({{\itbf r}},{\boldsymbol{\xi}})$, which is the mean Wigner transform of
the source $(f_{j},f_{l})$:

$W_{{\rm m}0}({{\itbf
r}},{\boldsymbol{\xi}}):=\int_{\mathbb{R}^{2}}\exp\big{(}-i{\boldsymbol{\xi}}\cdot{{\itbf
q}}\big{)}\mathbb{E}\Big{[}f_{j}\big{(}{{\itbf r}}+\frac{{{\itbf
q}}}{2}\big{)}\overline{f}_{l}\big{(}{{\itbf r}}-\frac{{{\itbf
q}}}{2}\big{)}\Big{]}d{{\itbf q}}.$

The transport equation (33) can be solved and we find

$\displaystyle W_{\rm m}(z,{{\itbf r}},{\boldsymbol{\xi}})=$
$\displaystyle\frac{1}{(2\pi)^{2}}\iint_{\mathbb{R}^{2}\times\mathbb{R}^{2}}\exp\Big{(}i{\boldsymbol{\zeta}}\cdot\big{(}{{\itbf
r}}-{\boldsymbol{\xi}}\frac{z}{k_{o}}\big{)}-i{\boldsymbol{\xi}}\cdot{{\itbf
q}}\Big{)}\hat{W}_{{\rm m}0}\big{(}{\boldsymbol{\zeta}},{{\itbf q}}\big{)}$
$\displaystyle\times\exp\Big{(}\frac{k_{o}^{2}}{4}\int_{0}^{z}{C}\big{(}{{\itbf
q}}+{\boldsymbol{\zeta}}\frac{z^{\prime}}{k_{o}}\big{)}-{C}({\bf
0})dz^{\prime}\Big{)}d{\boldsymbol{\zeta}}d{{\itbf q}},$ (34)

where $\hat{W}_{{\rm m}0}$ is defined in terms of the source $(f_{j},f_{l})$
as:

$\hat{W}_{{\rm m}0}({\boldsymbol{\zeta}},{{\itbf
q}})=\int_{\mathbb{R}^{2}}\exp\big{(}-i{\boldsymbol{\zeta}}\cdot{{\itbf
r}}\big{)}\mathbb{E}\Big{[}f_{j}\big{(}{{\itbf r}}+\frac{{{\itbf
q}}}{2}\big{)}\overline{f}_{l}\big{(}{{\itbf r}}-\frac{{{\itbf
q}}}{2}\big{)}\Big{]}d{{\itbf r}}.$ (35)

The mean Wigner transform gives an equation for the mutual coherence function
and we next discuss this in a situation with strong scattering.

##  Appendix C The Mutual Coherence Function in the Strongly Scattering Regime

We consider a Gauss-Schell model for the source, which is a field with
Gaussian statistics, mean zero, and covariance function (14). By taking the
inverse Fourier transform of the mean Wigner transform, we find that the
covariance function of the transmitted field has the form (3). We discuss next
this mutual coherence function in more detail and to find explicit expressions
we assume in the rest of this section that scattering is strong and smooth, in
the sense that

$\displaystyle k_{o}^{2}C({\bf 0})z\gg 1,$ (36) $\displaystyle C({{\itbf
x}})=C({\bf 0})-\bar{\gamma}_{2}|{{\itbf x}}|^{2}+o(|{{\itbf x}}|^{2}).$ (37)

From (16) written in Itô form [34] it follows that the scattering mean free
path $\ell_{\rm mfp}$ (that is the typical propagation distance over which a
coherent wave becomes incoherent) is

$\ell_{\rm mfp}^{-1}=\frac{k_{o}^{2}C({\bf
0})}{8}=\frac{k_{o}^{2}\sigma^{2}\ell_{z}}{8}.$ (38)

Thus, in the regime (36) the propagation distance is large compared to the
scattering mean free path. Note moreover that $\bar{\gamma}_{2}$ can be
interpreted as

$\bar{\gamma}_{2}=\frac{\sigma^{2}\ell_{z}}{\ell_{\perp}^{2}}$ (39)

where $\sigma^{2}$ is the variance of medium fluctuations as above, $\ell_{z}$
is the longitudinal correlation length of the medium (such that $C({\bf
0})=\sigma^{2}\ell_{z}$), and $\ell_{\perp}$ is its transverse correlation
radius of the medium defined by:

$\ell_{\perp}^{-2}=-\frac{\Delta C({\bf 0})}{4C({\bf 0})}.$ (40)

If the medium is isotropic (as assumed in Appendix A), for instance such that
$\mathbb{E}[(m_{\epsilon}+m_{\mu})({{\itbf x}}^{\prime}+z,{{\itbf
x}}^{\prime}+z)(m_{\epsilon}+m_{\mu})({{\itbf
x}}^{\prime},z^{\prime})]=\sigma^{2}\exp(-|{{\itbf x}}|^{2}/\ell_{\rm
c}^{2}-z^{2}/\ell_{\rm c}^{2})$, then we have $\ell_{z}=\sqrt{\pi}\ell_{\rm
c}$ and $\ell_{\perp}=\ell_{\rm c}$.

If (36-37) hold, then the mean intensity is given by (26). This is found via a
Gaussian calculation after inserting (36) in (4). In this expression we can
identify the original beam radius $r_{o}^{2}$, the spreading due to
diffraction $\frac{z^{2}}{k_{o}^{2}\rho_{o}^{2}}$, and the spreading due to
random scattering $\frac{\bar{\gamma}_{2}z^{3}}{3}$. The covariance function
of the field (or mutual coherence function) is given by

$\displaystyle\mathbb{E}\left[{u}_{j}\big{(}z,{{\itbf r}}+\frac{{{\itbf
q}}}{2}\big{)}\overline{u_{j}}\big{(}z,{{\itbf r}}-\frac{{{\itbf
q}}}{2}\big{)}\right]=$ $\displaystyle
A_{j}^{2}\frac{r_{o}^{2}}{R(z;\rho_{o})^{2}}\exp\Big{(}-\frac{|{{\itbf
r}}|^{2}}{R^{2}(z;\rho_{o})}-\frac{|{{\itbf
q}}|^{2}}{4\rho^{2}(z;\rho_{o})}+i\theta(z;\rho_{o}){{\itbf r}}\cdot{{\itbf
q}}\Big{)},$

where the correlation radius of the beam, $\rho(z;\rho_{o})$, is given by
(29), the beam radius, $R(z;\rho_{o})$, by (27) and

$\theta(z;\rho_{o})=\frac{\frac{z}{k_{o}\rho_{o}^{2}}+\frac{k_{o}\bar{\gamma}_{2}z^{2}}{2}}{r_{o}^{2}+\frac{\bar{\gamma}_{2}z^{3}}{3}+\frac{z^{2}}{k_{o}^{2}\rho_{o}^{2}}}.$
(41)

If the source is coherent $\rho_{o}=r_{o}$, then we recover the classical
result obtained in [17], while if the medium is homogeneous and the source is
partially coherent $\rho_{o}<r_{o}$, then we recover the result obtained in
[14]. Note that if the medium is homogeneous then
$1/\sqrt{\theta(z;\rho_{o})}\ll\rho(z;\rho_{o})\leq
R(z;\rho_{o}),~{}\hbox{as}~{}z\to\infty$, while in the random case
$\rho(z;\rho_{o})\ll 1/\sqrt{\theta(z;\rho_{o})}\ll
R(z;\rho_{o}),~{}\hbox{as}~{}z\to\infty$. In fact in both cases we have
$1/\sqrt{\theta}\sim\sqrt{\lambda_{o}z}$, up to a constant. Thus, the coherent
phase modulation is slow relative to the field decorrelation scale for deep
probing in the random medium.

For large propagation distance so that the spreading due to the random medium
dominates,
$z\gg\max\big{(}(r_{o}^{2}/\bar{\gamma}_{2})^{1/3},1/(\bar{\gamma}_{2}k_{o}^{2}\rho_{o}^{2})\big{)}$,
we have:

$\displaystyle R(z;\rho_{o})$ $\displaystyle\sim\mathcal{A}_{\rm
tr}(z):=\sqrt{\frac{\bar{\gamma}_{2}z^{3}}{3}},$ (42)
$\displaystyle\rho(z;\rho_{o})$ $\displaystyle\sim\rho_{\rm
tr}(z):=\frac{1}{k_{o}\sqrt{\bar{\gamma}_{2}z}},$ (43)
$\displaystyle\theta(z;\rho_{o})$ $\displaystyle\sim\frac{3k_{o}}{2z}.$ (44)

Note that these parameters are independent of the parameters $r_{o}$ and
$\rho_{o}$ of the partially coherent source, so that information about the
source is “forgotten” in the case of deep probing. We refer to the parameters
$\mathcal{A}_{\rm tr},\rho_{\rm tr}$ as the time reversal aperture and
resolution respectively. We note that we have the Rayleigh resolution relation

$\displaystyle\rho_{\rm tr}(L)=\frac{\lambda_{o}L}{\mathcal{A}_{\rm tr}(L)}.$
(45)

where $\rho_{\rm tr}$ corresponds to the refocusing resolution one obtains
when a point source emits a wave which is captured on a time reversal mirror
at distance $L$ and reemitted (after time reversal) toward the source. Then it
will refocus at the original source location with a resolution of the order of
the lateral correlation range $\rho_{\rm tr}(L)$ essentially independently of
the actual physical aperture [22]. This can be understood in that the
propagator of the transmitted wave decorrelates (laterally) on this scale and
the refocused wave is essentially the convolution of the propagator with
itself.

We remark that the field correlation radius $\rho(z;\rho_{o})$ is important in
determining statistical stability. If we average a field quantity over an
aperture then the signal-to-noise ratio will in general depend on the ratio of
the aperture to the field correlation radius. We remark finally that we have
$\rho_{\rm tr}\approx\lambda_{o}$ when $z\approx 1/{\bar{\gamma}_{2}}$, which
means that the paraxial approximation is not valid beyond this propagation
distance.

##  Appendix D The Intensity Covariance Function

In this appendix we derive the expression for the intensity covariance
function in the scintillation regime. We start by introducing

$\mu_{4,j_{1}l_{1}j_{2}l_{2}}(z,{{\itbf x}}_{1},{{\itbf x}}_{2},{{\itbf
y}}_{1},{{\itbf y}}_{2})=\mathbb{E}\big{[}u_{j_{1}}(z,{{\itbf
x}}_{1})\overline{u_{l_{1}}(z,{{\itbf y}}_{1})}u_{j_{2}}(z,{{\itbf
x}}_{2})\overline{u_{l_{2}}(z,{{\itbf y}}_{2})}\big{]}.$ (46)

We are interested in the second-order moment of the intensity:

$\mathbb{E}\big{[}I(z,{{\itbf x}}_{1})I(z,{{\itbf
x}}_{2})\big{]}=\sum_{j_{1},j_{2}=1}^{2}\mu_{4,j_{1}j_{1}j_{2}j_{2}}(z,{{\itbf
x}}_{1},{{\itbf x}}_{2},{{\itbf x}}_{1},{{\itbf x}}_{2}).$ (47)

We find using Eq. (16) that the general fourth-order moment
$\mu_{4,j_{1}l_{1}j_{2}l_{2}}$ satisfies the equation

$\displaystyle\frac{\partial\mu_{4,j_{1}l_{1}j_{2}l_{2}}}{\partial
z}=\frac{i}{2k_{o}}\Big{(}\Delta_{{{\itbf x}}_{1}}+\Delta_{{{\itbf
x}}_{2}}-\Delta_{{{\itbf y}}_{1}}-\Delta_{{{\itbf
y}}_{2}}\Big{)}\mu_{4,j_{1}l_{1}j_{2}l_{2}}$
$\displaystyle\hbox{}+\frac{k_{o}^{2}}{4}{U}_{4}\big{(}{{\itbf x}}_{1},{{\itbf
x}}_{2},{{\itbf y}}_{1},{{\itbf y}}_{2}\big{)}\mu_{4,j_{1}l_{1}j_{2}l_{2}},$
(48)

with the generalized potential

$\displaystyle{U}_{4}\big{(}{{\itbf x}}_{1},{{\itbf x}}_{2},{{\itbf
y}}_{1},{{\itbf y}}_{2}\big{)}$ $\displaystyle=$ $\displaystyle{C}({{\itbf
x}}_{1}-{{\itbf y}}_{1})+{C}({{\itbf x}}_{2}-{{\itbf y}}_{2})+{C}({{\itbf
x}}_{1}-{{\itbf y}}_{2})+{C}({{\itbf x}}_{2}-{{\itbf y}}_{1})$ (49)
$\displaystyle-{C}({{\itbf x}}_{1}-{{\itbf x}}_{2})-{C}({{\itbf
y}}_{1}-{{\itbf y}}_{2})-2{C}({\bf 0}),$

and the initial condition:

$\mu_{4,j_{1}l_{1}j_{2}l_{2}}(z=0,{{\itbf x}}_{1},{{\itbf x}}_{2},{{\itbf
y}}_{1},{{\itbf y}}_{2})=\mathbb{E}\big{[}f_{j_{1}}({{\itbf
x}}_{1})\overline{f_{l_{1}}({{\itbf y}}_{1})}f_{j_{2}}({{\itbf
x}}_{2})\overline{f_{l_{2}}({{\itbf y}}_{2})}\big{]}.$

Using the Gaussian property of the source, the initial condition for the
fourth-order moment is:

$\displaystyle\mu_{4,j_{1}l_{1}j_{2}l_{2}}(z=0,{{\itbf x}}_{1},{{\itbf
x}}_{2},{{\itbf y}}_{1},{{\itbf y}}_{2})$
$\displaystyle=\mathbb{E}\big{[}f_{j_{1}}({{\itbf
x}}_{1})\overline{f_{l_{1}}}({{\itbf
y}}_{1})\big{]}\mathbb{E}\big{[}f_{j_{2}}({{\itbf
x}}_{2})\overline{f_{l_{2}}}({{\itbf
y}}_{2})\big{]}+\mathbb{E}\big{[}f_{j_{1}}({{\itbf
x}}_{1})\overline{f_{l_{2}}}({{\itbf
y}}_{2})\big{]}\mathbb{E}\big{[}f_{j_{2}}({{\itbf
x}}_{2})\overline{f_{l_{1}}}({{\itbf y}}_{1})\big{]},$

where the covariance function of the source is given by (2).

We parameterize the four points ${{\itbf x}}_{1},{{\itbf x}}_{2},{{\itbf
y}}_{1},{{\itbf y}}_{2}$ in (46) in the special way:

$\displaystyle{{\itbf x}}_{1}=\frac{{{\itbf r}}_{1}+{{\itbf r}}_{2}+{{\itbf
q}}_{1}+{{\itbf q}}_{2}}{2},\quad\quad{{\itbf y}}_{1}=\frac{{{\itbf
r}}_{1}+{{\itbf r}}_{2}-{{\itbf q}}_{1}-{{\itbf q}}_{2}}{2},$ (50)
$\displaystyle{{\itbf x}}_{2}=\frac{{{\itbf r}}_{1}-{{\itbf r}}_{2}+{{\itbf
q}}_{1}-{{\itbf q}}_{2}}{2},\quad\quad{{\itbf y}}_{2}=\frac{{{\itbf
r}}_{1}-{{\itbf r}}_{2}-{{\itbf q}}_{1}+{{\itbf q}}_{2}}{2}.$ (51)

In particular ${{\itbf r}}_{1}/2$ is the barycenter of the four points
${{\itbf x}}_{1},{{\itbf x}}_{2},{{\itbf y}}_{1},{{\itbf y}}_{2}$:

$\displaystyle{{\itbf r}}_{1}=\frac{{{\itbf x}}_{1}+{{\itbf x}}_{2}+{{\itbf
y}}_{1}+{{\itbf y}}_{2}}{2},\quad\quad{{\itbf q}}_{1}=\frac{{{\itbf
x}}_{1}+{{\itbf x}}_{2}-{{\itbf y}}_{1}-{{\itbf y}}_{2}}{2},$
$\displaystyle{{\itbf r}}_{2}=\frac{{{\itbf x}}_{1}-{{\itbf x}}_{2}+{{\itbf
y}}_{1}-{{\itbf y}}_{2}}{2},\quad\quad{{\itbf q}}_{2}=\frac{{{\itbf
x}}_{1}-{{\itbf x}}_{2}-{{\itbf y}}_{1}+{{\itbf y}}_{2}}{2}.$

We denote by $\mu$ the fourth-order moment in these new variables (without
writing the dependence on $j_{1},j_{2},l_{1},l_{2}$):

$\mu(z,{{\itbf q}}_{1},{{\itbf q}}_{2},{{\itbf r}}_{1},{{\itbf
r}}_{2})=\mu_{4,j_{1}l_{1}j_{2}l_{2}}(z,{{\itbf x}}_{1},{{\itbf
x}}_{2},{{\itbf y}}_{1},{{\itbf y}}_{2})$ (52)

with ${{\itbf x}}_{1},{{\itbf x}}_{2},{{\itbf y}}_{1},{{\itbf y}}_{2}$ given
by (50-51) in terms of ${{\itbf q}}_{1},{{\itbf q}}_{2},{{\itbf
r}}_{1},{{\itbf r}}_{2}$.

In the variables $({{\itbf q}}_{1},{{\itbf q}}_{2},{{\itbf r}}_{1},{{\itbf
r}}_{2})$ the function ${\mu}$ satisfies the system:

$\frac{\partial{\mu}}{\partial z}=\frac{i}{k_{o}}\big{(}\nabla_{{{\itbf
r}}_{1}}\cdot\nabla_{{{\itbf q}}_{1}}+\nabla_{{{\itbf
r}}_{2}}\cdot\nabla_{{{\itbf
q}}_{2}}\big{)}{\mu}+\frac{k_{o}^{2}}{4}{U}({{\itbf q}}_{1},{{\itbf
q}}_{2},{{\itbf r}}_{1},{{\itbf r}}_{2}){\mu},$ (53)

with the generalized potential

$\displaystyle{U}({{\itbf q}}_{1},{{\itbf q}}_{2},{{\itbf r}}_{1},{{\itbf
r}}_{2})$ $\displaystyle=$ $\displaystyle{C}({{\itbf q}}_{2}+{{\itbf
q}}_{1})+{C}({{\itbf q}}_{2}-{{\itbf q}}_{1})+{C}({{\itbf r}}_{2}+{{\itbf
q}}_{1})+{C}({{\itbf r}}_{2}-{{\itbf q}}_{1})$ (54) $\displaystyle-{C}({{\itbf
q}}_{2}+{{\itbf r}}_{2})-{C}({{\itbf q}}_{2}-{{\itbf r}}_{2})-2{C}({\bf 0}).$

Note in particular that the generalized potential does not depend on the
barycenter ${{\itbf r}}_{1}$, and this comes from the fact that the medium is
statistically homogeneous. The Fourier transform (in ${{\itbf q}}_{1}$,
${{\itbf q}}_{2}$, ${{\itbf r}}_{1}$, and ${{\itbf r}}_{2}$) of the fourth-
order moment is defined by:

$\displaystyle\hat{\mu}(z,{\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2})$
$\displaystyle=$
$\displaystyle\iint_{\mathbb{R}^{2}\times\mathbb{R}^{2}}{\mu}(z,{{\itbf
q}}_{1},{{\itbf q}}_{2},{{\itbf r}}_{1},{{\itbf r}}_{2})$ (55)
$\displaystyle\times\exp\big{(}-i{{\itbf
q}}_{1}\cdot{\boldsymbol{\xi}}_{1}-i{{\itbf
q}}_{2}\cdot{\boldsymbol{\xi}}_{2}-i{{\itbf
r}}_{1}\cdot{\boldsymbol{\zeta}}_{1}-i{{\itbf
r}}_{2}\cdot{\boldsymbol{\zeta}}_{2}\big{)}d{{\itbf q}}_{1}d{{\itbf
q}}_{2}d{{\itbf r}}_{1}d{{\itbf r}}_{2}.\hskip 21.68121pt$

It satisfies

$\displaystyle\frac{\partial\hat{\mu}}{\partial
z}+\frac{i}{k_{o}}\big{(}{\boldsymbol{\xi}}_{1}\cdot{\boldsymbol{\zeta}}_{1}+{\boldsymbol{\xi}}_{2}\cdot{\boldsymbol{\zeta}}_{2}\big{)}\hat{\mu}=\frac{k_{o}^{2}}{4(2\pi)^{2}}\int_{\mathbb{R}^{2}}\hat{C}({{\itbf
k}})\bigg{[}\hat{\mu}({\boldsymbol{\xi}}_{1}-{{\itbf
k}},{\boldsymbol{\xi}}_{2}-{{\itbf
k}},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2})$
$\displaystyle\quad+\hat{\mu}({\boldsymbol{\xi}}_{1}+{{\itbf
k}},{\boldsymbol{\xi}}_{2}-{{\itbf
k}},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2})-2\hat{\mu}({\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2})$
$\displaystyle\quad+\hat{\mu}({\boldsymbol{\xi}}_{1}+{{\itbf
k}},{\boldsymbol{\xi}}_{2},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2}-{{\itbf
k}})+\hat{\mu}({\boldsymbol{\xi}}_{1}-{{\itbf
k}},{\boldsymbol{\xi}}_{2},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2}-{{\itbf
k}})$
$\displaystyle\quad-\hat{\mu}({\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2}-{{\itbf
k}},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2}-{{\itbf
k}})-\hat{\mu}({\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2}+{{\itbf
k}},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2}-{{\itbf
k}})\bigg{]}d{{\itbf k}}.$ (56)

If $j_{1}=j_{2}=l_{2}=l_{2}\equiv j\in\\{1,2\\}$, then the initial condition
is

$\displaystyle\hat{\mu}(z=0,{\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2})$
$\displaystyle=$
$\displaystyle(2\pi)^{4}A_{j}^{4}\phi_{\rho_{o}}^{1}({\boldsymbol{\xi}}_{1})\phi_{\rho_{o}}^{1}({\boldsymbol{\xi}}_{2})\phi_{r_{o}}^{1}({\boldsymbol{\zeta}}_{1})\phi_{r_{o}}^{1}({\boldsymbol{\zeta}}_{2})$
(57)
$\displaystyle+(2\pi)^{4}A_{j}^{4}\phi_{\rho_{o}}^{1}({\boldsymbol{\xi}}_{1})\phi_{\rho_{o}}^{1}({\boldsymbol{\zeta}}_{2})\phi_{r_{o}}^{1}({\boldsymbol{\zeta}}_{1})\phi_{r_{o}}^{1}({\boldsymbol{\xi}}_{2}),$

with

$\phi^{1}_{\rho}({\boldsymbol{\xi}})=\frac{\rho^{2}}{2\pi}\exp\Big{(}-\frac{\rho^{2}|{\boldsymbol{\xi}}|^{2}}{2}\Big{)}.$
(58)

Similar Gaussian expressions hold for the initial condition in the other cases
for $(j_{1},j_{2},l_{1},l_{2})$, we only address the case
$j_{1}=j_{2}=l_{2}=l_{2}\equiv j$ in the following.

We cannot solve the problem for the fourth moment $\mu$ explicitly and
consider next a secondary scaling limit where we can identify an explicit
solution. We consider the scintillation regime, discussed in more detail in
Appendix A, where the correlation radius of the source is of the same order as
the correlation radius of the medium, but the beam radius of the source is
much larger:

$\rho_{o}\to{\rho_{o}},\quad C({{\itbf x}})\to{\varepsilon}C({{\itbf
x}}),\quad r_{o}\to\frac{r_{o}}{{\varepsilon}},\quad
z\to\frac{z}{{\varepsilon}}.$ (59)

We introduce the rescaled function

$\tilde{\mu}^{\varepsilon}(z,{\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2})=\mu\Big{(}\frac{z}{{\varepsilon}},{\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2}\Big{)}\exp\Big{(}i\frac{z}{k_{o}{\varepsilon}}({\boldsymbol{\xi}}_{1}\cdot{\boldsymbol{\zeta}}_{1}+{\boldsymbol{\xi}}_{2}\cdot{\boldsymbol{\zeta}}_{2})\Big{)}.$
(60)

Then the limit ${\varepsilon}\to 0$ is applied, corresponding to the
scintillation regime.

In the scintillation regime the rescaled function $\tilde{\mu}^{\varepsilon}$
satisfies the equation with fast phases

$\displaystyle\frac{\partial\tilde{\mu}^{\varepsilon}}{\partial z}$
$\displaystyle=$
$\displaystyle\frac{k_{o}^{2}}{4(2\pi)^{2}}\int_{\mathbb{R}^{2}}\hat{C}({{\itbf
k}})\bigg{[}-2\tilde{\mu}^{\varepsilon}({\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2})$
(61) $\displaystyle+\tilde{\mu}^{\varepsilon}({\boldsymbol{\xi}}_{1}-{{\itbf
k}},{\boldsymbol{\xi}}_{2}-{{\itbf
k}},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2})e^{i\frac{c_{o}z}{{\varepsilon}\omega_{o}}{{\itbf
k}}\cdot({\boldsymbol{\zeta}}_{2}+{\boldsymbol{\zeta}}_{1})}$
$\displaystyle+\tilde{\mu}^{\varepsilon}({\boldsymbol{\xi}}_{1}-{{\itbf
k}},{\boldsymbol{\xi}}_{2},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2}-{{\itbf
k}})e^{i\frac{c_{o}z}{{\varepsilon}\omega_{o}}{{\itbf
k}}\cdot({\boldsymbol{\xi}}_{2}+{\boldsymbol{\zeta}}_{1})}$
$\displaystyle+\tilde{\mu}^{\varepsilon}({\boldsymbol{\xi}}_{1}+{{\itbf
k}},{\boldsymbol{\xi}}_{2}-{{\itbf
k}},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2})e^{i\frac{c_{o}z}{{\varepsilon}\omega_{o}}{{\itbf
k}}\cdot({\boldsymbol{\zeta}}_{2}-{\boldsymbol{\zeta}}_{1})}$
$\displaystyle+\tilde{\mu}^{\varepsilon}({\boldsymbol{\xi}}_{1}+{{\itbf
k}},{\boldsymbol{\xi}}_{2},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2}-{{\itbf
k}})e^{i\frac{c_{o}z}{{\varepsilon}\omega_{o}}{{\itbf
k}}\cdot({\boldsymbol{\xi}}_{2}-{\boldsymbol{\zeta}}_{1})}$
$\displaystyle-\tilde{\mu}^{\varepsilon}({\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2}-{{\itbf
k}},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2}-{{\itbf
k}})e^{i\frac{c_{o}z}{{\varepsilon}\omega_{o}}({{\itbf
k}}\cdot({\boldsymbol{\zeta}}_{2}+{\boldsymbol{\xi}}_{2})-|{{\itbf k}}|^{2})}$
$\displaystyle-\tilde{\mu}^{\varepsilon}({\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2}-{{\itbf
k}},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2}+{{\itbf
k}})e^{i\frac{c_{o}z}{{\varepsilon}\omega_{o}}({{\itbf
k}}\cdot({\boldsymbol{\zeta}}_{2}-{\boldsymbol{\xi}}_{2})+|{{\itbf
k}}|^{2})}\bigg{]}d{{\itbf k}},$

starting from

$\displaystyle\tilde{\mu}^{\varepsilon}(z=0,{\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2})$
$\displaystyle=$
$\displaystyle(2\pi)^{8}A_{j}^{4}\phi^{1}_{\rho_{o}}({\boldsymbol{\xi}}_{1})\phi^{1}_{\rho_{o}}({\boldsymbol{\xi}}_{2})\phi^{\varepsilon}_{r_{o}}({\boldsymbol{\zeta}}_{1})\phi^{\varepsilon}_{r_{o}}({\boldsymbol{\zeta}}_{2})$
(62)
$\displaystyle+(2\pi)^{8}A_{j}^{4}\phi^{1}_{\rho_{o}}({\boldsymbol{\xi}}_{1})\phi^{1}_{\rho_{o}}({\boldsymbol{\zeta}}_{2})\phi^{\varepsilon}_{r_{o}}({\boldsymbol{\zeta}}_{1})\phi^{\varepsilon}_{r_{o}}({\boldsymbol{\xi}}_{2}),$

where $\phi^{\varepsilon}_{\rho}$ is defined by:

$\phi^{\varepsilon}_{\rho}({\boldsymbol{\xi}})=\frac{\rho^{2}}{2\pi{\varepsilon}^{2}}\exp\Big{(}-\frac{\rho^{2}}{2{\varepsilon}^{2}}|{\boldsymbol{\xi}}|^{2}\Big{)}.$
(63)

The following result shows that $\tilde{\mu}^{\varepsilon}$ exhibits a multi-
scale behavior as ${\varepsilon}\to 0$, with some components evolving at the
scale ${\varepsilon}$ and some components evolving on the order one scale
[20].

###### Proposition D.1

The function
$\tilde{\mu}^{\varepsilon}(z,{\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2})$
can be expanded as

$\displaystyle\tilde{\mu}^{\varepsilon}\big{(}z,{\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2}\big{)}$
$\displaystyle=$
$\displaystyle(2\pi)^{8}A_{j}^{4}\phi^{\varepsilon}_{r_{o}}({\boldsymbol{\zeta}}_{1})\phi^{\varepsilon}_{r_{o}}({\boldsymbol{\zeta}}_{2})B\big{(}z,{\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2},\frac{{\boldsymbol{\zeta}}_{1}}{{\varepsilon}},\frac{{\boldsymbol{\zeta}}_{2}}{{\varepsilon}}\big{)}$
(64)
$\displaystyle+(2\pi)^{8}A_{j}^{4}\phi^{\varepsilon}_{r_{o}}({\boldsymbol{\zeta}}_{1})\phi^{\varepsilon}_{r_{o}}({\boldsymbol{\xi}}_{2})B\big{(}z,{\boldsymbol{\xi}}_{1},{\boldsymbol{\zeta}}_{2},\frac{{\boldsymbol{\zeta}}_{1}}{{\varepsilon}},\frac{{\boldsymbol{\xi}}_{2}}{{\varepsilon}}\big{)}$
$\displaystyle+R^{\varepsilon}(z,{\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2}),$

where

$\displaystyle
B(z,{\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2})=\frac{1}{(2\pi)^{4}}\iint_{\mathbb{R}^{2}\times\mathbb{R}^{2}}d{{\itbf
x}}d{{\itbf y}}\exp\Big{(}-\frac{|{{\itbf x}}|^{2}+|{{\itbf
y}}|^{2}}{2\rho_{o}^{2}}-i{\boldsymbol{\xi}}_{1}\cdot{{\itbf
x}}-i{\boldsymbol{\xi}}_{2}\cdot{{\itbf y}}$
$\displaystyle+\frac{k_{o}^{2}}{4}\int_{0}^{z}C\big{(}{{\itbf x}}+{{\itbf
y}}+\frac{z^{\prime}}{k_{o}}({\boldsymbol{\zeta}}_{1}+{\boldsymbol{\zeta}}_{2})\big{)}+C\big{(}{{\itbf
x}}-{{\itbf
y}}+\frac{z^{\prime}}{k_{o}}({\boldsymbol{\zeta}}_{1}-{\boldsymbol{\zeta}}_{2})\big{)}-2C({\bf
0})dz^{\prime}\Big{)},$ (65)

and the function $R^{\varepsilon}$ goes to $0$ as ${\varepsilon}\to 0$.

As a consequence, the second-order moment of the intensity is

$\displaystyle\mathbb{E}\big{[}\big{|}u\big{(}\frac{z}{{\varepsilon}},\frac{{{\itbf
r}}}{{\varepsilon}}+\frac{{{\itbf
q}}}{2}\big{)}\big{|}^{2}\big{|}u\big{(}\frac{z}{{\varepsilon}},\frac{{{\itbf
r}}}{{\varepsilon}}-\frac{{{\itbf
q}}}{2}\big{)}\big{|}^{2}\big{]}=A_{j}^{4}{\cal I}_{\rho_{o}}(z,2{{\itbf
r}},{\bf 0})+A_{j}^{4}{\cal J}_{\rho_{o}}(z,2{{\itbf r}},{{\itbf q}}),$ (66)

with

$\displaystyle{\cal I}_{\rho_{o}}(z,{{\itbf r}}_{1},{{\itbf r}}_{2})$
$\displaystyle=\iint_{\mathbb{R}^{2}\times\mathbb{R}^{2}}d{\boldsymbol{\xi}}_{1}d{\boldsymbol{\xi}}_{2}d{\boldsymbol{\zeta}}_{1}d{\boldsymbol{\zeta}}_{2}\phi_{r_{o}}^{1}({\boldsymbol{\zeta}}_{1})\phi_{r_{o}}^{1}({\boldsymbol{\zeta}}_{2})B(z,{\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\zeta}}_{2})$
$\displaystyle\quad\times\exp\Big{(}-i\frac{z}{k_{o}}({\boldsymbol{\xi}}_{1}\cdot{\boldsymbol{\zeta}}_{1}+{\boldsymbol{\xi}}_{2}\cdot{\boldsymbol{\zeta}}_{2})+i{\boldsymbol{\zeta}}_{1}\cdot{{\itbf
r}}_{1}+i{\boldsymbol{\zeta}}_{2}\cdot{{\itbf r}}_{2}\Big{)}$
$\displaystyle=\Big{(}\int_{\mathbb{R}^{2}}\phi^{1}_{\frac{r_{o}}{\sqrt{2}}}({\boldsymbol{\zeta}})\exp\Big{(}i{\boldsymbol{\zeta}}\cdot\frac{{{\itbf
r}}_{1}+{{\itbf
r}}_{2}}{2}+\frac{k_{o}^{2}}{4}\int_{0}^{z}C(\frac{{\boldsymbol{\zeta}}z^{\prime}}{k_{o}})-C({\bf
0})dz^{\prime}-\frac{|{\boldsymbol{\zeta}}|^{2}z^{2}}{4k_{o}^{2}\rho_{o}^{2}}\Big{)}d{\boldsymbol{\zeta}}\Big{)}$
$\displaystyle\quad\times\Big{(}\int_{\mathbb{R}^{2}}\phi^{1}_{\frac{r_{o}}{\sqrt{2}}}({\boldsymbol{\zeta}})\exp\Big{(}i{\boldsymbol{\zeta}}\cdot\frac{{{\itbf
r}}_{1}-{{\itbf
r}}_{2}}{2}+\frac{k_{o}^{2}}{4}\int_{0}^{z}C(\frac{{\boldsymbol{\zeta}}z^{\prime}}{k_{o}})-C({\bf
0})dz^{\prime}-\frac{|{\boldsymbol{\zeta}}|^{2}z^{2}}{4k_{o}^{2}\rho_{o}^{2}}\Big{)}d{\boldsymbol{\zeta}}\Big{)}$
$\displaystyle={\cal H}_{\rho_{o}}(z,\frac{{{\itbf r}}_{1}+{{\itbf
r}}_{2}}{2},{\bf 0}){\cal H}_{\rho_{o}}(z,\frac{{{\itbf r}}_{1}-{{\itbf
r}}_{2}}{2},{\bf 0})$ (67)

and

$\displaystyle{\cal J}_{\rho_{o}}(z,{{\itbf r}}_{1},{{\itbf r}}_{2})$
$\displaystyle=\iint_{\mathbb{R}^{2}\times\mathbb{R}^{2}}d{\boldsymbol{\xi}}_{1}d{\boldsymbol{\xi}}_{2}d{\boldsymbol{\zeta}}_{1}d{\boldsymbol{\zeta}}_{2}\phi_{r_{o}}^{1}({\boldsymbol{\zeta}}_{1})\phi_{r_{o}}^{1}({\boldsymbol{\xi}}_{2})B(z,{\boldsymbol{\xi}}_{1},{\boldsymbol{\zeta}}_{2},{\boldsymbol{\zeta}}_{1},{\boldsymbol{\xi}}_{2})$
$\displaystyle\quad\times\exp\Big{(}-i\frac{z}{k_{o}}({\boldsymbol{\xi}}_{1}\cdot{\boldsymbol{\zeta}}_{1}+{\boldsymbol{\xi}}_{2}\cdot{\boldsymbol{\zeta}}_{2})+i{\boldsymbol{\zeta}}_{1}\cdot{{\itbf
r}}_{1}+i{\boldsymbol{\zeta}}_{2}\cdot{{\itbf r}}_{2}\Big{)}$
$\displaystyle=\Big{|}\int_{\mathbb{R}^{2}}\phi^{1}_{\frac{r_{o}}{\sqrt{2}}}({\boldsymbol{\zeta}})\exp\Big{(}i{\boldsymbol{\zeta}}\cdot\frac{{{\itbf
r}}_{1}}{2}+\frac{k_{o}^{2}}{4}\int_{0}^{z}C(\frac{{\boldsymbol{\zeta}}z^{\prime}}{k_{o}}-{{\itbf
r}}_{2})-C({\bf
0})dz^{\prime}-\frac{|\frac{{\boldsymbol{\zeta}}z}{k_{o}}-{{\itbf
r}}_{2}|^{2}}{4\rho_{o}^{2}}\Big{)}d{\boldsymbol{\zeta}}\Big{|}^{2}$
$\displaystyle=\big{|}{\cal H}_{\rho_{o}}(z,\frac{{{\itbf r}}_{1}}{2},{{\itbf
r}}_{2})\big{|}^{2},$ (68)

where ${\cal H}_{\rho_{o}}$ is defined by (4). We finally remark that for far
away points the second-order moment of the intensity is

$\displaystyle\mathbb{E}\left[\big{|}u\big{(}\frac{z}{{\varepsilon}},\frac{{{\itbf
r}}}{{\varepsilon}}+\frac{{{\itbf
q}}}{2{\varepsilon}}\big{)}\big{|}^{2}\big{|}u\big{(}\frac{z}{{\varepsilon}},\frac{{{\itbf
r}}}{{\varepsilon}}-\frac{{{\itbf
q}}}{2{\varepsilon}}\big{)}\big{|}^{2}\right]=\big{(}A_{1}^{2}+A_{2}^{2}\big{)}^{2}{\cal
I}_{\rho_{o}}(z,2{{\itbf r}},{{\itbf q}})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}\hbox{}=\big{(}A_{1}^{2}+A_{2}^{2}\big{)}^{2}{\cal
H}_{\rho_{o}}(z,{{\itbf r}}+\frac{{{\itbf q}}}{2},{\bf 0}){\cal
H}_{\rho_{o}}(z,{{\itbf r}}-\frac{{{\itbf q}}}{2},{\bf 0})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}=\mathbb{E}\left[\big{|}u\big{(}\frac{z}{{\varepsilon}},\frac{{{\itbf
r}}}{{\varepsilon}}+\frac{{{\itbf
q}}}{2{\varepsilon}}\big{)}\big{|}^{2}\right]\mathbb{E}\left[\big{|}u\big{(}\frac{z}{{\varepsilon}},\frac{{{\itbf
r}}}{{\varepsilon}}-\frac{{{\itbf
q}}}{2{\varepsilon}}\big{)}\big{|}^{2}\right],$

so that the intensities then indeed are uncorrelated.

## Acknowledgements

JG was supported by the Agence Nationale pour la Recherche under Grant No.
ANR-19-CE46-0007 (project ICCI), and Air Force Office of Scientific Research
under grant FA9550-18-1-0217.
KS was supported by the Air Force Office of Scientific Research under grant
FA9550-18-1-0217, and the National Science Foundation under grant DMS-2010046.

## References

  * [1] L.C. Andrews and R.L. Philipps, Laser Beam Propagation Through Random Media, SPIE Press, Bellingham, 2005.
  * [2] S. Avramov-Zamurovic, C. Nelson, R. Malek-Madanib, and O. Korotkova, Polarization-induced reduction in scintillation of optical beams propagating in simulated turbulent atmospheric channels, Waves in Random and Complex Media 24:4 452–462 (2014).
  * [3] S. Avramov-Zamurovic, C. Nelson, S. Guth, and O. Korotkova, Flatness parameter influence on scintillation reduction for multi-Gaussian Schell-model beams propagating in turbulent air, Applied Optics 55:13 3442–3446 (2016).
  * [4] G. Baker, Gaussian beam weak scintillation: low-order turbulence effects and applicability of the Rytov method, J. Opt. Soc. Am. A 23:2 395–417 (2006).
  * [5] G. Bal and O. Pinaud, Dynamics of wave scintillation in random media, Comm. Partial Differential Equations 35 1176–1235 (2010).
  * [6] A. Beckus, A. Tamasan, A. Dogariu, A.F. Abouraddy, and G.K. Atia, On the inverse problem of source reconstruction from coherence measurements, J. Opt. Soc. Am. A 35:6 959–968 (2018).
  * [7] L. Borcea, J. Garnier, G. Papanicolaou, and C. Tsogka, Enhanced statistical stability in coherent interferometric imaging, Inverse Problems 27:8 085004 (2011).
  * [8] L. Borcea, G. Papanicolaou and C. Tsogka, Time and direction of arrival detection and filtering for imaging in strongly scattering random media, Waves in Random and Complex Media 27:4 664–689 (2017).
  * [9] M. Born and E. Wolf, Principles of Optics, Cambridge University Press, Cambridge, 1999.
  * [10] M. Charnotskii, Extended Huygens-Fresnel principle and optical waves propagation in turbulence: discussion, J. Opt. Soc. Am. A 32:7 1357–1365 (2015).
  * [11] D. Dawson and G. Papanicolaou, A random wave process, Appl. Math. Optim. 12 97–114 (1984).
  * [12] J.-P. Fouque, J. Garnier, G. Papanicolaou, and K. Sølna, Wave Propagation and Time Reversal in Randomly Layered Media, Springer, New York, 2007.
  * [13] J.-P. Fouque, G. Papanicolaou, and Y. Samuelides, Forward and Markov approximation: the strong-intensity-fluctuations regime revisited, Waves in Random Media 8 303–314 (1998).
  * [14] J. Garnier, C. Gouédard, and L. Videau, Propagation of a partially coherent beam under the interaction of small and large scales, Opt. Commun. 176 281–297 (2000).
  * [15] J. Garnier and G. Papanicolaou, Passive Imaging with Ambient Noise, Cambridge University Press, Cambridge, 2016.
  * [16] J. Garnier and K. Sølna, Coupled paraxial wave equations in random media in the white-noise regime, Ann. Appl. Probab. 19 318–346 (2009).
  * [17] J. Garnier and K. Sølna, Scaling limits for wave pulse transmission and reflection operators, Wave Motion 46 122–143 (2009).
  * [18] J. Garnier and K. Sølna, Paraxial coupling of electromagnetic waves in random media, SIAM J. on Multiscale Modeling and Simulation 7 1928–1955 (2009).
  * [19] J. Garnier and K. Sølna, White-noise paraxial approximation for a general random hyperbolic system, SIAM J. Multiscale Model. Simul. 13:3 1022–1060 (2015).
  * [20] J. Garnier and K. Sølna, Fourth-moment analysis for beam propagation in the white-noise paraxial regime, Arch. Rational Mech. Anal. 220 37–81 (2016).
  * [21] J. Garnier and K. Sølna, Scintillation in the white-noise paraxial regime, Comm. Part. Differ. Equat. 39:4 626–650 (2014).
  * [22] J. Garnier and K. Sølna, Focusing waves through a randomly scattering medium in the white-noise paraxial SIAM Journal of Applied Mathematics 77:2 pp. 500-519 (2017).
  * [23] G. Gbur and E. Wolf, Spreading of partially coherent beams in random media, J. Opt. Soc. Am. A 19:8 1592–1598 (2002).
  * [24] G. Gbur, Partially coherent beam propagation in atmospheric turbulence, J. Opt. Soc. Am. A 31:9 2038–2045 (2014).
  * [25] Y. Gu and T. Komorowski, Gaussian fluctuations from random Schrödinger equation, Comm. Part. Differ. Equat. 46:2 201–232 (2021).
  * [26] A. Ishimaru, Wave Propagation and Scattering in Random Media, IEEE Press, Piscataway, 1997.
  * [27] T. Komorowski and L. Ryzhik, Fluctuations of solutions to Wigner equation with an Ornstein-Uhlenbeck potential, Discrete and Continuous Dynamical Systems-Series B 17:3 871–914 (2012).
  * [28] O. Korotkova, Theoretical Statistical Optics, World Scientific, 2022.
  * [29] O. Korotkova and G. Gbur, Applications of optical coherence theory, Prog. Opt. 65 43–104 (2020).
  * [30] M. Moscoso, A. Novikov, G. Papanicolaou, and C. Tsogka. Synthetic aperture imaging with intensity-only measurements, IEEE Transactions on Computational Imaging 6 87–94 (2020).
  * [31] Y. Liu, Y. Chen, F. Wang, Y. Cai, C. Liang, and O. Korotkova, Robust far-field imaging by spatial coherence engineering, Advances in Opto-Electronics, in press (2022).
  * [32] V.I. Tatarskii, Waves Propagation in a Turbulent Medium, McGraw-Hill, 1961.
  * [33] V.I. Tatarskii, The Effects of Turbulent Atmosphere on Wave Propagation, U.S. Department of Commerce, TT-68-50464, Springfield, 1971.
  * [34] B. Øksendal, Stochastic Differential Equations, An Introduction with Applications, Springer, 2007.
  * [35] D. Voelz and K. Fitzhenry, Pseudo-partially coherent beam for free-space laser communication, Proc. of SPIE 5550 218–224 (2004).
  * [36] F. Wang, X. Liu, and Y. Cai, Propagation of partially coherent beam in turbulent atmosphere: A review, Progress In Electromagnetics Research 150 123–143 (2015).