# A coherent full wave scattering formulation for random layered media

E. Calzetta  Universidad de Buenos Aires, Facultad de Ciencias Exactas y
Naturales, Departamento de Física. Ciudad Universitaria, Int. Güiraldes 2160,
Buenos Aires, CABA, C1428EGA, Argentina  CONICET-Instituto de Física de Buenos
Aires (IFIBA). Buenos Aires, Argentina.  Corresponding author:
<EMAIL_ADDRESS>M. Franco  CONICET-Universidad de Buenos Aires, Instituto
de Astronomía y Física del Espacio (IAFE), Ciudad Universitaria, Av. Cantilo
S/N, Buenos Aires, C1428ZAA, Argentina.

###### Abstract

We present a fully coherent, analytic model of the backscattering intensity in
all HH, HV, VH and VV channels, for the volume scattering of radiation from a
layer of finite thickness, such as a vegetation layer over bare soil. We aim
for a simple, not numerically intensive model which could be used either as
forward model in a Bayesian estimation scheme, or else as a preliminary means
to identify key features of a concrete problem, for its further analysis by
more sophisticated theoretical and numerical approaches.

††journal: josaa

##  1 Introduction

The scattering of electromagnetic waves in random dielectric media is a
subject that continues to be challenging at present. This scenario can be
found both at the microscopic scale where the scattering medium is illuminated
by wavelengths in the visible or infrared spectrum and at regional scales
where the earth’s surface is monitored with orbital or airborne radars.
Whatever the case, the typical lengths at which a medium can be considered
homogeneous or have anisotropies are always defined with respect to the
wavelength of the incident signal. Moreover, to extract useful information
from the data it is necessary to have a theoretical model connecting the
measured quantities to the physical parameters of the target.

Our goal in this paper is to develop such a theoretical model in order to be
able to compute, from the assumed geometric and statistical properties of a
vegetation layer, the intensity of backscattered radiation in all HH, VH, VV
and HV channels.

The two considerations that guide our work are that we want a coherent model,
and a numerically non-intensive model which could be used as an input in a
Bayesian estimation scheme [1, 2, 3]. For this reason we aim for an analytic
model, instead of a model based on intensive numerical simulations. We include
in this latter category approaches where the forward model is used to build a
data cube which is then the input for the estimation scheme. For an extensive
review of both analytical and numerical models see [4, 5, 6, 7].

In the context of remote sensing, where the earth’s surface is monitored at
scales of the order of kilometers, the randomness of the scattering medium can
be described mainly by three factors: the roughness of the surface, the
stratification of the underlying medium, and the vegetation on the surface.
Regarding the first two scattering mechanisms, there is an extensive
literature describing how to calculate the amplitude and phase of the
scattered wave under different types of approaches [5, 8, 9, 10]. The
different types of theories take into account the wavelength of the incident
signal and the parameters indicating the randomness of the different types of
scenarios, the most important being the roughness of the scattering surface or
the inhomogeneities of the material medium. The third scattering mechanism
turns out to be the most complicated to calculate analytically. The vegetation
on the surface is usually modeled as a dielectric medium with random
inhomogeneities or, alternatively as a collection of randomly oriented
cylinders. The propagation and scattering of electromagnetic waves in such
media is a problem in itself. To deal with these issues, two of the most
commonly used approaches used in remote sensing are the Radiative Transfer
Equation (RTE) [11, 12, 4, 10] and the Distorted Born Approximation (DBA) [13,
14, 15, 16]. A well-known result is that for both RTE and DBA cross-
polarization effects in the backscattered signal are only obtained if the
calculations are performed to second order [4, 13, 7].

Although the RTE method may in principle be simpler to implement than DBA-
based models, the disadvantage of the former is that it is not coherent: it
only yields the backscattered energy of the incident signal. Nowadays this is
a limitation since in the last decades full polarimetric orbital missions have
been developed, i.e. radars that can measure both the amplitude and the phase
of the received signals. Furthermore, they can measure in the different modes
of incident and received polarization (i.e. four channels: HH, VH, VV and HV).
For example, the RADARSAT-2 mission operates in C-band
($f=5.3\leavevmode\nobreak\ \text{GHz}$) and measures backscattered power on
all four channels [17]. The SAOCOM mission operates in L-band
($f=1.23\leavevmode\nobreak\ \text{GHz}$) and measures amplitude and phase on
all four channels [18]. Therefore, in order to make full use of the
information obtained by this type of mission, it is necessary to have direct
models that can account for all the observed parameters (namely, the
covariance matrix $C$ or the coherence matrix $T$ [19]). By making full use of
the information obtained from polarimetric missions, the inference of
biophysical parameters of interest, such as soil moisture, can be
substantially improved. For example, in [20] the dielectric constant of non
vegetated soils is estimated from the $T$ matrix using the Small Perturbation
Method in layered media as a direct model [21, 22]. Following this line of
work, the next step is to consider the vegetation on the surface. This leads
to have a theoretical model that allows to calculate the amplitude and phase
of the backscattered signal in such a scenario. This article is oriented
towards this goal.

Because of the same reasons that have guided our preference for an analytic
model, namely, that its use in a Bayesian scheme involves running the model a
large number of times, we have aimed for the simplest possible theoretical
framework. We have opted for a model based on field theory techniques, namely
we compute the scattered intensities and phases from a Schwinger-Dyson
equation where the characteristics of the vegetation layer are introduced
through its thickness and the two-point correlation in the dielectric constant
fluctuations. The complexity of the vegetation layer makes an accurate
modeling of these correlations elusive, and a truly realistic model would
necessitate a large number of parameters, thereby detracting from predictive
power. Our experience is that given the large uncertainties in relevant
parameters a qualitatively correct ansatz for the two point correlations, of
which we will give examples below, is all that is required for a successful
model. Also by reducing the statistical features of the dielectric constant
fluctuations to their two point correlations we are assuming these
fluctuations are Gaussian; for a more general treatment of non Gaussian
fluctuations see [23, 24].

Once the problem has been formulated as a set of Schwinger-Dyson equations, we
shall once again choose simplicity by seeking a solution through a second
order Born approximation [13, 14, 15, 16]. This is the most straigthforward
way to find fully analytical expressions for the scattered fields, whereby we
find the scattered energy by taking the second moments of the first and second
order fields. Actually at this order in powers of the dielectric constant
fluctuation we may find interference terms between the zeroeth and fourth, and
first and third, order solutions, but these are known to be negligible [13].
Given the order of magnitude of the scattering amplitudes, computing higher
orders of the perturbative expansion is unwarranted.

Still the Born approximation yields the desired intensities as multiple
integrals over the vegetation layer. For a simple enough form of the two-point
fluctuation correlation these integrals are elementary but their numerical
evaluation is nontrivial, not least because the integrands are strongly
oscillatory. However, this same feature makes a stationary phase approximation
accurate and much simpler to implement. As a check, we recover a factor two
enhancement of the backscattered flux, as expected for coherent scattering
[14, 25, 26].

In summary, we present in this contribution a fully coherent, analytic model
of the backscattering intensity in all HH, HV, VH and VV channels, which is
simple enough to allow repeated evaluation as a forward model in a Bayesian
inference scheme. This set up may be used to obtain relevant parameters such
as soil humidity, or simply as a fast way to identify the best model to be
employed as input in more sophisticated approaches [27, 28, 29, 30, 31], among
which we include the consideration of non Gaussian statistical fluctuations.
Moreover, the model only yields the volume scattering from the vegetation
layer, to which we must add the surface scattering. Given the intensities
involved, these two effects simply add and no cross correlation needs to be
considered.

This article is organized as follows. In Section 2 we established the system
of equations that must be solved to find the amplitude and phase of the wave
scattered by a stratified and inhomogeneous medium. Then, in Section 3 we
develop the perturbative scheme that we will use to solve the equations given
above. Section 4 and Section 5 are devoted to compute up to second order the
scattering amplitude when the incident field corresponds to Transverse
Electric (TE) or Magnetic (TM) case. In Section 6 we give the mean values of
the scattered energy at each order in perturbations and we discuss how to
extend these results to compute the polarimetric $T$ or $C$ matrix. In Section
7 we show the main results obtained from the proposed model. Finally, Section
8 contains our conclusions and future developments arising from what has been
show here.

##  2 Preliminaries

Our first objective is to establish the basic scheme that we will use later to
calculate, through a perturbative scheme, the fields scattered by a stratified
non-homogeneous target as we show in Figure 1.

Figure 1: Geometry of the problem: the dielectric permittivity is piece-wise
constant, with a discontinuity at $z=0$ and at $z=-h$.

In this configuration there are three regions of interest: region I, from
where the incident wave comes from; region II with flat interfaces at $z=0$
(with region I) and at $z=-d$ with region III; and region III, a semi-infinite
medium starting at $z=-d$. Regions I and III are homogeneous with dielectric
permittivity $\epsilon_{0I}$ and $\epsilon_{0III}$, respectively; the region
II has a non-homogeneous permittivity that can be decomposed as its mean value
plus fluctuation, i.e.,
$\epsilon_{II}=\epsilon_{0II}+\varepsilon_{II}(\mathbf{r})$.

###  2.1 Problem Definition

Our starting point is the Maxwell equations without sources in each region

$\mathbf{\Delta}\mathbf{E}-\nabla\left(\nabla\mathbf{E}\right)+\epsilon(\mathbf{r})\,\left(\frac{\omega}{c}\right)^{2}\mathbf{E}=0$
(1)

where the dielectric permittivity $\epsilon(\mathbf{r})=\epsilon(\vec{x},z)$
indicates a piece-wise medium. We look for solutions of the form

$\mathbf{E}(\mathbf{r})=\int\frac{d^{2}p}{(2\pi)^{2}}\,e^{i\vec{p}\vec{x}}\,\mathbf{E}_{\vec{p}}(z)$
(2)

where $\vec{x}=\left(x,y\right)$ and $\vec{p}=\left(p_{x},p_{y}\right)$. Also
we call $\mathbf{E}_{\vec{p}}=\left(E_{\vec{p}},E_{\vec{p}z}\right)$ where
further $E_{\vec{p}}=\left(E_{\vec{p}x},E_{\vec{p}y}\right)$. The transverse
component $\vec{E}_{\vec{p}}$ can be further decomposed into $H$ polarized
modes where
$\mathbf{E}_{H\vec{p}}=\left(\mathbf{K}\times\hat{p}\right)\varphi\left(z\right)$
and $V$ polarized modes where $\vec{E}_{V\vec{p}}=\hat{p}\psi\left(z\right)$,
where $\hat{p}=\vec{p}/p$.

For the $H$ modes we have $E_{H\vec{p}z}=\vec{p}\vec{E}_{H\vec{p}}=0$, and
then

$\frac{\partial^{2}}{\partial z^{2}}\varphi+P^{2}\varphi=0$ (3)

where $P^{2}=\epsilon(z)\left(\frac{\omega}{c}\right)^{2}-p^{2}$ is the
vertical component of the wave number in each space region.

For $V$ modes

$E_{z}=\frac{\imath p}{P^{2}}\frac{\partial\psi}{\partial z}$ (4)

and

$\frac{\partial^{2}\psi}{\partial z^{2}}+\frac{\partial}{\partial
z}\frac{p^{2}}{P^{2}}\frac{\partial\psi}{\partial
z}+\epsilon\left(\frac{\omega}{c}\right)^{2}\psi=0$ (5)

Observe that $E_{z}$ is defined from $\psi^{\prime}$ evaluated at the same
height.

Along with the dynamic equations for each mode we must give the matching
conditions for each of them at the interfaces where the dielectric
permittivity is not continuous. For the $H$ case we ask that $\varphi$ and
$\varphi^{\prime}$ are continuous over each interface; for the $V$ mode we
require that $\psi$ and $\epsilon E_{z}$ are continuous.

###  2.2 Coupled equations and cross-pol effect

We now consider the case where in each region of the space the dielectric
permittivity is written es
$\epsilon=\epsilon_{0}\left(z\right)+\varepsilon\left(\vec{x},z\right)$, with
$\epsilon_{0}\left(z\right)$ piece-wise constant. In the case of fluctuations
$\varepsilon$ we also perform a two-dimensional Fourier transform

$\displaystyle\varepsilon_{\mathbf{p}}(z)=\int\frac{d^{2}x}{(2\pi)^{2}}\,e^{\imath\mathbf{p}\cdot\mathbf{x}}\,\varepsilon(\mathbf{x},z)$
$\displaystyle\varepsilon^{*}_{\mathbf{p}}(z)=\varepsilon_{-\mathbf{p}}(z)\,\,,\,\,\langle\varepsilon_{\mathbf{p}}(z)\rangle=0$
(6)

The Fourier amplitudes are a Gaussian process with zero mean and correlation

$\displaystyle\langle\varepsilon(\mathbf{r})\,\varepsilon^{*}(\mathbf{r}^{\prime})\rangle=C(\mathbf{r}-\mathbf{r}^{\prime})$
$\displaystyle\langle\varepsilon_{\mathbf{p}}(z)\,\varepsilon^{*}_{\mathbf{q}}(z^{\prime})\rangle=\delta(\mathbf{p}-\mathbf{q})\,C_{\mathbf{p}}(z-z^{\prime})$
(7)

where $C_{\mathbf{p}}(z)$ is the 2D Fourier transform of correlation function
$C(\mathbf{r}-\mathbf{r}^{\prime})$. Again we decompose the Fourier expansion
of the fields in terms of $H$ and $V$ modes,
$\vec{E}_{\vec{p}}=\left(\mathbf{K}\times\hat{p}\right)\varphi_{\vec{p}}+\hat{p}\psi_{\vec{p}}$
and write these modes in terms of the orthogonal component of the electric
field, so

$\varphi_{\vec{p}}=\mathbf{K}\cdot\left(\hat{p}\times\vec{E}_{\vec{p}}\right)=-\hat{p}_{y}\,E_{\vec{p}x}+\hat{p}_{x}\,E_{\vec{p}y}$
(8a)
$\psi_{\vec{p}}=\hat{p}\cdot\vec{E}_{\vec{p}}=\hat{p}_{x}\,E_{\vec{p}x}+\hat{p}_{y}\,E_{\vec{p}y}$
(8b)

Using the above relations, we get a Klein-Gordon type equation for each field

$\displaystyle\frac{\partial^{2}}{\partial
z^{2}}\varphi_{\vec{p}}+P^{2}\varphi_{\vec{p}}\equiv-J_{XH\vec{p}}$
$\displaystyle
J_{XH\vec{p}}={\left(\frac{\omega}{c}\right)^{2}}\int\frac{d^{2}q}{\left(2\pi\right)^{2}}\epsilon_{1\left(\vec{p}-\hat{q}\right)}$
$\displaystyle\left[\hat{p}\cdot\hat{q}\;\varphi_{\vec{q}}+\mathbf{K}\cdot\left(\hat{p}\times\hat{q}\right)\psi_{\vec{q}}\right]$
(9a) $\displaystyle\frac{\partial^{2}}{\partial
z^{2}}\psi_{\vec{p}}+\frac{\partial}{\partial
z}\frac{p^{2}}{P^{2}}\frac{\partial}{\partial
z}\psi_{\vec{p}}+\epsilon_{0}\left(\frac{\omega}{c}\right)^{2}\psi_{\vec{p}}$
$\displaystyle\equiv-J_{XV\vec{p}}-\frac{\partial}{\partial
z}J^{\prime}_{XV\vec{p}}$ $\displaystyle
J_{XV\vec{p}}={\left(\frac{\omega}{c}\right)^{2}}\int\frac{d^{2}q}{\left(2\pi\right)^{2}}\epsilon_{1\left(\vec{p}-\vec{q}\right)}$
$\displaystyle\left[\mathbf{K}\cdot\left(\hat{p}\times\hat{q}\right)\varphi_{\vec{q}}-\hat{p}\cdot\hat{q}\;\psi_{\vec{q}}\right]$
$\displaystyle J^{\prime}_{XV\vec{p}}=\imath
p\left(\frac{\omega}{c}\right)^{2}\frac{1}{P^{2}}\int\frac{d^{2}q}{\left(2\pi\right)^{2}}\epsilon_{1\left(\vec{p}-\vec{q}\right)}E_{\vec{q}z}$
(9b)

Equation (9b) shows that it is possible to generate $V$ polarization from $H$
polarization, though not in backscattering to first order.

##  3 Perturbative expansion

As we discussed above, we consider the case where we have three regions. In
region $I$, $z>0$, $\epsilon_{0}=\epsilon_{0I}$. In region $II$, $0>z>-h$,
$\epsilon_{0}=\epsilon_{0II}$ . In region $III$, $-h>z$,
$\epsilon_{0}=\epsilon_{0III}$ . The dielectric constant is not fluctuating in
regions $I$ and $III$, and shows fluctuations $\varepsilon_{1}$ in region
$II$. We shall seek an expansion in powers of $\varepsilon_{1}$. To this
effect we write

$\displaystyle\varphi$ $\displaystyle=$
$\displaystyle\varphi^{\left(0\right)}+\varphi^{\left(s\right)}$
$\displaystyle\psi$ $\displaystyle=$
$\displaystyle\psi^{\left(0\right)}+\psi^{\left(s\right)}$ $\displaystyle
E_{z}$ $\displaystyle=$ $\displaystyle
E_{z}^{\left(0\right)}+E_{z}^{\left(s\right)}$ (10)

The zero superscript fields are solutions to the homogeneous problem with
horizontal wave number $\vec{k}$ which in region III are purely downgoing,
$\mathbf{E}=\mathbf{E}_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{x}}$. Therefore
the magnetic field (recall that $\left|\mathbf{k}\right|=\omega/c$) is
$\mathbf{B}=\hat{\mathbf{k}}\times\mathbf{E}$ and the Poynting vector results
$\mathbf{S}=\hat{\mathbf{k}}\left|\mathbf{E}\right|^{2}$.

We have two situations, the horizontal polarization case (TE mode) where

$\varphi^{\left(0\right)}_{H\vec{p}}=\left(2\pi\right)^{2}\delta\left(\vec{p}-\vec{k}\right)\sqrt{2K_{I}}F^{\downarrow}_{Hk}\,\,,\,\,\psi^{\left(0\right)}_{H\vec{p}}=0$
(11)

and the vertical polarization case (TM mode) where we have

$\varphi^{\left(0\right)}_{V\vec{p}}=0\,\,,\,\,\psi^{\left(0\right)}_{V\vec{p}}=\left(2\pi\right)^{2}\delta\left(\vec{p}-\vec{k}\right)\sqrt{2K_{I}}F^{\downarrow}_{Vk}$
(12)

In both modes $E^{\left(0\right)}_{zV\vec{p}}$ is obtained from equation (4).

For either $X=H,V$, the scattered fields $\varphi^{\left(s\right)}_{X}$
$\psi^{\left(s\right)}_{X}$ obey equations of the form of (9a) or (9b). Also,
the vertical component of the electric field satisfies,

$\displaystyle
P^{2}E^{\left(s\right)}_{X\vec{p}z}+\left(\frac{\omega}{c}\right)^{2}\int\frac{d^{2}q}{\left(2\pi\right)^{2}}\epsilon_{1\left(\vec{p}-\vec{q}\right)}E^{\left(s\right)}_{X\vec{q}z}=\imath
p\frac{d}{dz}\psi^{\left(s\right)}_{X\vec{p}}$
$\displaystyle-\left(\frac{\omega}{c}\right)^{2}\int\frac{d^{2}q}{\left(2\pi\right)^{2}}\epsilon_{1\left(\vec{p}-\vec{q}\right)}E^{\left(0\right)}_{X\vec{q}z}$
(13)

The sources $J_{XH\vec{p}}$ and $J_{XV\vec{p}}$ in the right-hand sides of
equations (9a) and (9b) are non vanishing in region $II$ only. The solutions
can be expressed in terms of the Green functions of the homogeneous problem

$G_{\left(H,V\right)\vec{p}}\left(z,z^{\prime}\right)=\imath\,F^{UP}_{\left(H,V\right)p}\left(z_{>}\right)F^{DOWN}_{\left(H,V\right)p}\left(z_{<}\right)$
(14)

The $UP\,(\uparrow)$ or $DOWN(\downarrow)$ modes are purely outgoing solutions
normalized to unit Wronskian. For $H$ and $V$ modes the Wronskian of two
solutions is

$\displaystyle
W_{H1,2}=\left(-i\right)\left[\varphi_{2}\varphi^{\prime}_{1}-\varphi_{1}\varphi^{\prime}_{2}\right]$
(15) $\displaystyle
W_{V1,2}=\left(-i\right)\left[1+\frac{p^{2}}{P^{2}}\right]\left\\{\psi^{\prime}_{1}\psi_{2}-\psi_{1}\psi^{\prime}_{2}\right\\}$
(16)

For both modes $m=H\,\text{or}\,V$ the outgoing waves, in each region of the
space, can be written as follows,

$\begin{cases}F^{\uparrow}_{m\,pI}=f_{m}(P_{I})\,e^{\imath iP_{I}z}\\\
F^{\uparrow}_{m\,pII}=f_{m}(P_{II})\,\left[r^{\uparrow}_{m\,pII}e^{\imath
P_{II}z}+t^{\uparrow}_{m\,pII}e^{-\imath P_{II}z}\right]\\\
F^{\uparrow}_{m\,pIII}=f_{m}(P_{III})\,\left[\frac{1}{t^{\downarrow}_{m\,pIII}}e^{\imath
P_{III}z}+t^{\uparrow}_{m\,pIII}e^{-\imath P_{III}z}\right]\end{cases}$ (17)

$\begin{cases}F^{\downarrow}_{m\,pI}=f_{m}(P_{I})\,\left[r^{\downarrow}_{m\,pI}e^{\imath
P_{I}z}+e^{-\imath P_{I}z}\right]\\\
F^{\downarrow}_{m\,p_{II}}=f_{m}(P_{II})\,\left[r^{\downarrow}_{m\,pII}e^{\imath
P_{II}z}+t^{\downarrow}_{m\,pII}e^{-\imath P_{II}z}\right]\\\
F^{\downarrow}_{m\,p_{III}}=f_{m}(P_{III})\,t^{\downarrow}_{m\,pIII}e^{-\imath
P_{III}z}\end{cases}$ (18)

where the coefficients $r$ and $t$ for each mode are the Fresnel coefficients
of reflection and transmission, and
$f_{m}(P_{\alpha})=\left\\{1/\sqrt{2P_{\alpha}},\sqrt{P_{\alpha}/(2\epsilon_{\alpha}(\omega/c)^{2})}\right\\}$
is a normalization factor enforcing the unit Wronskian condition. Note that
the Green functions obey the reciprocity condition [32],

$G_{\left(H,V\right)\vec{p}}\left(z,z^{\prime}\right)=G_{\left(H,V\right)\left(-\vec{p}\right)}\left(z^{\prime},z\right)$
(19)

which, since by rotation symmetry
$G_{\left(H,V\right)\vec{p}}\left(z,z^{\prime}\right)=G_{\left(H,V\right){p}}\left(z,z^{\prime}\right)$,
reduces to symmetry under the exchange of $z$ and $z^{\prime}$.

The solutions for the scattered fields in region $I$ are

$\varphi^{\left(s\right)}_{X\vec{p}}\left(z\right)=-F^{\uparrow}_{HpI}\left(z\right)\int_{-h}^{0}dz^{\prime}\;F^{\downarrow}_{HpII}\left(z^{\prime}\right)J_{XH\vec{p}}\left(z^{\prime}\right)$
(20a)
$\displaystyle\psi^{\left(s\right)}_{X\vec{p}}\left(z\right)=-F^{\uparrow}_{VpI}\left(z\right)\int_{-h}^{0}dz^{\prime}\;F^{\downarrow}_{VpII}\left(z^{\prime}\right)J_{XV\vec{p}}\left(z^{\prime}\right)$
$\displaystyle+F^{\uparrow}_{VpI}\left(z\right)\int_{-h}^{0}dz^{\prime}\;\frac{\partial
F^{\downarrow}_{VpII}}{\partial
z^{\prime}}\left(z^{\prime}\right)J^{\prime}_{XV\vec{p}}\left(z^{\prime}\right)$
(20b)

In the last term we have made an integration by parts. Since we are assuming
there are no fluctuations at the boundary of region II the integrated term
vanishes.

In region $II$ the fields result

$\displaystyle\varphi^{\left(s\right)}_{X\vec{p}}\left(z\right)=$
$\displaystyle-F^{\uparrow}_{HpII}\left(z\right)\int_{-h}^{z}\,dz^{\prime}\;F^{\downarrow}_{HpII}\left(z^{\prime}\right)J_{XH\vec{p}}\left(z^{\prime}\right)\,\Theta(z>z^{\prime})$
$\displaystyle-F^{\downarrow}_{HpII}\left(z\right)\int_{z}^{0}dz^{\prime}\;F^{\uparrow}_{HpII}\left(z^{\prime}\right)J_{XH\vec{p}}\left(z^{\prime}\right)\,\Theta(z<z^{\prime})$
$\displaystyle\psi^{\left(s\right)}_{X\vec{p}}\left(z\right)=$
$\displaystyle-F^{\uparrow}_{VpII}\left(z\right)\int_{-h}^{z}\,dz^{\prime}\;F^{\downarrow}_{VpII}\left(z^{\prime}\right)J_{XV\vec{p}}\left(z^{\prime}\right)\,\Theta(z>z^{\prime})$
$\displaystyle+F^{\uparrow}_{VpII}\left(z\right)\int_{-h}^{z}\,dz^{\prime}\;\frac{\partial
F^{\downarrow}_{VpII}}{\partial
z^{\prime}}J^{\prime}_{XV\vec{p}}\left(z^{\prime}\right)\,\Theta(z>z^{\prime})$
$\displaystyle-F^{\downarrow}_{VpII}\left(z\right)\int_{-h}^{0}\,dz^{\prime}\;F^{\uparrow}_{VpII}\left(z^{\prime}\right)J_{XV\vec{p}}\left(z^{\prime}\right)\,\Theta(z<z^{\prime})$
$\displaystyle+F^{\downarrow}_{VpII}\left(z\right)\int_{-h}^{0}\,dz^{\prime}\;\frac{\partial
F^{\downarrow}_{VpII}}{\partial
z^{\prime}}\left(z^{\prime}\right)J^{\prime}_{XV\vec{p}}\left(z^{\prime}\right)\,\Theta(z<z^{\prime})$

##  4 TE incident mode

To study the TE incident mode we use as source field equations (11). We will
define two scattering amplitudes: $S_{HH}$ and $S_{VH}$ for the $H$-mode and
$V$-mode scattered fields, respectively. Each of these amplitudes will be
computed order by order in terms of the perturbations in the dielectric
permittivity of region II as well as the lower order $\varphi$ and $\psi$
fields. Namely, the solution to first order will be proportional to the
product between the zero-order fields and the dielectric fluctuations, the
second order solution will be proportional to the first-order fields and the
dielectric fluctuations, and so on.

In order to find the scattered field in $H$ mode at first order we use
equation (20a) evaluating the source with the field given by (11). Similarly,
the $V$ mode is obtained using (20) with the incident field (11). The first-
order results evaluated in region $II$, which are given by equations (21) and
(21), will be used as sources to find the scattered fields to second order.

To simplify the upcoming formulae, we introduce the following notations. An
index
$\mu=\\{\uparrow,\downarrow\\}(\Rightarrow\bar{\mu}=\\{\downarrow,\uparrow\\})$
indicates the propagation mode of waves in region II,
$\Theta_{\uparrow}(z,z^{\prime})=\Theta(z-z^{\prime})(\Rightarrow\Theta_{\downarrow}(z,z^{\prime})=\Theta(z^{\prime}-z))$
shows where the second scattering process occurs. Index $m=\\{H,V\\}$
indicates the polarization of waves in medium II; associated with this index
we have the geometric coefficients
$\beta_{H}=(\hat{p}\cdot\hat{q})\,(\hat{q}\cdot\hat{k})$ and
$\beta_{V}=\mathbf{K}\cdot\left(\hat{p}\times\hat{q}\right)\,\mathbf{K}\cdot\left(\hat{q}\times\hat{k}\right)$.
We also define

$\frac{dF^{\mu}_{mQ}}{dz}(z)=F^{\prime\mu}_{mQ}(z)=\imath
Q\tilde{F}^{\mu}_{mQ}(z)$ (22)

and coefficients
$\gamma_{H}=\mathbf{K}\cdot\left(\hat{p}\times\hat{q}\right)\,(\hat{q}\cdot\hat{k})$,
$\gamma_{V}=-\left(\hat{p}\cdot\hat{q}\right)\,\mathbf{K}\cdot(\hat{q}\times\hat{k})$.
We obtain

$\displaystyle\varphi_{\mathbf{p}}^{(2)}(z)=-\sqrt{2K_{I}}\,\left(\frac{\omega}{c}\right)^{4}\,F^{\uparrow}_{HPI}(z)$
$\displaystyle\int\,\frac{d^{2}q}{(2\pi)^{2}}\,\int^{0}_{-h}\,dz_{1}\,\int^{L}_{l}\,dz_{2}\,\,\varepsilon_{\mathbf{p}-\mathbf{q}}(z_{1})\,\varepsilon_{\mathbf{q}-\mathbf{k}_{i}}(z_{2})$
$\displaystyle
F^{\downarrow}_{HPII}(z_{1})\,\left[\sum_{\mu\,m}\beta_{m}\,F^{\mu}_{mQII}(z_{1})\right.$
$\displaystyle\left.F^{\hat{\mu}}_{mQII}(z_{2})\,\Theta_{\mu}(z_{1},z_{2})\right]\,F^{\downarrow}_{HKII}(z_{2})$
(23)

$\displaystyle\psi^{(2)}_{\mathbf{p}}(z)=-\sqrt{2KI}\left(\frac{\omega}{c}\right)^{4}\,F^{\uparrow}_{VPI}(z)$
$\displaystyle\int^{0}_{-h}\,dz_{1}\,\int^{0}_{-h}\,dz_{2}\,\int\frac{d^{2}q}{(2\pi)^{2}}\varepsilon_{\mathbf{p}-\mathbf{q}}(z_{1})\,\varepsilon_{\mathbf{q}-\mathbf{k}_{i}}(z_{2})$
$\displaystyle
F^{\downarrow}_{VPII}(z_{1})\left[\sum_{\mu\,m}\gamma_{m}F^{\mu}_{mQII}(z_{1})\right.$
$\displaystyle\left.F^{\bar{\mu}}_{mQII}(z_{2})\Theta_{\mu}(z_{1},z_{2})\right]F^{D}_{HKII}(z_{2})$
(24)

where the index $m$ now takes three values $\left\\{H,V,\tilde{V}1\right\\}$.
The last value indicates that the two functions $F_{V}$ evaluated at $z_{1}$
have a tilde.

##  5 TM incident mode

To study the TM incident mode we shall proceed in the same way than we did in
Section 4. Now we use the as source the field given by equations (12). We will
define two new scattering amplitudes: $S_{HV}$ and $S_{VV}$, which correspond
to $H$ mode and $V$ mode scattered fields, respectively.

We use the incident field given by (12) to evaluate the currents. To obtain
the scattering amplitude for the HV channel to first order we use (20a). We
then use the first-order fields as sources to write the propagating waves in
region II with equations (21) and (21) and compute the currents. We define the
coefficients

$\displaystyle\left\\{\delta_{H};\delta_{V},\delta_{\tilde{V}2}\right\\}=\left\\{(\hat{p}\cdot\hat{q})\,\mathbf{K}\cdot(\hat{q}\times\hat{k}),\right.$
$\displaystyle\left.\mathbf{K}\cdot(\hat{p}\times\hat{q})\,(\hat{q}\cdot\hat{k}),-\mathbf{K}\cdot\left(\hat{p}\times\hat{q}\right)\,\frac{q\,k}{QII\,KII}\right\\}$

corresponding to the index $m=\left\\{H,V,\tilde{V}2\right\\}$. The amplitudes
for the scattered fields result

$\displaystyle\varphi_{\mathbf{p}}^{(2)}=-\sqrt{2K_{I}}\,\left(\frac{\omega}{c}\right)^{4}\,F^{\uparrow}_{HPI}(z)$
$\displaystyle\int\,\frac{d^{2}q}{(2\pi)^{2}}\,\int^{0}_{-h}\,dz_{1}\,\int^{L)}_{l)}\,dz_{2}\,\varepsilon_{\mathbf{p}-\mathbf{q}}(z_{1})\,\varepsilon_{\mathbf{q}-\mathbf{k}_{i}}(z_{2})$
$\displaystyle
F^{\downarrow}_{HPII}(z_{1})\left[\sum_{\mu\,m}\,\delta_{m}\,F^{\mu}_{mQII}(z_{1})\right.$
$\displaystyle\left.F^{\hat{\mu}}_{mQII}(z_{2})\,F^{\downarrow}_{VKII}(z_{2})\,\Theta_{\mu}(z_{1},z_{2})\right]$
(25)

$\displaystyle\psi^{(2)}_{\mathbf{p}}(z)=\sqrt{2K_{I}}\,\left(\frac{\omega}{c}\right)^{4}\,F^{\uparrow}_{VPI}(z)$
$\displaystyle\int\,\frac{d^{2}q}{(2\pi)^{2}}\,\int^{0}_{-h}\,dz_{1}\,\int^{L}_{l}\,dz_{2}\,\,\varepsilon_{\mathbf{p}-\mathbf{q}}(z_{1})\,\varepsilon_{\mathbf{q}-\mathbf{k}_{i}}(z_{2})$
$\displaystyle\sum_{\mu\,m}\,\eta_{m}\,F^{\downarrow}_{VPII}(z_{1})\,F^{\mu}_{mQII}(z_{1})$
$\displaystyle
F^{\hat{\mu}}_{mQII}(z_{2})\,F^{\downarrow}_{VKII}(z_{2})\,\Theta_{\mu}(z_{1},z_{2})$
(26)

Here the index $m$ takes five values:
$\left\\{H,V,\tilde{V}_{1},\tilde{V}_{2},\tilde{V}_{12}\right\\}$; the last
index means that the four functions $F$ with polarization $V$ have tilde. The
geometric coefficients are:

$\displaystyle\left\\{\eta_{H},\eta_{V},\eta_{\tilde{V}1},\eta_{\tilde{V}2},\eta_{\tilde{V}12}\right\\}=$
$\displaystyle\left\\{-\mathbf{K}\cdot(\hat{p}\times\hat{q})\,\mathbf{K}\cdot(\hat{q}\times\hat{k}),(\hat{p}\cdot\hat{q})\,(\hat{q}\cdot\hat{k}),\right.$
$\displaystyle\left.-\frac{(\hat{p}\cdot\hat{q})\,p\,k}{PII\,KII},\frac{(\hat{q}\cdot\hat{k})\,p\,q}{PII\,QII},-\frac{p\,q^{2}\,k}{PII\,Q^{2}II\,KII}\right\\}.$

##  6 Mean values

###  6.1 First Order

The scattering amplitude for HH, VH and HV channels can be written as

$\displaystyle
S^{(1)}_{ab}(z)=\sqrt{2KI}\,\left(\frac{\omega}{c}\right)^{2}\,\mathcal{A}_{ab}(PI,PII,KII)$
$\displaystyle
F^{\uparrow}_{aPI}(z)\,\int^{0}_{-h}\,dz^{\prime}\,\varepsilon_{\mathbf{p}-\mathbf{k}}(z^{\prime})\,\left[r_{a}\,e^{\imath
PIIz^{\prime}}+t_{a}\,e^{-\imath PIIz^{\prime}}\right]$
$\displaystyle\left[r_{b}\,e^{\imath KIIz^{\prime}}+t_{b}\,e^{-\imath
KIIz^{\prime}}\right]$ (27)

where $a\,(b)$ is the polarization of the scattered (incident) wave,
$\mathcal{A}_{ab}$ is a factor coming from the product between the $F$
functions involved in the scattering amplitude and $r$ ($t$) is the reflection
(transmission) coefficient of downward waves with polarization $a$ or $b$
evaluated at the $PII$ or $KII$ mode, respectively. From now on we will
replace $\omega/c$ by $k_{0}$, the incident wave number in region I.

The mean scattered fields vanish in virtue of (2.2). For the scattered mean
energy we must compute:

$\displaystyle\langle\big{|}S^{(1)}_{ab}(z)\big{|}^{2}\rangle=\delta_{\mathbf{k}}(\mathbf{0})\,2KI\,k^{4}_{0}\,|\mathcal{A}_{ab}|^{2}\,|F^{\uparrow}_{aPI}|^{2}$
$\displaystyle\int^{0}_{-h}\,dz_{1}\,\int^{0}_{-h}\,dz_{2}\,C_{\mathbf{p}-\mathbf{k}}(z_{1}-z_{2})$
$\displaystyle\left[|r_{a}|^{2}\,e^{\imath(z_{1}-z_{2})PII}+|t_{a}|^{2}\,e^{-\imath(z_{1}-z_{2})PII}\right.$
$\displaystyle\left.+2\Re\left\\{r_{a}\,t^{*}_{a}e^{\imath(z_{1}+z_{2})PII}\right\\}\right]$
$\displaystyle\left[|r_{b}|^{2}\,e^{\imath(z_{1}-z_{2})KII}+|t_{b}|^{2}\,e^{-\imath(z_{1}-z_{2})KII}\right.$
$\displaystyle\left.+2\Re\left\\{r_{b}\,t^{*}_{b}e^{\imath(z_{1}+z_{2})KII}\right\\}\right]$
(28)

Here $\delta_{\mathbf{k}}(\mathbf{0})$ indicates the area illuminated by the
incident wave from region I.

Finally, for VV case we must note from (5) there is an extra term coming from
the term proportional to $E_{z}$ given by (4). Besides that, the mean energy
can be written similarly to (6.1) with the difference that, instead of making
the product of two functions $F$ evaluated at $z_{1}$ by two functions $F$
evaluated at $z_{2}$, we must make the product between
$\mathcal{A}_{ab}\,F_{a}\,F_{b}+\tilde{\mathcal{A}}_{ab}\tilde{F}_{a}\,\tilde{F}_{b}$
at $z_{1}$ and the corresponding one evaluated at $z_{2}$. This will produce
an expression analogous to (6.1) but with more terms in the double integral
and without the global factor $|\mathcal{A}_{ab}|^{2}$ but with terms
dependent on $|\mathcal{A}_{ab}|^{2}$, $|\tilde{\mathcal{A}}_{ab}|^{2}$,
$\mathcal{A}_{ab}\,\tilde{\mathcal{A}}^{*}_{ab}$ and their complex conjugate.

###  6.2 Second Order

The four channels can be arranged as

$\displaystyle
S^{(2)}_{ab}(z)=-\sqrt{2K_{I}}\,k^{4}_{0}\,F^{\uparrow}_{aPI}(z)$
$\displaystyle\int\,\frac{d^{2}q}{(2\pi)^{2}}\,\int^{0}_{-h}\,dz_{1}\,\int^{L}_{l}\,dz_{2}\,\,\varepsilon_{\mathbf{p}-\mathbf{q}}(z_{1})\,\varepsilon_{\mathbf{q}-\mathbf{k}_{i}}(z_{2})$
$\displaystyle\sum_{\mu,m}\,F^{\downarrow}_{aPII}(z_{1})\,\mathcal{H}^{\mu,m}_{\mathbf{p},\mathbf{q}}(z_{1},z_{2})\,F^{\downarrow}_{bKII}(z_{2})$
(29)

where the function $\mathcal{H}^{\mu,m}_{\mathbf{p},\mathbf{q}}(z_{1},z_{2})$
contains all dependence on geometrical factors and on the products between $F$
and/or $\tilde{F}$ functions as well as the Heaviside functions for each
propagation mode.

At second order, the scattering amplitudes depend on the product between
permittivity fluctuations evaluated at different points of region II,
representing multiple scattering processes. Therefore, the mean field in no
longer zero, but the mean energy is non vanishing in the specular direction
only.

To compute the average energy to second order we start from

$\displaystyle\langle|S^{(2)}_{ab}(z)|^{2}\rangle=2K_{I}\,k^{8}_{0}\,|F^{\uparrow}_{aPI}(z)|^{2}$
$\displaystyle\int\,\frac{d^{2}q}{(2\pi)^{2}}\int\frac{d^{2}q^{\prime}}{(2\pi)^{2}}\int^{0}_{-h}dz_{1}\int^{L}_{l}dz_{2}\int^{0}_{-h}dz^{\prime}_{1}\int^{L^{\prime}}_{l^{\prime}}dz^{\prime}_{2}$
$\displaystyle\langle\varepsilon_{\mathbf{p}-\mathbf{q}}(z_{1})\,\varepsilon_{\mathbf{q}-\mathbf{k}_{i}}(z_{2})\,\epsilon^{*}_{\mathbf{p}-\mathbf{q}^{\prime}}(z^{\prime}_{1})\,\epsilon^{*}_{\mathbf{q}^{\prime}-\mathbf{k}_{i}}(z^{\prime}_{2})\rangle$
$\displaystyle\left[\sum_{\mu,m}\,F^{\downarrow}_{aPII}(z_{1})\,\mathcal{H}^{\mu,m}_{\mathbf{p},\mathbf{q}}(z_{1},z_{2})\,F^{\downarrow}_{bKII}(z_{2})\right]$
$\displaystyle\left[\sum_{\mu^{\prime},m^{\prime}}\,F^{\downarrow}_{aPII}(z^{\prime}_{1})\,\mathcal{H}^{\mu^{\prime},m^{\prime}}_{\mathbf{p},\mathbf{q}^{\prime}}(z^{\prime}_{1},z^{\prime}_{2})\,F^{\downarrow}_{bKII}(z^{\prime}_{2})\right]^{*}$
(30)

We apply Wick theorem in order to compute the expectation value in (6.2). The
mean energy at second order can be rearranges as follows,

$\displaystyle\langle|S^{(2)}_{ab}(z)|^{2}\rangle=\big{|}\langle
S^{(2)}_{ab}(z)\rangle\big{|}^{2}+2K_{I}\,k^{8}_{0}\,|F^{\uparrow}_{aPI}(z)|^{2}$
$\displaystyle\int\frac{d^{2}q}{(2\pi)^{2}}\int\frac{d^{2}q^{\prime}}{(2\pi)^{2}}\int^{0}_{-h}dz_{1}\int^{L}_{l}dz_{2}\int^{0}_{-h}dz^{\prime}_{1}\int^{L^{\prime}}_{l^{\prime}}dz^{\prime}_{2}$
$\displaystyle\sum_{\mu,m}\,\,\sum_{\mu^{\prime},m^{\prime}}$
$\displaystyle\left[\left(\delta\left(\mathbf{q}-\mathbf{q}^{\prime}\right)\right)^{2}\,C_{\mathbf{p}-\mathbf{q}}(z_{1}-z^{\prime}_{1})\,C_{\mathbf{q}-\mathbf{k}_{i}}(z_{2}-z^{\prime}_{2})\right.$
$\displaystyle\left.+\left(\delta\left(\mathbf{p}+\mathbf{k}_{i}-\mathbf{q}-\mathbf{q}^{\prime}\right)\right)^{2}\right.$
$\displaystyle\left.C_{\mathbf{p}-\mathbf{q}}(z_{1}-z^{\prime}_{2})\,C_{\mathbf{q}-\mathbf{k}_{i}}(z_{2}-z^{\prime}_{1})\right]$
$\displaystyle\left[F^{D}_{aPII}(z_{1})\,\mathcal{H}^{\mu,m}_{\mathbf{p},\mathbf{q}}(z_{1},z_{2})\,F^{D}_{bKII}(z_{2})\right]$
$\displaystyle\left[F^{D}_{aPII}(z^{\prime}_{1})\,\mathcal{H}^{\mu^{\prime},m^{\prime}}_{\mathbf{p},\mathbf{q}^{\prime}}(z^{\prime}_{1},z^{\prime}_{2})\,F^{D}_{bKII}(z^{\prime}_{2})\right]^{*}$
(31)

The incoherent scattered field is

$\displaystyle\langle|S^{(2)}_{ab}(z)|^{2}\rangle-\big{|}\langle
S^{(2)}_{ab}(z)\rangle\big{|}^{2}=\delta_{\mathbf{k}}(\mathbf{0})\,2K_{I}\,k^{8}_{0}|F^{\uparrow}_{aPI}(z)|^{2}$
$\displaystyle\int\frac{d^{2}q}{(2\pi)^{2}}\int^{0}_{-h}dz_{1}\int^{L}_{l}dz_{2}\int^{0}_{-h}dz^{\prime}_{1}\int^{L^{\prime}}_{l^{\prime}}dz^{\prime}_{2}$
$\displaystyle\left[\mathcal{L}_{ab}(z1-z^{\prime}_{1},z_{2}-z^{\prime}_{2})+\mathcal{C}_{ab}(z_{1}-z^{\prime}_{2},z_{2}-z^{\prime}_{1})\right]$
(32)

In the above expression we have defined the ladder ($\mathcal{L}$) and cross
($\mathcal{C}$) terms, listed below.

$\displaystyle\mathcal{L}_{ab}(z_{1}-z^{\prime}_{1},z_{2}-z^{\prime}_{2})=$
$\displaystyle=C_{\mathbf{p}-\mathbf{q}}(z_{1}-z^{\prime}_{1})\,C_{\mathbf{q}-\mathbf{k}_{i}}(z_{2}-z^{\prime}_{2})$
$\displaystyle\sum_{\mu,m}\,\sum_{\mu^{\prime},m^{\prime}}\,\left[F^{\downarrow}_{aPII}(z_{1})\,\mathcal{H}^{\mu,m}_{\mathbf{p},\mathbf{q}}(z_{1},z_{2})\,F^{\downarrow}_{bKII}(z_{2})\right]$
$\displaystyle\left[F^{\downarrow}_{aPII}(z^{\prime}_{1})\,\mathcal{H}^{\mu^{\prime},m^{\prime}}_{\mathbf{p},\mathbf{q}}(z^{\prime}_{1},z^{\prime}_{2})\,F^{\downarrow}_{bKII}(z^{\prime}_{2})\right]^{*}$
(33a)
$\displaystyle\mathcal{C}_{ab}(z_{1}-z^{\prime}_{2},z_{2}-z^{\prime}_{1})=$
$\displaystyle=C_{\mathbf{p}-\mathbf{q}}(z_{1}-z^{\prime}_{2})\,C_{\mathbf{q}-\mathbf{k}_{i}}(z_{2}-z^{\prime}_{1})$
$\displaystyle\sum_{\mu,m}\,\sum_{\mu^{\prime},m^{\prime}}\,\left[F^{\downarrow}_{aPII}(z_{1})\,\mathcal{H}^{\mu,m}_{\mathbf{p},\mathbf{q}}(z_{1},z_{2})\,F^{\downarrow}_{bKII}(z_{2})\right]$
$\displaystyle\left[F^{\downarrow}_{aPII}(z^{\prime}_{1})\,\mathcal{H}^{\mu^{\prime},m^{\prime}}_{\mathbf{p},\mathbf{r}}(z^{\prime}_{1},z^{\prime}_{2})\,F^{\downarrow}_{bKII}(z^{\prime}_{2})\right]^{*}$
(33b)

where $\mathbf{r}=\mathbf{p}+\mathbf{k}-\mathbf{q}$ (see Fig. (2)).

Figure 2: Ladder and cross diagrams.

With the expression (6.2) we complete, at least formally, the calculation of
the mean energy scattered to the region I. This expression could be evaluated
numerically, at a high computational cost. In the following, we shall present
an approximation which is nevertheless accurate enough for comparison against
field data.

To this end we will use a stationary phase approximation to the integrals over
oscillatory functions, reducing the problem to a single integral over the
intermediate impulse variable $\mathbf{q}$.

####  6.2.1 Stationary phase approximation

To get a compact expression for ladder a cross terms we introduce the
following notation: each function $F$ propagating with mode $\mu$ through
region $II$, has polarization $m$ and depends of mode $Q$, is writen as:

$\displaystyle F^{\mu}_{m\,Q}(z)=f_{m}(Q)\,\left[r^{\mu}_{m}(Q)\,e^{\imath
Qz}+t^{\mu}_{m}(Q)\,e^{-\imath Qz}\right]$ $\displaystyle
F^{\mu}_{m\,Q}(z)=f_{m}(Q)\,\sum_{s=\pm 1}\,X^{\mu s}_{m\,Q}\,e^{\imath sQz}$
$\displaystyle\tilde{F}^{\mu}_{m\,Q}(z)=f_{m}(Q)\,\sum_{s=\pm 1}\,s\,X^{\mu
s}_{m\,Q}\,e^{\imath sQz}$ $\displaystyle
X^{\mu\,+1}_{m\,Q}=r^{\mu}_{m}(Q)\,\,,\,\,X^{\mu\,-1}_{m\,Q}=t^{\mu}_{m}(Q)$
(34)

We shall carry on the calculation in detail for the HH channel only, where
indexes $m$ and $m^{\prime}$ just take values $H$ or $V$. However, the
following development is valid for the other three channels, provided the
possible values of these indexes are extended as appropriate for each case.

Using (34) the ladder term given by (40) becomes

$\displaystyle\mathcal{L}_{ab}=2K_{I}\,k^{8}_{0}\,|F^{\uparrow}_{a\,PI}(z)|^{2}$
$\displaystyle\int\,\frac{d^{2}q}{(2\pi)^{2}}\,\int^{0}_{-h}\,dz_{1}\,\int^{L}_{l}dz_{2}\,\int^{0}_{-h}dz^{\prime}_{1}\,\int^{L^{\prime}}_{l^{\prime}}dz^{\prime}_{2}$
$\displaystyle\sum_{\mu,m,s_{\mu},s_{\bar{\mu}}}\sum_{s_{1}\,s_{2}}\,\sum_{\mu^{\prime},m^{\prime},s^{\prime}_{\mu},s^{\prime}_{\bar{\mu}}}\sum_{s^{\prime}_{1}\,s^{\prime}_{2}}C_{\mathbf{p}-\mathbf{q}}(z_{1}-z^{\prime}_{1})$
$\displaystyle
C_{\mathbf{q}-\mathbf{k}_{i}}(z_{2}-z^{\prime}_{2})\Theta\left(s_{\mu}(z_{1}-z_{2})\right)\,\Theta\left(s^{\prime}_{\mu}(z^{\prime}_{1}-z^{\prime}_{2})\right)$
$\displaystyle\beta_{m}\,f_{a}(PII)\,f_{m}(QII)\,f_{m}(QII)\,f_{b}(KII)$
$\displaystyle\left[\beta_{m^{\prime}}\,f_{a}(PII)\,f_{m^{\prime}}(QII)\,f_{m^{\prime}}(QII)\,f_{b}(KII)\right]^{*}$
$\displaystyle
X^{D\,s_{1}}_{a\,PII}\,X^{\mu\,s_{\mu}}_{m\,QII}\,X^{\bar{\mu}\,s_{\bar{\mu}}}_{m\,QII}\,X^{D\,s_{2}}_{b\,KII}$
$\displaystyle\left[X^{D\,s^{\prime}_{1}}_{a\,PII}\,X^{\mu^{\prime}\,s^{\prime}_{\mu}}_{m^{\prime}\,QII}\,X^{\bar{\mu}^{\prime}\,s^{\prime}_{\bar{\mu}}}_{m^{\prime}\,QII}\,X^{D\,s^{\prime}_{2}}_{b\,KII}\right]^{*}$
$\displaystyle
e^{\imath(s_{1}\,PII+s_{\mu}\,QII)z_{1}}\,e^{\imath(s_{\bar{\mu}}\,QII+s_{2}\,KII)z_{2}}$
$\displaystyle
e^{-\imath(s^{\prime}_{1}\,PII+s^{\prime}_{\mu}\,QII)z^{\prime}_{1}}\,e^{-\imath(s^{\prime}_{\bar{\mu}}\,QII+s^{\prime}_{2}\,KII)z^{\prime}_{2}}$
(35)

We look for the stationary phase points: to eliminate the oscillation on
$z_{1}$ we need $s_{\mu}=-s_{1}$ (then $PII=QII$) and to eliminate the
oscillation on $z_{2}$ we need $s_{2}=-s_{\bar{\mu}}$ (then $QII=KII$);
analogous conditions arise on the primed variables. We use the above
assumptions to simplify the sums over the indices
$\left\\{s_{\mu},s_{\bar{\mu}},s^{\prime}_{\mu},s^{\prime}_{\bar{\mu}}\right\\}$.
After this we make the following change of variables:
$z_{j}=\xi_{j}+\frac{1}{2}\psi_{j}$,
$z^{\prime}_{j}=\xi_{j}-\frac{1}{2}\psi_{j}$. Then, we get

$\displaystyle\mathcal{L}_{ab}=2K_{I}\,k^{8}_{0}\,|F^{\uparrow}_{a\,PI}(z)|^{2}$
$\displaystyle\int\,\frac{d^{2}q}{(2\pi)^{2}}\,\int^{0}_{-h}\,d\xi_{1}\,\int^{L}_{l}\,d\psi_{1}\,\int\,d\xi_{2}\,\int\,d\psi_{2}$
$\displaystyle\sum_{s_{1}\,s_{2}}\,\sum_{\mu,m}\,\,\sum_{\mu^{\prime},m^{\prime}}\,\sum_{s^{\prime}_{1}\,s^{\prime}_{2}}\,\,C_{\mathbf{p}-\mathbf{q}}(\psi_{1})\,C_{\mathbf{q}-\mathbf{k}_{i}}(\psi_{2})$
$\displaystyle\Theta\left(s_{1}(\xi_{1}-\xi_{2}+\frac{1}{2}(\psi_{1}-\psi_{2}))\right)$
$\displaystyle\Theta\left(s^{\prime}_{1}(\xi_{1}-\xi_{2}-\frac{1}{2}(\psi_{1}-\psi_{2}))\right)$
$\displaystyle\beta_{m}\,\beta^{*}_{m^{\prime}}\,\,|f_{a}(PII)|^{2}\,|f_{m}(QII)|^{2}\,|f_{m^{\prime}}(QII)|^{2}$
$\displaystyle
X^{D\,s_{1}}_{a\,PII}\,X^{\mu\,-s_{1}}_{m\,QII}\,X^{\bar{\mu}\,-s_{2}}_{m\,QII}\,X^{D\,s_{2}}_{b\,KII}$
$\displaystyle\left[X^{D\,s^{\prime}_{1}}_{a\,PII}\,X^{\mu^{\prime}\,-s^{\prime}_{1}}_{m^{\prime}\,QII}\,X^{\bar{\mu}^{\prime}\,-s^{\prime}_{2}}_{m^{\prime}\,QII}\,X^{D\,s^{\prime}_{2}}_{b\,KII}\right]^{*}\,|f_{b}(KII)|^{2}$
$\displaystyle e^{\imath(PII-
QII)(s_{1}-s^{\prime}_{1})\xi_{1}}\,e^{\imath\frac{1}{2}(PII-
QII)(s_{1}+s^{\prime}_{1})\psi_{1}}$ $\displaystyle
e^{\imath(QII-\,KII)(s_{2}-s^{\prime}_{2})\xi_{2}}\,e^{\imath\frac{1}{2}(QII-\,KII)(s_{2}+s^{\prime}_{2})\psi_{2}}$
(36)

Here it is worth noting that:

  * •

integration over $\xi_{j}$ takes its maximum value when
$s_{j}=s^{\prime}_{j}$;

  * •

when $\mu=\uparrow\Rightarrow z_{1}-z_{2}>0\Rightarrow
2(\xi_{1}-\xi_{2})>\psi_{2}-\psi_{1}$; and when
$\mu^{\prime}=\uparrow\Rightarrow z^{\prime}_{1}-z^{\prime}_{2}>0\Rightarrow
2(\xi_{1}-\xi_{2})>\psi_{1}-\psi_{2}$. A similar condition occurs when
$\mu=\mu^{\prime}=\downarrow$. In both cases the condition
$\psi_{1}=\psi_{2}=0$ belongs to the integration domain and gives the maximum
contribution to the integral since the oscillations on the $\psi$ variables
cancel out;

  * •

when $\mu=\uparrow\Rightarrow z_{1}-z_{2}>0\Rightarrow
2(\xi_{1}-\xi_{2})>\psi_{2}-\psi_{1}$; and when
$\mu^{\prime}=\downarrow\Rightarrow z^{\prime}_{1}-z^{\prime}_{2}<0\Rightarrow
2(\xi_{1}-\xi_{2})<\psi_{1}-\psi_{2}$. A similar condition occurs when
$\mu=\downarrow$ and $\mu^{\prime}=\uparrow$. In both cases the condition
$\psi_{1}=\psi_{2}=0$ it does not belong to the integration domain and
therefore this combination of propagation modes will be subdominant in the
total contribution to the ladder term.

Because of the above two points, from the equation (36) we will only keep half
of the terms: all those that are $\uparrow\uparrow$ and $\downarrow\downarrow$
adding up all possible combinations of polarization for propagating waves in
medium II, which will give us a total of eight combinations.

Using that $\mu=\mu^{\prime}$, $s_{1}=s^{\prime}_{1}$, $s_{2}=s^{\prime}_{2}$
and that one of the integrals over $\psi$ gives a factor $2h$, the ladder term
reduces to

$\displaystyle\mathcal{L}_{ab}=2K_{I}\,k^{8}_{0}\,|F^{\uparrow}_{a\,PI}(z)|^{2}$
$\displaystyle\int\,\frac{d^{2}q}{(2\pi)^{2}}\,\int^{h}_{-h}\,d\psi\,\int^{0}_{-h}\,d\xi_{1}\,\int^{0}_{-h}\,d\xi_{2}\,$
$\displaystyle\sum_{s_{1}\,s_{2}}\,\sum_{\mu,m,m^{\prime}}\,\,C_{\mathbf{p}-\mathbf{q}}(\psi)\,C_{\mathbf{q}-\mathbf{k}_{i}}(\psi)\,\Theta\left(s_{1}(\xi_{1}-\xi_{2})\right)$
$\displaystyle\beta_{m}\,\beta^{*}_{m^{\prime}}\,\,|f_{a}(PII)|^{2}\,|f_{m}(QII)|^{2}\,|f_{m^{\prime}}(QII)|^{2}\,|f_{b}(KII)|^{2}$
$\displaystyle|X^{D\,s_{1}}_{a\,PII}|^{2}\,X^{\mu\,-s_{1}}_{m\,QII}\,X^{\bar{\mu}\,-s_{2}}_{m\,QII}\,\,X^{*\mu\,-s_{1}}_{m^{\prime}\,QII}\,X^{*\bar{\mu}\,-s_{2}}_{m^{\prime}\,QII}\,|X^{D\,s_{2}}_{b\,KII}|^{2}$
$\displaystyle e^{\imath s_{1}(PII-QII)\psi}\,e^{\imath s_{2}(QII-\,KII)\psi}$
(37)

The integrals over $d\xi_{1},d\xi_{2}$ give half the area $[-h,0]\times[-h,0]$
due to $\Theta\left(s_{1}(\xi_{1}-\xi_{2})\right)$. The integral over $\psi$
is the product between an oscillating factor and a simple function (Gaussian,
exponential decaying or step), and it is easily done analytically. From this
integration comes the term $C_{I}$,

$C_{I}(P)=\int^{h}_{-h}d\psi\,e^{\imath\psi P}\,C_{\mathbf{p}}(\psi)$ (38)

Thus, we arrive at a compact expression for the ladder term:

$\displaystyle\mathcal{L}_{ab}=2K_{I}\,k^{8}_{0}\,|F^{\uparrow}_{a\,PI}(z)|^{2}\,h^{3}\,|f_{a}(PII)|^{2}\,|f_{b}(KII)|^{2}$
$\displaystyle\int\,\frac{d^{2}q}{(2\pi)^{2}}\,\hat{C}_{\perp\,\mathbf{p}-\mathbf{q}}\,\hat{C}_{\perp\,\mathbf{q}-\mathbf{k}_{i}}$
$\displaystyle\sum_{\mu,m\,m^{\prime}}\,\beta_{m}\,\beta^{*}_{m^{\prime}}\,|f_{m}(QII)|^{2}\,|f_{m^{\prime}}(QII)|^{2}$
$\displaystyle\sum_{s_{1}\,s_{2}}\,|X^{D\,s_{1}}_{a\,PII}|^{2}\,X^{\mu\,-s_{1}}_{m\,QII}\,X^{\bar{\mu}\,-s_{2}}_{m\,QII}\,\,X^{*\mu\,-s_{1}}_{m^{\prime}\,QII}\,X^{*\bar{\mu}\,-s_{2}}_{m^{\prime}\,QII}\,|X^{D\,s_{2}}_{b\,KII}|^{2}$
$\displaystyle C_{I}(s_{1}(PII-QII)+s_{2}(QII-KII))$ (39)

The integral in the momentum variable $\mathbf{q}$ is computed numerically.
Here we see that the ladder term is real: when $m=m^{\prime}$ all the
coefficients $X$ appear as square module; when $m\neq m^{\prime}$ we see that
adding the complementary cases gives twice the real part of those terms. So we
can reduce the expression by eliminating the sum over $m^{\prime}$ and in the
case where $m^{\prime}\neq m$ we use $m^{\prime}=\bar{m}$ 111in the HH channel
it is easy: if $m=H\Rightarrow m^{\prime}=V$ and vice versa; for the VH we
have to be more careful: if $m=H\Rightarrow m^{\prime}=V$ or
$m^{\prime}=\tilde{V}1$, making it more practical to use the expression (39);
the same is true for the HV and VV channels. . Then,

$\displaystyle\mathcal{L}_{ab}=2K_{I}\,k^{8}_{0}\,|F^{\uparrow}_{a\,PI}(z)|^{2}\,h^{3}\,|f_{a}(PII)|^{2}\,|f_{b}(KII)|^{2}$
$\displaystyle\int\,\frac{d^{2}q}{(2\pi)^{2}}\,\hat{C}_{\perp\,\mathbf{p}-\mathbf{q}}\,\hat{C}_{\perp\,\mathbf{q}-\mathbf{k}_{i}}$
$\displaystyle\sum_{\mu,m}\,\sum_{s_{1}\,s_{2}}\,C_{I}(s_{1}(PII-
QII)+s_{2}(QII-KII))$
$\displaystyle|X^{\downarrow\,s_{1}}_{a\,PII}|^{2}\,\left[|f_{m}(QII)|^{4}\,|\beta_{m}|^{2}\,|X^{\mu\,-s_{1}}_{m\,QII}|^{2}\,|X^{\bar{\mu}\,-s_{2}}_{m\,QII}|^{2}\right.$
$\displaystyle\left.+2\Re\left\\{\beta_{m}\,\beta^{*}_{\bar{m}}\,|f_{m}(QII)|^{2}\,|f_{\bar{m}}(QII)|^{2}\right.\right.$
$\displaystyle\left.\left.X^{\mu\,-s_{1}}_{m\,QII}\,X^{*\mu\,-s_{1}}_{\bar{m}\,QII}\,X^{\bar{\mu}\,-s_{2}}_{m\,QII}\,X^{*\bar{\mu}\,-s_{2}}_{\bar{m}\,QII}\right\\}\right]\,|X^{\downarrow\,s_{2}}_{b\,KII}|^{2}$
(40)

The cross term is analyzed in a similar way. The details to note are that: one
of the auxiliary modes is $\mathbf{r}=\mathbf{k}+\mathbf{p}-\mathbf{q}$ and
that the vertical coordinates in the correlation functions are grouped as
$z_{1}-z^{\prime}_{2}$ and $z_{2}-z^{\prime}_{1}$. Because of this, the
propagating modes contributing to the stationary phase are
$z_{1}-z^{\prime}_{2}$. The other detail to note is that oscillations are
suppressed only under the backscattering condition ($\mathbf{p}=-\mathbf{k}$),
which implies $\mathbf{r}=-\mathbf{q}$ and $RII=QII$. Thus, the cross term
reduces to,

$\displaystyle\mathcal{C}_{ab}=2K_{I}\,k^{8}_{0}\,|F^{\uparrow}_{a\,PI}(z)|^{2}\,h^{3}\,|f_{a}(PII)|^{2}\,|f_{b}(KII)|^{2}$
$\displaystyle\int\,\frac{d^{2}q}{(2\pi)^{2}}\,\hat{C}_{\perp\,\mathbf{p}-\mathbf{q}}\,\hat{C}_{\perp\,\mathbf{q}-\mathbf{k}_{i}}$
$\displaystyle\sum_{\mu,m}\,\sum_{s_{1}\,s_{2}}\,C_{I}((s_{1}-s_{2})(KII-
QII))$
$\displaystyle|X^{\downarrow\,s_{1}}_{a\,PII}|^{2}\,\left[|f_{m}(QII)|^{4}\,|\beta_{m}|^{2}\,|X^{\mu\,-s_{1}}_{m\,QII}|^{2}\,|X^{\bar{\mu}\,-s_{2}}_{m\,QII}|^{2}\right.$
$\displaystyle\left.+2\Re\left\\{\beta_{m}\,\beta^{*}_{\bar{m}}\,|f_{m}(QII)|^{2}\,|f_{\bar{m}}(QII)|^{2}\right.\right.$
$\displaystyle\left.\left.X^{\mu\,-s_{1}}_{m\,QII}\,X^{*\mu\,-s_{1}}_{\bar{m}\,QII}\,X^{\bar{\mu}\,-s_{2}}_{m\,QII}\,X^{*\bar{\mu}\,-s_{2}}_{\bar{m}\,QII}\right\\}\right]\,|X^{\downarrow\,s_{2}}_{b\,KII}|^{2}$
(41)

By comparing (40) with (41), in backscattering, i.e. $PII=KII$, the cross and
ladder terms are equal, leading to the so-called backscattering enhancement
[14].

With the expressions (40) and (41) the ladder and cross terms of the four
channels can be computed. We only have to take into account that, depending on
the channel, the index $m$ can take one or three more values (VH adds
$\tilde{V}1$, HV $\tilde{V}2$ and VV
$\tilde{V}_{1},\tilde{V}_{2},\tilde{V}_{12}$)); see 34)).

###  6.3 Bistatic scattering coefficients

So far we have presented the scattering amplitude at first and second order
for both incident and scattered polarization. These amplitudes have been
called $S_{HH}$ and $S_{VH}$ when we have a TE incident mode and the scattered
field has $H$ ($\varphi$) or $V$ ($\psi$) polarization; the other two
combinations, $S_{HV}$ and $S_{VV}$, occur with a TM incident mode. The
Poynting vector is $\mathbf{S}=\hat{\mathbf{k}}\,|\mathbf{E}|^{2}$. Now we
introduce the bistatic scattering coefficient as,

$\sigma_{ab}(\mathbf{p},\mathbf{k})=\frac{1}{\mathcal{A}}\frac{\langle\mathbf{p}\cdot\mathbf{S}_{a}(\mathbf{p})\rangle_{\text{Incoh}}}{|\mathbf{k}\cdot\mathbf{E}_{b}(\mathbf{k})|^{2}}$
(42)

This magnitude represents the incoherent scattered energy with polarization
$a$ in direction $\mathbf{p}$ normalized by unit area and incident energy
where the incident wave has polarization $b$ and propagates in direction
$\mathbf{k}$.

As we discuss in Section 66.1 at first order the scattered energy is already
incoherent. At second order we defined the incoherent scattered energy through
equation (6.2), which we then write as the sum of the ladder and cross terms.
Both first order mean energy, ladder and cross term are proportional to
$\delta_{\mathbf{k}}(\mathbf{0})$ which represents the illuminated area by the
incident wave. In addition to this, we have developed all our calculation
using a unit incident electric field, so the normalization factor in (42) only
introduces a factor $\cos(\theta_{i})$ in the denominator in the ratio between
scattered and incident energies. Also, we are interested in analyzing the
results in the backscattering condition ($\mathbf{p}=-\mathbf{k}$) where
ladder and cross terms are equal. Therefore, the bistatic scattering
coefficients will be written using expressions (6.1) and (40). Both
expressions depend, respectively, on the correlation function
$C_{\mathbf{p}-\mathbf{k}}(z_{1}-z_{2})$ and its 2D Fourier transform
$\hat{C}_{\perp\mathbf{p}}\,C_{I}(P)$. The first of these quantities
represents the fluctuations of dielectric permittivity in region II and was
defined in (2.2) and the second one was defined in (38).

To get a closed form to (6.1) and (40) we choose an explicit expression for
the correlation function $C_{\mathbf{p}}(z)$. To illustrate our results we
will use three kinds of correlations:

$C(\mathbf{r}-\mathbf{r}^{\prime})=s\begin{cases}\exp\\{-\frac{|\mathbf{x}-\mathbf{x}^{\prime}|^{2}}{l^{2}_{r}}\\}\,\exp\\{-\frac{(z-z^{\prime})^{2}}{l^{2}_{z}}\\}\\\
\exp\\{-\frac{|\mathbf{x}-\mathbf{x}^{\prime}|^{2}}{l^{2}_{r}}\\}\,\exp\\{-\frac{|z-z^{\prime}|}{l_{z}}\\}\\\
\Theta\left[1-l^{-2}_{r}|\mathbf{x}-\mathbf{x}^{\prime}|^{2}-l^{-2}_{z}(z-z^{\prime})^{2}\right]\end{cases}$
(43)

which correspond to Gaussian, Exponential and ellipsoidal Heaviside-like,
respectively. In the above $s=\delta\,\Re\\{\epsilon_{0}\\}^{2}$ represents
the fluctuation intensity, $l_{r}$ ($l_{z}$) is the radial (vertical)
correlation length. These three functions have the following 2D Fourier
transforms:

$C_{\mathbf{p}}(z-z^{\prime})=s\begin{cases}\frac{l^{2}_{r}}{(2\pi)^{2}}\,\exp\\{-\frac{l^{2}_{r}}{4}\mathbf{p}^{2}\\}\,\exp\\{-\frac{(z-z^{\prime})^{2}}{l^{2}_{z}}\\}\\\
\frac{l^{2}_{r}}{(2\pi)^{2}}\,\exp\\{-\frac{l^{2}_{r}}{4}\mathbf{p}^{2}\\}\,\exp\\{-\frac{|z-z^{\prime}|}{l_{z}}\\}\\\
\frac{l_{r}}{2\pi\,|\mathbf{p}|}\,g(z,z^{\prime})\,J_{1}\left(l_{r}|\mathbf{p}|g(z,z^{\prime})\right)\end{cases}$
(44)

where in the last line $g(z,z^{\prime})=\sqrt{1-(z-z^{\prime})^{2}/l^{2}_{z}}$
and $J_{1}(x)$ is the Bessel function of the first kind. From this expression
it is straightforward to find $C_{I}(P)$.

Also we consider a correlation function for modeling overlapping randomly
distributed cylinders which have radius $R$ and height $H$ and low space
density $\rho$. To define this correlation we assume a complete vertical
overlap through the entire thickness of region II and that there is no
overlapping if the distance between center of the two cylinders are bigger
than $2R$; when this distance is lower than $2R$ we write the correlation as
the intersecting area between the circular sections of the two cylinders.
Then, this correlation can be written as follows,

$C(\mathbf{r}-\mathbf{r}^{\prime})=\rho\,Z\,\Theta\left(Z\right)\,P(\chi)\,\Theta(1-\chi)$
(45)

where $\chi=\frac{|r-r^{\prime}|}{2R}$,
$P(\chi)=\arccos(\chi)-\chi\sqrt{1-\chi^{2}}$ and
$Z=1-\frac{|z-z^{\prime}|}{2H}$. Its 2D Fourier transform is,

$C_{\mathbf{p}}(z-z^{\prime})=\frac{2}{\pi}\,\rho\,Z\,\Theta(Z)\left[\frac{J_{1}(R|\mathbf{p}|)}{|\mathbf{p}|}\right]^{2}$
(46)

At this point it is important to mention that what has been developed so far
can easily be readapted to calculate both the coherence matrix $T$ or the
covariance matrix $C$. Both matrices involve the calculation of mean values
between the different $S_{qp}$ amplitudes. For example, to calculate $\langle
S_{HH}\,S^{*}_{VV}\rangle$ we can use the results given by equations (6.1) and
(6.2) by simply replacing one of the amplitudes by Sab by Srs and obtain the
desired mean value. This change will only modify geometrical coefficients or
reflection or transmission amplitudes, but will not affect in any way the
development of the calculations made in the two previous subsections.

In the next section we will show the results obtained from the set up we have
presented so far.

##  7 Results

We illustrate our results through numerical simulations of different
backscattering scenarios, where we either fix the dielectric and geometric
parameters and vary the incidence angle, or else we consider layered media
with different geometric o dielectric parameters at a given incidence angle.

First we show in Figure 3 the backscattering coefficient as a function of
incident angle. We use the following parameters: incident wavelength
$\lambda=0.24\leavevmode\nobreak\ \text{m}$; homogeneous layered media
composed by relative dielectric permittivity $\epsilon_{0I}=1$,
$\epsilon_{0II}=1.06+\imath\,0.016$ (which corresponds to a low density
vegetated medium), $\epsilon_{0III}=6+\imath\,0.6$ and height
$h=1.05\leavevmode\nobreak\ \text{m}$ (describing a dry subsurface soil);
correlation functions with $l_{r}=8\leavevmode\nobreak\ \text{mm}$,
$l_{z}=8.5\leavevmode\nobreak\ \text{cm}$, $s=2.5\times
10^{-3}\,(\Re\\{\epsilon_{0II}\\})^{2}$. In the case of using two overlapping
cylinders we consider their height equal to the thickness of region II and
also that their the density is $\rho=2.5\times 10{-3}$. These results were
computed using a desktop computer which has a processor Intel i7-9700 with
eight kernels. The complete Figure 3 took 31.84 seconds. Each of the four
curves in the four panels has 45 points, therefore, each point requires 0.043
seconds of calculation. These times make the results presented here plausible
for use in a Bayesian inference scheme. the computation times for the figures
presented below were similar.

Figure 3: Normalized backscattering coefficient as a function of incident
angle. Solid, dash-dot and dashed line correspond to Gaussian, Exponential or
Heaviside-like correlation function, respectively. Solid-dot line correspond
to the correlation function describing overlapping cylinders.

Here we see that, although the trend of the four channels is the same with
respect to the angle of incidence, the level of the backscattered signal
increases considerably when using the correlation given by the cylinder
overlap (see equation (46)).

Then, in Figure 4 we show the backscattering coefficient as a function of
layer thickness of region II for three different incident angles. We use the
same parameters as in Figure 3 as corresponds.

Figure 4: Normalized backscattering coefficient as a function of layer
thickness for three different incident angles. Line style indicates the same
correlation functions as in Figure 3.

Later, in Figure 5 we show the backscattering coefficient as a function of
incident frequency for three different values for the thickness of the region
II. We use the same values for the dielectric permittivity as in Figure 3,
correlation functions with $l_{r}=8\leavevmode\nobreak\ \text{cm}$,
$l_{z}=5\leavevmode\nobreak\ \text{mm}$, $s=2.5\times
10^{-3}\,(\Re\\{\epsilon_{0II}\\})^{2}$.

Figure 5: Normalized backscattering coefficient as a function of incident
frequency for three different thickness of region II. Line style indicates the
same correlation functions as in Figure 3.

##  8 Conclusions

In this paper we present a fully coherent, analytic model of the
backscattering intensity in all HH, HV, VH and VV channels, for the volume
scattering of radiation from a layer of finite thickness, such as a vegetation
layer over bare soil. We aim for a simple, not numerically intensive model
which could be used either as forward model in a Bayesian estimation scheme,
or else as a preliminary means to identify key features of a concrete problem,
for its further analysis by more sophisticated theoretical and numerical
approaches. Coherence is required to make full use of available observational
capabilities, while simplicity is desirable given the large number of
iterations of the model required by Bayesian estimation.

We have checked our model by computing the scattering intensities for four
simple functional forms of the dielectric constant fluctuations correlation
function, 43 and 45. The results agree well with expectations from actual
observations and show overall a common pattern. However there is a clear
intensity gap between the results for the correlation function of overlapping
cylinders 45 and the other ones 43, and it would be an easy exercise to
imagine further analytical expressions for the correlation function giving
intermediate results between these extremes. This shows that a proper choice
of the two-point correlation of the dielectric constant fluctuation is the
single most important factor which determines the accuracy of the model, with
an increased complexity in modeling (such as computing further orders in
perturbation theory, including non Gaussian fluctuations, using a numerically
intensive approach or replacing the analytical model by a numerical solution
of Maxwell’s equations from scratch) bringing only incremental gains.

Moreover, while in this paper we have focused on the backscattered
intensities, the model yields full information about the scattered field phase
and thus it can also be used in polarimetric analysis, which we will discuss
in a separate contribution.

##  9 Founding Information

Work supported in part by CONICET, Universidad de Buenos Aires UBACYT
20020220300204BA, CONICET PIP 11220210100595CO, ANPCYT PICT 2018-03684 and
CONICET PICT 2020-01830.

## References

  * [1] P. M. Lee, _Bayesian statistics_ (Oxford University Press London:, 1989).
  * [2] P. Gregory, _Bayesian logical data analysis for the physical sciences: a comparative approach with mathematica® support_ (Cambridge University Press, 2005).
  * [3] N. Heard _et al._ , _An introduction to Bayesian inference, methods and computation_ (Springer, 2021).
  * [4] L. Tsang, J. A. Kong, and K.-H. Ding, _Scattering of electromagnetic waves: theories and applications_ , vol. 15 (John Wiley & Sons, 2000).
  * [5] L. Tsang and J. A. Kong, _Scattering of electromagnetic waves: advanced topics_ (John Wiley & Sons, 2004).
  * [6] L. Tsang, J. A. Kong, K.-H. Ding, and C. O. Ao, _Scattering of electromagnetic waves: numerical simulations_ (John Wiley & Sons, 2004).
  * [7] D. Long and F. Ulaby, _Microwave radar and radiometric remote sensing_ (Artech, 2015).
  * [8] A. G. Voronovich, _Wave scattering from rough surfaces_ , vol. 17 (Springer Science & Business Media, 2013).
  * [9] N. Pinel and C. Boulier, _Electromagnetic wave scattering from random rough surfaces: Asymptotic models_ (John Wiley & Sons, 2013).
  * [10] A. Ishimaru, _Electromagnetic wave propagation, radiation, and scattering: from fundamentals to applications_ (John Wiley & Sons, 2017).
  * [11] F. T. Ulaby, R. K. Moore, and A. K. Fung, _Microwave remote sensing: Active and passive. volume 1-microwave remote sensing fundamentals and radiometry_ (Addison-Wesley, 1981).
  * [12] R. Shin and J. A. Kong, “Radiative transfer theory for active remote sensing of two-layer random medium,” Progress In Electromagnetics Research 1, 359–417 (1989).
  * [13] M. Zuniga, J. Kong, and L. Tsang, “Depolarization effects in the active remote sensing of random media,” Journal of Applied Physics 51, 2315–2325 (1980).
  * [14] Y. Jin and J. Kong, “Ladder and cross terms in second-order distorted born approximation,” Journal of mathematical physics 26, 994–1011 (1985).
  * [15] M. Borgeaud, R. T. Shin, and J. A. Kong, “Theoretical models for polarimetric radar clutter,” Journal of Electromagnetic Waves and Applications 1, 73–89 (1987).
  * [16] M. Borgeaud, S. Nghiem, R. Shin, and J. Kong, “Theoretical models for polarimetric microwave remote sensing of earth terrain,” Journal of Electromagnetic Waves and Applications 3, 61–81 (1989).
  * [17] L. Morena, K. James, and J. Beck, “An introduction to the radarsat-2 mission,” Canadian Journal of Remote Sensing 30, 221–234 (2004).
  * [18] M. Azcueta, J. P. C. Gonzalez, T. Zajc, J. Ferreyra, and M. Thibeault, “External calibration results of the saocom-1a commissioning phase,” IEEE Transactions on Geoscience and Remote Sensing 60, 1–8 (2022).
  * [19] J.-S. Lee and E. Pottier, _Polarimetric radar imaging: from basics to applications_ (CRC press, 2017).
  * [20] J. Arellana, M. Franco, and F. Grings, “Using saocom data and bayesian inference to estimate soil dielectric constant in agricultural soils,” IEEE Geoscience and Remote Sensing Letters 20, 1–5 (2023).
  * [21] A. Tabatabaeenejad and M. Moghaddam, “Bistatic scattering from three-dimensional layered rough surfaces,” IEEE Transactions on Geoscience and Remote Sensing 44, 2102–2114 (2006).
  * [22] M. Sanamzadeh, L. Tsang, J. T. Johnson, R. J. Burkholder, and S. Tan, “Scattering of electromagnetic waves from 3d multilayer random rough surfaces based on the second-order small perturbation method: energy conservation, reflectivity, and emissivity,” JOSA A 34, 395–409 (2017).
  * [23] M. Franco and E. Calzetta, “Wave propagation in non-gaussian random media,” Journal of Physics A: Mathematical and Theoretical 48, 045206 (2015).
  * [24] F. Lamagna and E. Calzetta, “A functional renormalization method for wave propagation in random media,” Journal of Physics A: Mathematical and Theoretical 50, 315102 (2017).
  * [25] B. Van Tiggelen, D. Wiersma, and A. Lagendijk, “Self-consistent theory for the enhancement factor in coherent backscattering,” Europhysics Letters 30, 1 (1995).
  * [26] A. Knothe and T. Wellens, “Flux conservation in coherent backscattering and weak localization of light,” Journal of Physics A: Mathematical and Theoretical 46, 315101 (2013).
  * [27] S.-B. Kim, M. Moghaddam, L. Tsang, M. Burgin, X. Xu, and E. G. Njoku, “Models of l-band radar backscattering coefficients over global terrain for soil moisture retrieval,” IEEE Transactions on Geoscience and Remote Sensing 52, 1381–1396 (2013).
  * [28] H. Huang, T.-H. Liao, S.-B. Kim, X. Xu, L. Tsang, T. J. Jackson, and S. H. Yueh, “L-band radar scattering and soil moisture retrieval of wheat, canola and pasture fields for smap active algorithms,” Progress In Electromagnetics Research 170, 129–152 (2021).
  * [29] L. Tsang, T.-H. Liao, R. Gao, H. Xu, W. Gu, and J. Zhu, “Theory of microwave remote sensing of vegetation effects, soop and rough soil surface backscattering,” Remote Sensing 14, 3640 (2022).
  * [30] J. Jeong, L. Tsang, X. Xu, S. Yueh, and S. A. Margulis, “Full-wave electromagnetic model simulations of p-band radio signal propagation through forest using the fast hybrid method,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing (2023).
  * [31] J. Jeong, L. Tsang, W. Gu, A. Colliander, and S. H. Yueh, “Wave propagation in vegetation field by combining fast multiple scattering theory and numerical electromagnetics in a hybrid method,” IEEE Transactions on Antennas and Propagation 71, 3598–3610 (2023).
  * [32] B. A. van Tiggelen and R. Maynard, “Reciprocity and coherent backscattering of light,” in _Wave Propagation in Complex Media,_ (Springer, 1998), pp. 247–271.