# Wigner measures of electromagnetic waves in bianisotropic heterogeneous
media

Jean-Luc Akian  ONERA–The French Aerospace Lab, France  jean-
<EMAIL_ADDRESS>and  Éric Savin  ONERA–The French Aerospace Lab, France
<EMAIL_ADDRESS>
###### Abstract.

We study the propagation of high-frequency electromagnetic waves in randomly
heterogeneous bianisotropic media with dissipative properties. For that
purpose we consider randomly fluctuating optical responses of such media with
correlation lengths comparable to the typical wavelength of the waves.
Although the fluctuations are weak, they induce multiple scattering over long
propagation times and/or distances such that the waves end up travelling in
many different directions with mixed polarizations. We derive the dispersion
and evolution properties of the Wigner measure of the electromagnetic fields,
which describes their angularly-resolved energy density in this propagation
regime. The analysis starts from Maxwell’s equations with general constitutive
equations. We first ignore the random fluctuations of the optical response and
obtain uncoupled transport equations for the components of the Wigner measure
on the different propagation modes (polarizations). Then we use a multi-scale
expansion of the Wigner mesure to obtain the radiative transfer equations
satisfied by these components when the fluctuations are no longer ignored. The
radiative transfer equations are coupled through their collisional parts,
which account for the scattering of waves by the random fluctuations and their
possible changes in polarization. The collisional kernels describing these
processes depend on the power and cross-power spectral densities of the
fluctuations at the wavelength scale. The overall derivation is based on the
interpretation of Wigner transforms and Wigner measures in terms of
semiclassical pseudo-differential operators in their standard quantization.

###### Key words and phrases:

Maxwell’s equations, Bianisotropic dielectric media, Kinetic model, Transport
equation, Radiative transfer

Corresponding author: É. Savin, ONERA–The French Aerospace Lab, 6 chemin de la
Vauve aux Granges, FR-91123 Palaiseau cedex, France (Eric.Savin@onera.fr).

##  1\. Introduction

The analysis of electromagnetic wave propagation in complex, heterogeneous
media has relevance to a wide range of engineering applications, including
radar imaging and remote sensing [48], optical imaging [13], atmospheric
physics [30], or communications [25], among others. As the complexity of the
propagation media can not be known exactly but only assessed in a statistical
sense, we model their inhomogeneities (or imperfections) as random fields and
thus consider Maxwell’s equations with randomly fluctuating electromagnetic
parameters. This approach pertains to waveguides filled with dielectric
materials having numerous imperfections, atmospheres and clouds, colloidal
suspensions, biological tissues, or plasmas, which is the type of media we
have more particularly in mind here.

Asymptotic methods can be used to characterize both quantitatively and
qualitatively the wave fields when some explicit scales can be identified and
separated to describe the corresponding propagation regimes [3, 45, 47].
Recent mathematical works on electromagnetic waves in random media have
considered layered half-spaces [22], waveguides [4], or beam propagation in
open environments [11], for example. In these works a preferred direction of
propagation can be identified, and the electromagnetic field can be decomposed
into plane wave components transverse to that direction. Their amplitudes are
random fields of which evolution is driven by the random fluctuations of the
dielectric parameters and can be studied in different asymptotic regimes. For
example, the cumulative scattering effects induced by weak inhomogeneities at
large distances of propagation in a waveguide or a beam are typically
manifested in the exponential decay of the average fields, the enhancement of
their random fluctuations, and depolarization.

We rather investigate in this paper the situation where no preferred direction
of propagation can be identified as waves travel in open random media. The
wavelength and the correlation length of the random fluctuations of the
electromagnetic parameters are small and comparable, while the propagation
distance and/or time lapse are large. The amplitude of the random fluctuations
is also weak. This high-frequency scaling is of interest if the objective is
for instance to probe some inhomogeneities from remote sensors. In this
setting, waves are multiply scattered by the inhomogeneities as they travel in
all directions, and their polarization and phase get randomized. It is then
known that the energy density is the relevant quantity to focus on, and that
it can be evaluated from the Wigner transform of the wave field and its non-
negative high-frequency limit measure–the Wigner measure [23, 35, 37, 53]. The
evolution of the Wigner measure is described by radiative transfer equations
[15, 30, 38, 51] with linear collisional operators of which kernels depend on
the correlation structure of the inhomogeneities [9, 46]. Actually, the sole
information needed to account for the fluctuations of the electromagnetic
parameters is their power spectral densities, _i.e._ the spatial Fourier
transforms of their auto-correlation functions. These radiative transfer
equations for electromagnetic vector waves were already derived in [39, 46]
for random fluctuations in isotropic media, and subsequently in [18, 19] for
random fluctuations in bianisotropic media with homogeneous mean background
parameters and no dissipation. Our aim is to extend these results to more
general constitutive equations, for heterogeneous materials with dispersive
and dissipative properties.

Radiative transfer equations have a geometrical interpretation [34] notably in
terms of bicharacteristic rays that is, rays in phase space, which parallels
classical asymptotic ray methods [7, 32]. The latter are applicable to
synthetic aperture radar (SAR) simulation for example [6], including the
investigation of multipath and vibration signatures in radar imaging and
remote sensing or the computation of antenna transfer functions in
heterogeneous media whenever the approximations of geometrical optics hold.
Ray tracing solvers are needed in this respect [24]. One of the possible
applications of this research concerns the development of ray tracing
algorithms in phase space for complex media such as plasmas [50]. Ray tracing
has been intensively developed in the last decade in the context of computer
graphics [41, 28], notably, and one also expects that the application
programming interfaces that have been implemented in the video game industry,
among others, could be used in the context of computational electromagnetics
alike.

The rest of the paper is organized as follows. In Sect. 2 we introduce the
basic physical framework and notations used throughout, and we summarize our
main results. The relevance of considering a Wigner transform and its non-
negative limit measure for the analysis of high-frequency wave propagation
phenomena is emphasized. We first consider the transport regime in the absence
of fluctuations of the optical response that models instantaneous effects, and
then the radiative transfer regime when weak random fluctuations are taken
into account. The detailed derivation of the transport equations is outlined
in Sect. 3, and the detailed derivation of the radiative transfer equations is
outlined in Sect. 4. We use pseudo-differential calculus and the
interpretation of Wigner transforms in terms of semiclassical pseudo-
differential operators, choosing the same quantization for both of them. We
argue that this choice clarifies the analyses as compared to other approaches
where the Wigner transforms and the operators have different quantizations. A
summary and outlook are finally drawn in Sect. 5.

##  2\. Maxwell’s equations in bianisotropic media and main results

###  2.1. Physical setting

We wish to analyze high-frequency electromagnetic wave propagation phenomena
in general heterogeneous, bianisotropic dielectric media. The Ampère and
Faraday laws read:

(1) $\frac{\partial{\boldsymbol{D}}}{\partial
t}={\boldsymbol{\nabla}}\times{\boldsymbol{H}}-{\boldsymbol{J}}\,,\quad\frac{\partial{\boldsymbol{B}}}{\partial
t}=-{\boldsymbol{\nabla}}\times{\boldsymbol{E}}\,,$

where ${\boldsymbol{B}}$ is the magnetic flux density, ${\boldsymbol{D}}$ is
the electric flux density, or electric displacement field, ${\boldsymbol{E}}$
is the electric field, ${\boldsymbol{H}}$ is the magnetic field, and
${\boldsymbol{J}}$ is the density of electric current. Taking the divergence
of Eq. (1) the fluxes are subjected to Gauss laws:

(2)
${\boldsymbol{\nabla}}\cdot{\boldsymbol{D}}=\rho\,,\quad{\boldsymbol{\nabla}}\cdot{\boldsymbol{B}}=0\,,$

for $\rho$ being the density of electric charge, which is then related to the
current by the equation of continuity of charge:

$\frac{\partial\rho}{\partial
t}+{\boldsymbol{\nabla}}\cdot{\boldsymbol{J}}=0\,.$

Here we consider open media in $\mathbb{R}^{3}$ with no external charge
($\rho=0$) and no external current (${\boldsymbol{J}}={\bf 0}$). Besides,
Gauss laws hold for all times provided that they hold at some initial time.
The Ampère and Faraday laws are supplemented with constitutive equations for
general linear dispersive media with the properties of causality, time-
invariance, and locality (see _e.g._ [14, 29, 31]), which read:

(3)
$\begin{split}{\boldsymbol{D}}({\boldsymbol{x}},t)&={\boldsymbol{\epsilon}}_{0}({\boldsymbol{x}}){\boldsymbol{E}}({\boldsymbol{x}},t)+({\boldsymbol{\epsilon}}_{d}\star{\boldsymbol{E}})({\boldsymbol{x}},t)+{\boldsymbol{\xi}}_{0}({\boldsymbol{x}}){\boldsymbol{H}}({\boldsymbol{x}},t)+({\boldsymbol{\xi}}_{d}\star{\boldsymbol{H}})({\boldsymbol{x}},t)\,,\\\
{\boldsymbol{B}}({\boldsymbol{x}},t)&={\boldsymbol{\zeta}}_{0}({\boldsymbol{x}}){\boldsymbol{E}}({\boldsymbol{x}},t)+({\boldsymbol{\zeta}}_{d}\star{\boldsymbol{E}})({\boldsymbol{x}},t)+{\boldsymbol{\mu}}_{0}({\boldsymbol{x}}){\boldsymbol{H}}({\boldsymbol{x}},t)+({\boldsymbol{\mu}}_{d}\star{\boldsymbol{H}})({\boldsymbol{x}},t)\,,\end{split}$

where ${\boldsymbol{\epsilon}}_{0}({\boldsymbol{x}})$ is the spatially
variable permittivity tensor, ${\boldsymbol{\mu}}_{0}({\boldsymbol{x}})$ is
the spatially variable permeability tensor, and
${\boldsymbol{\xi}}_{0}({\boldsymbol{x}})$ and
${\boldsymbol{\zeta}}_{0}({\boldsymbol{x}})$ are spatially variable
magnetoelectric tensors. Also $\star$ stands for the time convolution product:

$({\boldsymbol{\epsilon}}_{d}\star{\boldsymbol{E}})({\boldsymbol{x}},t)=\int_{0}^{t}{\boldsymbol{\epsilon}}_{d}({\boldsymbol{x}},\tau){\boldsymbol{E}}({\boldsymbol{x}},t-\tau){\mathrm{d}}\tau$

with similar expressions for the other terms. In Eq. (3)
${\boldsymbol{P}}=\smash{{\boldsymbol{\epsilon}}_{d}\star{\boldsymbol{E}}}$
stands for the polarization field and
${\boldsymbol{M}}=\smash{{\boldsymbol{\mu}}_{0}^{-1}({\boldsymbol{\mu}}_{d}\star{\boldsymbol{H}})}$
stands for the magnetization field, but these constitutive equations also
account for coupled magnetoelectric effects through the magnetoelectric
tensors. Here we consider the subclass of bianisotropic materials for which
the following symmetry relationships hold:
$\smash{{\boldsymbol{\zeta}}_{0}}=\smash{{{\boldsymbol{\xi}}}^{*}_{0}}$,
_i.e._ the conjugate transpose of $\smash{{\boldsymbol{\xi}}_{0}}$; and
$\smash{{\boldsymbol{\epsilon}}_{0}}=\smash{{{\boldsymbol{\epsilon}}}^{*}_{0}}$,
$\smash{{\boldsymbol{\mu}}_{0}}=\smash{{{\boldsymbol{\mu}}}^{*}_{0}}$. These
symmetry conditions reduce the number of independent constitutive parameters
in ${\boldsymbol{\epsilon}}_{0}$, $\smash{{\boldsymbol{\mu}}_{0}}$, and
$\smash{{\boldsymbol{\xi}}_{0}}$ to 21.

Now plugging Eq. (3) into Eq. (1), we arrive at Maxwell’s $6\times 6$ system
for the electric field ${\boldsymbol{E}}$ and magnetic field
${\boldsymbol{H}}$:

(4) $\frac{\partial}{\partial
t}\left(\begin{bmatrix}{\boldsymbol{\epsilon}}_{0}&{\boldsymbol{\xi}}_{0}\\\
{{\boldsymbol{\xi}}}^{*}_{0}&{\boldsymbol{\mu}}_{0}\end{bmatrix}+\begin{bmatrix}{\boldsymbol{\epsilon}}_{d}&{\boldsymbol{\xi}}_{d}\\\
{\boldsymbol{\zeta}}_{d}&{\boldsymbol{\mu}}_{d}\end{bmatrix}\star\right)\begin{pmatrix}{\boldsymbol{E}}\\\
{\boldsymbol{H}}\end{pmatrix}+{\boldsymbol{\nabla}}\times\begin{bmatrix}{\bf
0}&-{\boldsymbol{I}}\\\ {\boldsymbol{I}}&{\bf
0}\end{bmatrix}\begin{pmatrix}{\boldsymbol{E}}\\\
{\boldsymbol{H}}\end{pmatrix}={\bf 0}\,,$

where ${\boldsymbol{I}}$ is the $3\times 3$ identity matrix. We introduce the
following matrices of electromagnetic tensors:

(5)
${\boldsymbol{K}}_{0}({\boldsymbol{x}})=\begin{bmatrix}{\boldsymbol{\epsilon}}_{0}({\boldsymbol{x}})&{\boldsymbol{\xi}}_{0}({\boldsymbol{x}})\\\
{{\boldsymbol{\xi}}}^{*}_{0}({\boldsymbol{x}})&{\boldsymbol{\mu}}_{0}({\boldsymbol{x}})\end{bmatrix}\,,\quad{\boldsymbol{K}}_{d}({\boldsymbol{x}},t)=\begin{bmatrix}{\boldsymbol{\epsilon}}_{d}({\boldsymbol{x}},t)&{\boldsymbol{\xi}}_{d}({\boldsymbol{x}},t)\\\
{\boldsymbol{\zeta}}_{d}({\boldsymbol{x}},t)&{\boldsymbol{\mu}}_{d}({\boldsymbol{x}},t)\end{bmatrix}\,.$

$\smash{{\boldsymbol{K}}_{0}}$ is called the optical response, or
(generalized) susceptibility tensor, and models instantaneous effects, while
$\smash{{\boldsymbol{K}}_{d}}$ is called the (generalized) susceptibility
kernel, and models memory as well as dissipation effects [29, 31]. We assume
throughout the paper that $\smash{{\boldsymbol{K}}_{0}}$, which is Hermitian,
is also positive definite, as in [18, 19]. A negative definite optical
response gives rise to a negative index of refraction, as in metamaterials
[52].

######  Example 1.

An isotropic, lossless medium is
$\smash{{\boldsymbol{K}}_{0}({\boldsymbol{x}})}=\smash{\operatorname{diag}(\epsilon_{0}({\boldsymbol{x}}){\boldsymbol{I}},\mu_{0}({\boldsymbol{x}}){\boldsymbol{I}})}$
and $\smash{{\boldsymbol{K}}_{d}}={\bf 0}$. A chiral medium is a reciprocal
bi-isotropic medium with
${\boldsymbol{\epsilon}}_{0}=\epsilon_{0}{\boldsymbol{I}}$,
${\boldsymbol{\mu}}_{0}=\mu_{0}{\boldsymbol{I}}$, and a purely imaginary
magnetoelectric tensor ${\boldsymbol{\xi}}_{0}=\mathrm{i}\xi{\boldsymbol{I}}$,
where $\smash{\epsilon_{0}}$ and $\smash{\mu_{0}}$ are the permittivity and
permeability constants, and $\xi\in\mathbb{R}$ is the magnetoelectric
constant. An example of a dissipative, local homogeneous material (in the
terminology of [14]) is the Lorentz model with damping in [44] with:

(6) ${\boldsymbol{K}}_{0}=\begin{bmatrix}\epsilon_{0}{\boldsymbol{I}}&{\bf
0}\\\ {\bf
0}&\mu_{0}{\boldsymbol{I}}\end{bmatrix}\,,\quad\widehat{{\boldsymbol{K}}}_{d}(\omega)=\begin{bmatrix}\widehat{\epsilon}_{d}(\omega){\boldsymbol{I}}&{\bf
0}\\\ {\bf 0}&{\bf 0}\end{bmatrix}\,,$

where
$\widehat{\epsilon}_{d}(\omega)=\smash{\frac{\epsilon_{0}\omega_{p}^{2}}{-\omega^{2}+\mathrm{i}\omega\Gamma+\omega_{0}^{2}}}$
is the susceptibility (in frequency domain), $\omega_{p}$ is the plasma
frequency such that
$\smash{\omega_{p}^{2}}=\smash{\frac{Ne^{2}}{m\epsilon_{0}}}$, $\omega_{0}$ is
a characteristic frequency for the motion of an electron with mass $m$ and
charge $e$, $N$ is the electron density, and $\Gamma$ is a damping loss rate.
The Fourier transform in time domain is
$\smash{\widehat{{\boldsymbol{K}}}_{d}}(\omega)=\smash{\int_{\mathbb{R}}\operatorname{e}^{-\mathrm{i}\omega
t}{\boldsymbol{K}}_{d}(t){\mathrm{d}}t}$. Drude material is $\Gamma=0$ and
$\smash{\omega_{0}=0}$. Other examples are described in e.g. [14, 18, 19, 31].

We then introduce the scalar product
$\left\langle{\boldsymbol{u}},{\boldsymbol{v}}\right\rangle={({\boldsymbol{K}}_{0}{\boldsymbol{u}})}^{*}{\boldsymbol{v}}={{\boldsymbol{u}}}^{*}{\boldsymbol{K}}_{0}{\boldsymbol{v}}$
(for $\smash{{{\boldsymbol{K}}}^{*}_{0}}=\smash{{\boldsymbol{K}}_{0}}$ where
$\smash{{{\boldsymbol{K}}}^{*}_{0}}$ stands for the conjugate transpose of
$\smash{{\boldsymbol{K}}_{0}}$), and the electromagnetic energy density
${\mathcal{E}}({\boldsymbol{x}},t)$ and the energy flow density (Poynting
vector) $\boldsymbol{{\mathcal{F}}}({\boldsymbol{x}},t)$:

(7)
$\begin{split}{\mathcal{E}}({\boldsymbol{x}},t)&=\frac{1}{2}\left\langle\begin{pmatrix}{\boldsymbol{E}}\\\
{\boldsymbol{H}}\end{pmatrix},\begin{pmatrix}{\boldsymbol{E}}\\\
{\boldsymbol{H}}\end{pmatrix}\right\rangle\\\
&=\frac{1}{2}{{\boldsymbol{E}}({\boldsymbol{x}},t)}^{*}{\boldsymbol{\epsilon}}_{0}({\boldsymbol{x}}){\boldsymbol{E}}({\boldsymbol{x}},t)+\frac{1}{2}{{\boldsymbol{H}}({\boldsymbol{x}},t)}^{*}{\boldsymbol{\mu}}_{0}({\boldsymbol{x}}){\boldsymbol{H}}({\boldsymbol{x}},t)\\\
&\quad\quad+\Re\mathfrak{e}\\{{{\boldsymbol{E}}({\boldsymbol{x}},t)}^{*}{\boldsymbol{\xi}}_{0}({\boldsymbol{x}}){\boldsymbol{H}}({\boldsymbol{x}},t)\\}\,,\\\
\boldsymbol{{\mathcal{F}}}({\boldsymbol{x}},t)&=\overline{{\boldsymbol{E}}({\boldsymbol{x}},t)}\times{\boldsymbol{H}}({\boldsymbol{x}},t)\,,\end{split}$

such that in the undamped case
$\smash{{\boldsymbol{K}}_{d}({\boldsymbol{x}},t)}\equiv{\bf 0}$:

$\frac{\partial{\mathcal{E}}}{\partial
t}+{\boldsymbol{\nabla}}\cdot\boldsymbol{{\mathcal{F}}}=0\,,\quad\frac{{\mathrm{d}}}{{\mathrm{d}}t}\int{\mathcal{E}}({\boldsymbol{x}},t){\mathrm{d}}{\boldsymbol{x}}=0\,.$

The conservation law above establishes how the energy density is spread in
space, but it does not describe how it propagates. Our main objective in this
paper is to describe how this quantity evolves both in space and direction for
high-frequency waves induced by, say, strongly oscillating initial conditions,
and for a randomly fluctuating optical response and a non-vanishing
susceptibility kernel $\smash{{\boldsymbol{K}}_{d}}$. In the next section we
show how this goal can be achieved using Wigner transforms of the wave fields
and their high-frequency limits. We illustrate it with the wave equation with
constant speed, before we turn to Maxwell’s system (4) and state our main
results in Sect. 2.3.

###  2.2. High-frequency limit and Wigner transform

Maxwell’s equations (4) are subjected to highly oscillating initial data for
the electromagnetic field
${\boldsymbol{u}}=({\boldsymbol{E}},{\boldsymbol{H}})$ at $t=0$. Consider for
example the WKB form [23, 21, 49]:

(8)
${\boldsymbol{u}}({\boldsymbol{x}},0)={\boldsymbol{u}}_{\varepsilon}^{\text{WKB}}({\boldsymbol{x}})=\Re\mathfrak{e}\\{{\boldsymbol{u}}_{I}({\boldsymbol{x}})\operatorname{e}^{\frac{\mathrm{i}}{\varepsilon}S_{I}({\boldsymbol{x}})}\\}\,,$

where $\smash{{\boldsymbol{u}}_{I}}$ is a square integrable function on
$\mathbb{R}^{3}$; $\smash{S_{I}}$ is an integrable scalar function on
$\mathbb{R}^{3}$ (at least locally), as well as its derivative
$\smash{{\boldsymbol{\nabla}}_{\boldsymbol{x}}S_{I}}$; and $\Re\mathfrak{e}$
stands for the real part. The small parameter $\varepsilon$ represents the
typical wavelength of oscillations of these data. These oscillations are
inherited by the actual electromagnetic field
${\boldsymbol{u}}({\boldsymbol{x}},t)$ satisfying (4) at all times, which
prevents it from converging nicely as the wavelength $\varepsilon$ gets
smaller and smaller in the high-frequency limit. A way to tackle this
shortcoming is to consider quadratic observables of ${\boldsymbol{u}}$
instead, as emphasized in _e.g._ [1, 8, 9, 10, 26, 27, 43, 46] for different
types of waves. In particular, the Wigner transform [23, 35, 37, 53] of the
solutions of Eq. (4) shall be considered since its evolution in the limit
$\varepsilon\to 0$ can be derived explicitly. It provides a phase space
description of how the associated energy density propagates in this very
limit. The advantages of this approach over more classical methods such as the
WKB method, which is suggested by the type of data (8), are discussed in [21,
49] for example. Notably, the Wigner transform needs much lower regularity
assumptions on $\smash{{\boldsymbol{u}}_{I}}$ and $\smash{S_{I}}$ than in the
WKB method, and evades the singularities (focal points or caustics) that can
develop in finite times with the latter method.

Let $\vartheta\in[0,1]$. For temperate distributions
${\boldsymbol{u}},{\boldsymbol{v}}$ defined on $\mathbb{R}^{m}$, the Wigner
transform is:

(9)
${\boldsymbol{W}}_{\varepsilon}^{\vartheta}[{\boldsymbol{u}},{\boldsymbol{v}}]({\boldsymbol{x}},{\boldsymbol{k}})=\frac{1}{(2\pi)^{m}}\int_{\mathbb{R}^{m}}\operatorname{e}^{\mathrm{i}{\boldsymbol{k}}\cdot{\boldsymbol{x}}}{\boldsymbol{u}}\left({\boldsymbol{x}}-\varepsilon(1-\vartheta){\boldsymbol{y}}\right){{\boldsymbol{v}}}^{*}\left({\boldsymbol{x}}+\varepsilon\vartheta{\boldsymbol{y}}\right)\,{\mathrm{d}}{\boldsymbol{y}}\,.$

We denote by
${\boldsymbol{W}}_{\varepsilon}^{\vartheta}[{\boldsymbol{u}}]:={\boldsymbol{W}}_{\varepsilon}^{\vartheta}[{\boldsymbol{u}},{\boldsymbol{u}}]$
the self Wigner transform of ${\boldsymbol{u}}$. Now let
$({\boldsymbol{u}}_{\varepsilon})$ be a bounded sequence in
$\smash{L^{2}(\mathbb{R}^{m})}$, the set of square integrable functions
defined on $\mathbb{R}^{m}$. Then it can be shown that, up to extracting a
subsequence if need be,
${\boldsymbol{W}}_{\varepsilon}^{\vartheta}[{\boldsymbol{u}}_{\varepsilon}]$
has a weak limit (in the sense of temperate distribution) as $\varepsilon\to
0$ which is independent of the choice of $\vartheta$ in Eq. (9). Also if
${\boldsymbol{W}}[{\boldsymbol{u}}_{\varepsilon}]$ is such a limit, it is a
non-negative, Hermitian matrix-valued measure. These results are detailed in
_e.g._ [23, 53]. Here the notation of [23, p. 330] is used for the limit
${\boldsymbol{W}}$ (independent of $\varepsilon$) of the family
$\smash{({\boldsymbol{u}}_{\varepsilon})}$ (dependent of $\varepsilon$) but
clearly the measure ${\boldsymbol{W}}[{\boldsymbol{u}}_{\varepsilon}]$ is
independent of $\varepsilon$. The scalar Wigner measure is the trace of the
latter
$W[{\boldsymbol{u}}_{\varepsilon}]=\operatorname{Tr}{\boldsymbol{W}}[{\boldsymbol{u}}_{\varepsilon}]$
and can be directly related to the energy density of the sequence
$({\boldsymbol{u}}_{\varepsilon})$ in the high-frequency limit. Taking the
Wigner transform of the sequence of WKB data (8) for example, one obtains:

$\begin{split}{\boldsymbol{W}}[{\boldsymbol{u}}_{\varepsilon}^{\text{WKB}}]({\boldsymbol{x}},{\boldsymbol{k}})&={\boldsymbol{u}}_{I}({\boldsymbol{x}}){{\boldsymbol{u}}}^{*}_{I}({\boldsymbol{x}})\delta({\boldsymbol{k}}-{\boldsymbol{\nabla}}_{\boldsymbol{x}}S_{I}({\boldsymbol{x}}))\,,\\\
W[{\boldsymbol{u}}_{\varepsilon}^{\text{WKB}}]({\boldsymbol{x}},{\boldsymbol{k}})&=|{\boldsymbol{u}}_{I}({\boldsymbol{x}})|^{2}\delta({\boldsymbol{k}}-{\boldsymbol{\nabla}}_{\boldsymbol{x}}S_{I}({\boldsymbol{x}}))\,.\end{split}$

If these data are propagated by the wave equation in $\mathbb{R}^{3}$ with
constant speed $c$:

$\begin{split}\frac{1}{c^{2}}\frac{\partial^{2}{\boldsymbol{u}}_{\varepsilon}}{\partial
t^{2}}-\Delta{\boldsymbol{u}}_{\varepsilon}={\bf
0}\,,&\quad{\boldsymbol{x}}\in\mathbb{R}^{3}\,,\;t>0\,,\\\
{\boldsymbol{u}}_{\varepsilon}({\boldsymbol{x}},0)={\bf
0}\,,\;\frac{\partial{\boldsymbol{u}}_{\varepsilon}}{\partial
t}({\boldsymbol{x}},0)={\boldsymbol{u}}_{\varepsilon}^{\text{WKB}}({\boldsymbol{x}})\,,&\quad{\boldsymbol{x}}\in\mathbb{R}^{3}\,,\end{split}$

then for $t>0$ Kirchhoff’s formula yields [17, §2.4]:

${\boldsymbol{u}}_{\varepsilon}({\boldsymbol{x}},t)=\int_{S^{2}}t{\boldsymbol{u}}_{\varepsilon}^{\text{WKB}}({\boldsymbol{x}}-c{\hat{\boldsymbol{p}}}t)\frac{{\mathrm{d}}\Omega({\hat{\boldsymbol{p}}})}{4\pi}\,,$

where $S^{2}$ is the unit sphere of $\mathbb{R}^{3}$ centered at ${\bf 0}$,
and ${\mathrm{d}}\Omega$ is its surface differential element. The scalar
Wigner measure of ${\boldsymbol{u}}_{\varepsilon}(\cdot,t)$ reads [23, §3.3]:

$W[{\boldsymbol{u}}_{\varepsilon}]({\boldsymbol{x}},{\boldsymbol{k}},t)=\frac{1}{2}\sum_{\alpha=\pm}|\smash{{\boldsymbol{u}}_{I}({\boldsymbol{x}}+\alpha
c{\hat{\boldsymbol{k}}}t)}|^{2}\delta({\boldsymbol{k}}-{\boldsymbol{\nabla}}_{\boldsymbol{x}}S_{I}({\boldsymbol{x}}+\alpha
c{\hat{\boldsymbol{k}}}t))\,,$

with the usual notation
${\hat{\boldsymbol{k}}}=\smash{\frac{{\boldsymbol{k}}}{\left|{\boldsymbol{k}}\right|}}$
for ${\boldsymbol{k}}\neq{\bf 0}$. It provides a characterization of the
evolution of the associated energy density both in space and direction.
Because of the assumed low regularity of $\smash{{\boldsymbol{u}}_{I}}$ and
$\smash{S_{I}}$, a classical WKB method could not be used with these data.
Note that more general data can be considered as well as shown in [23].

###  2.3. Main results

We derive the Wigner measure ${\boldsymbol{W}}$ of the high-frequency
electromagnetic field ${\boldsymbol{u}}$ which obeys Maxwell’s $6\times 6$
system (4) when the optical response ${\boldsymbol{K}}_{0}$ exhibits in
addition random fluctuations. That is:

(10) $\frac{\partial}{\partial
t}\left({\boldsymbol{K}}({\boldsymbol{x}}){\boldsymbol{u}}({\boldsymbol{x}},t)+\int_{0}^{t}{\boldsymbol{K}}_{d}({\boldsymbol{x}},\tau){\boldsymbol{u}}({\boldsymbol{x}},t-\tau){\mathrm{d}}\tau\right)+\boldsymbol{M}({\boldsymbol{\nabla}}_{\boldsymbol{x}}){\boldsymbol{u}}({\boldsymbol{x}},t)={\bf
0}$

for ${\boldsymbol{x}},t\in{\mathcal{O}}\times\mathbb{R}_{+}^{*}$ in some open
domain ${\mathcal{O}}\subseteq\smash{\mathbb{R}^{3}}$ and for highly
oscillating initial conditions at $t=0$ characterized by the small scale
$\varepsilon\ll 1$. Here:

(11) $\boldsymbol{M}({\boldsymbol{k}})=\begin{bmatrix}{\bf
0}&-{\boldsymbol{\Omega}}({\boldsymbol{k}})\\\
{\boldsymbol{\Omega}}({\boldsymbol{k}})&{\bf 0}\end{bmatrix}$

is the symmetric Maxwell operator for
${\boldsymbol{\Omega}}({\boldsymbol{k}})$ being the skew-symmetric matrix such
that ${\boldsymbol{\Omega}}({\boldsymbol{k}}){\bf
a}={\boldsymbol{k}}\times{\bf a}$ for any vector ${\bf a}\in\mathbb{C}^{3}$.
Also the actual optical response ${\boldsymbol{K}}({\boldsymbol{x}})$ reads:

${\boldsymbol{K}}({\boldsymbol{x}})={\boldsymbol{K}}_{0}({\boldsymbol{x}})\left[{\boldsymbol{I}}+\sigma{\boldsymbol{V}}\left(\frac{{\boldsymbol{x}}}{\ell_{c}}\right)\right]\,,$

where the dimensionless random matrix ${\boldsymbol{V}}$ represents random
fluctuations of the optical response ${\boldsymbol{K}}_{0}$ of the background
medium, and is such that
${\boldsymbol{K}}_{0}{\boldsymbol{V}}={{\boldsymbol{V}}}^{*}{\boldsymbol{K}}_{0}$
to preserve the Hermiticity of the actual optical response ${\boldsymbol{K}}$.
It is assumed that
$({\boldsymbol{V}}({\boldsymbol{y}}),\,{\boldsymbol{y}}\in\mathbb{R}^{3})$ is
a second-order, mean-square homogeneous (spatially stationary) random field
with zero mean
${\mathbb{E}}\left\\{{\boldsymbol{V}}({\boldsymbol{y}})\right\\}={\bf 0}$ and
integrable fourth-order autocorrelation tensor ${\boldsymbol{R}}$ such that:

${\boldsymbol{R}}({\boldsymbol{y}}-{\boldsymbol{y}}^{\prime}):={\mathbb{E}}\left\\{{\boldsymbol{V}}({\boldsymbol{y}})\otimes{\boldsymbol{V}}({\boldsymbol{y}}^{\prime})\right\\}=\int_{\mathbb{R}^{3}}\operatorname{e}^{\mathrm{i}{\boldsymbol{k}}\cdot({\boldsymbol{y}}-{\boldsymbol{y}}^{\prime})}\widehat{{\boldsymbol{R}}}({\boldsymbol{k}}){\mathrm{d}}{\boldsymbol{k}}\,,$

where ${\mathbb{E}}\left\\{\cdot\right\\}$ stands for the mathematical
expectation, or average. The second equality above stems from Bochner’s
theorem and the fact that
${\boldsymbol{y}}\mapsto{\boldsymbol{R}}({\boldsymbol{y}})$ is positive semi-
definite. The power spectral density of that random field is the matrix:

(12)
$\widehat{{\boldsymbol{R}}}({\boldsymbol{k}})=\frac{1}{(2\pi)^{3}}\int_{\mathbb{R}^{3}}\operatorname{e}^{-\mathrm{i}{\boldsymbol{k}}\cdot{\boldsymbol{y}}}{\boldsymbol{R}}({\boldsymbol{y}}){\mathrm{d}}{\boldsymbol{k}}$

which is positive. The autocorrelation tensor is normalized such that
$\int_{\mathbb{R}^{3}}{\boldsymbol{R}}({\boldsymbol{y}}){\mathrm{d}}{\boldsymbol{y}}=\mathrm{O}(1)$
and ${\boldsymbol{R}}({\bf 0})=\mathrm{O}(1)$. The length scale $\ell_{c}$ is
the correlation length and the dimensionless scalar $\sigma\geq 0$ quantifies
the amplitude of the fluctuations. These parameters are of the same order for
all correlations of the fluctuations in the proposed model.

At first, the random fluctuations of the optical response are disregarded,
setting $\sigma=0$. Let:

(13)
$\begin{split}{\boldsymbol{L}}_{0}({\boldsymbol{x}},{\boldsymbol{k}})&={\boldsymbol{K}}_{0}^{-1}({\boldsymbol{x}})\boldsymbol{M}({\boldsymbol{k}})\,,\\\
{\boldsymbol{L}}_{1}({\boldsymbol{x}},\omega)&=\mathrm{i}\omega{\boldsymbol{K}}_{0}^{-1}({\boldsymbol{x}})\widehat{{\boldsymbol{K}}}_{d}({\boldsymbol{x}},\omega)\,,\end{split}$

where $\widehat{{\boldsymbol{K}}}_{d}({\boldsymbol{x}},\omega)$ is the Fourier
transform of ${\boldsymbol{K}}_{d}({\boldsymbol{x}},t)$ with respect to $t$.
${\boldsymbol{L}}_{0}$ being symmetric for the scalar product
$\left\langle\cdot\right\rangle$, let us introduce the eigen-expansion:

(14)
${\boldsymbol{L}}_{0}=\sum_{\alpha}\omega_{\alpha}\boldsymbol{\Pi}_{\alpha}$

where $\boldsymbol{\Pi}_{\alpha}$ is the projector on the $\alpha$-th eigen-
subspace with associated real eigenvalue $\omega_{\alpha}$, such that
$\sum_{\alpha}\boldsymbol{\Pi}_{\alpha}={\boldsymbol{I}}$ and
${\boldsymbol{K}}_{0}\boldsymbol{\Pi}_{\alpha}={\boldsymbol{\Pi}}^{*}_{\alpha}{\boldsymbol{K}}_{0}$.
The right eigenvectors
${\boldsymbol{b}}_{\alpha}=({\boldsymbol{b}}_{\alpha}^{1},\dots\boldsymbol{b}_{\alpha}^{A})$
are such that
$\left\langle{\boldsymbol{b}}_{\alpha},{\boldsymbol{b}}_{\beta}\right\rangle={({\boldsymbol{K}}_{0}{\boldsymbol{b}}_{\alpha})}^{*}{\boldsymbol{b}}_{\beta}=\delta_{\alpha\beta}{\boldsymbol{I}}_{A}$,
where $A$ is the order of multiplicity of the $\alpha$-th eigenvalue and
${\boldsymbol{I}}_{A}$ is the $A\times A$ identity matrix; also $\sum A=6$.
Letting
${\boldsymbol{c}}_{\alpha}={\boldsymbol{K}}_{0}{\boldsymbol{b}}_{\alpha}$ be
the left eigenvectors, one has therefore:

(15)
$\begin{split}{\boldsymbol{L}}_{0}{\boldsymbol{b}}_{\alpha}&=\omega_{\alpha}{\boldsymbol{b}}_{\alpha}\,,\\\
{{\boldsymbol{c}}}^{*}_{\alpha}{\boldsymbol{L}}_{0}&=\omega_{\alpha}{{\boldsymbol{c}}}^{*}_{\alpha}\,,\\\
\boldsymbol{\Pi}_{\alpha}&={\boldsymbol{b}}_{\alpha}{{\boldsymbol{c}}}^{*}_{\alpha}\,.\end{split}$

Then it is shown that:

(16)
${\boldsymbol{W}}=\sum_{\alpha}\delta(\omega+\omega_{\alpha}){\boldsymbol{W}}_{\alpha}\,,$

where
${\boldsymbol{W}}_{\alpha}=\boldsymbol{\Pi}_{\alpha}{\boldsymbol{W}}={\boldsymbol{W}}{\boldsymbol{\Pi}}^{*}_{\alpha}=\boldsymbol{\Pi}_{\alpha}{\boldsymbol{W}}{\boldsymbol{\Pi}}^{*}_{\alpha}$
are $6\times 6$ matrix-valued measures for each mode of propagation $\alpha$.

Alternatively ${\boldsymbol{W}}$ may be written:

(17)
${\boldsymbol{W}}=\sum_{\alpha}\delta(\omega+\omega_{\alpha}){\boldsymbol{b}}_{\alpha}{\bf
w}_{\alpha}{{\boldsymbol{b}}}^{*}_{\alpha}\,,$

where ${\bf
w}_{\alpha}={{\boldsymbol{c}}}^{*}_{\alpha}{\boldsymbol{W}}_{\alpha}{\boldsymbol{c}}_{\alpha}$
are $A\times A$ matrix-valued measures. In Sect. 3, the latter are shown to
satisfy the transport equations:

(18) $\partial_{t}{\bf w}_{\alpha}+\\{\omega_{\alpha},{\bf
w}_{\alpha}\\}+{\boldsymbol{\ell}}_{\alpha}{\bf w}_{\alpha}+{\bf
w}_{\alpha}{{\boldsymbol{\ell}}}^{*}_{\alpha}+{\boldsymbol{n}}_{\alpha}{\bf
w}_{\alpha}-{\bf w}_{\alpha}{\boldsymbol{n}}_{\alpha}={\bf 0}\,,$

where $\\{\omega_{\alpha},{\bf
w}_{\alpha}\\}={\boldsymbol{\nabla}}_{\boldsymbol{k}}\omega_{\alpha}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{x}}{\bf
w}_{\alpha}-{\boldsymbol{\nabla}}_{\boldsymbol{x}}\omega_{\alpha}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}{\bf
w}_{\alpha}$ stands for the usual Poisson bracket, ${\boldsymbol{n}}_{\alpha}$
is a $A\times A$ skew-symmetric matrix given below by Eq. (43), and
$\smash{{\boldsymbol{\ell}}_{\alpha}}=\smash{{{\boldsymbol{c}}}^{*}_{\alpha}{\boldsymbol{L}}_{1}{\boldsymbol{b}}_{\alpha}}$.
Taking the trace of Eq. (18) yields:

$\partial_{t}w_{\alpha}+\\{\omega_{\alpha},w_{\alpha}\\}+2\operatorname{Tr}({\boldsymbol{\ell}}_{\alpha}^{s}{\bf
w}_{\alpha})=0\,,$

where $w_{\alpha}=\operatorname{Tr}{\bf
w}_{\alpha}=\operatorname{Tr}({\boldsymbol{K}}_{0}{\boldsymbol{W}}_{\alpha})$,
and
$\smash{{\boldsymbol{\ell}}_{\alpha}^{s}}=\smash{\frac{1}{2}({\boldsymbol{\ell}}_{\alpha}+{{\boldsymbol{\ell}}}^{*}_{\alpha})}$
is the symmetric part of $\smash{{\boldsymbol{\ell}}_{\alpha}}$. The high-
frequency energy density and Poynting vector (7) are then:

(19)
$\begin{split}{\mathcal{E}}({\boldsymbol{x}},t)&=\frac{1}{2}\int_{\mathbb{R}^{3}}\operatorname{Tr}({\boldsymbol{K}}_{0}({\boldsymbol{x}}){\boldsymbol{W}}({\boldsymbol{x}},{\boldsymbol{k}},t)){\mathrm{d}}{\boldsymbol{k}}=\frac{1}{2}\sum_{\alpha}\int_{\mathbb{R}^{3}}w_{\alpha}({\boldsymbol{x}},{\boldsymbol{k}},t){\mathrm{d}}{\boldsymbol{k}}\,,\\\
\boldsymbol{{\mathcal{F}}}({\boldsymbol{x}},t)&=\frac{1}{2}\int_{\mathbb{R}^{3}}\operatorname{Tr}({\boldsymbol{\nabla}}_{\boldsymbol{k}}\boldsymbol{M}({\boldsymbol{k}}){\boldsymbol{W}}({\boldsymbol{x}},{\boldsymbol{k}},t)){\mathrm{d}}{\boldsymbol{k}}=\frac{1}{2}\sum_{\alpha}\int_{\mathbb{R}^{3}}w_{\alpha}({\boldsymbol{x}},{\boldsymbol{k}},t){\boldsymbol{\nabla}}_{\boldsymbol{k}}\omega_{\alpha}{\mathrm{d}}{\boldsymbol{k}}\,.\end{split}$

The evolution of the energy density and flow in phase space are thus described
by $\sum_{\alpha}w_{\alpha}$ which can also select the different modes of
propagation. We observe however from Eq. (13) with Eq. (11), that
$\smash{\omega_{0}}=0$ ($\alpha=\text{"0"}$) is always an eigenvalue of
$\smash{{\boldsymbol{L}}_{0}}$ with multiplicity $A=2$ and eigenvectors
$\smash{({\hat{\boldsymbol{k}}},{\bf 0})}$ and $\smash{({\bf
0},{\hat{\boldsymbol{k}}})}$ for ${\boldsymbol{k}}\neq{\bf 0}$. This mode is
non propagative, though.

Secondly, we analyze the influence of random fluctuations of the optical
response $\smash{{\boldsymbol{K}}_{0}}$ of the background medium by setting
$\sigma>0$. We consider the scattering regime where $\ell_{c}$ is of the order
of the typical wavelength $\lambda$, which is the scale of variation of the
strongly oscillating initial data, _i.e._ it is small with respect to the
typical propagation distance $L$:

(20) $\frac{\ell_{c}}{L}\approx\frac{\lambda}{L}=\varepsilon\,.$

Besides, the amplitude of fluctuations is also small with the scaling:

(21) $\sigma^{2}\approx\frac{\ell_{c}}{L}\,.$

For a large distance of propagation $L$ the electromagnetic waves are multiply
scattered by the fluctuations of the background medium and their energy is
spread over many directions of propagation. Considering these fluctuations, it
is shown in Sect. 4 that the transport equations (18) are modified to the
coupled radiative transfer equations:

(22) $\partial_{t}{\bf w}_{\alpha}+\\{\omega_{\alpha},{\bf
w}_{\alpha}\\}+({\boldsymbol{\ell}}_{\alpha}+{\boldsymbol{\Sigma}}_{\alpha}){\bf
w}_{\alpha}+{\bf
w}_{\alpha}{({\boldsymbol{\ell}}_{\alpha}+{\boldsymbol{\Sigma}}_{\alpha})}^{*}+{\boldsymbol{n}}_{\alpha}{\bf
w}_{\alpha}-{\bf w}_{\alpha}{\boldsymbol{n}}_{\alpha}\\\
=\sum_{\beta}\int_{\mathbb{R}^{3}}\delta(\omega_{\beta}({\boldsymbol{x}},{\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{x}},{\boldsymbol{k}})){\boldsymbol{\sigma}}_{\alpha\beta}({\boldsymbol{x}},{\boldsymbol{k}},{\boldsymbol{p}}):{\bf
w}_{\beta}({\boldsymbol{x}},{\boldsymbol{p}},t){\mathrm{d}}{\boldsymbol{p}}\,,$

where ${\boldsymbol{\sigma}}_{\alpha\beta}$ is a linear operator called
differential scattering cross-section, and ${\boldsymbol{\Sigma}}_{\alpha}$ is
a $A\times A$ matrix called total scattering cross-section. They are
explicitly given in terms of the power spectral density matrix (12) of the
random fluctuations ${\boldsymbol{V}}$ (see Eq. (70) and Eq. (71),
respectively). The $A\times A\times B\times B$ tensor
${\boldsymbol{\sigma}}_{\alpha\beta}({\boldsymbol{x}},{\boldsymbol{k}},{\boldsymbol{p}})$
describes how, in the multiple scattering process, the $B\times B$ matrix
${\bf w}_{\beta}$ for the mode $\beta$ with multiplicity $B$ in the direction
${\boldsymbol{p}}$ is locally converted to the $A\times A$ matrix ${\bf
w}_{\alpha}$ for the mode $\alpha$ with multiplicity $A$ in the direction
${\boldsymbol{k}}$:

$[{\boldsymbol{\sigma}}_{\alpha\beta}({\boldsymbol{x}},{\boldsymbol{k}},{\boldsymbol{p}}):{\bf
w}_{\beta}({\boldsymbol{x}},{\boldsymbol{p}},t)]_{aa^{\prime}}=\sum_{1\leq
b,b^{\prime}\leq
B}[{\boldsymbol{\sigma}}_{\alpha\beta}({\boldsymbol{x}},{\boldsymbol{k}},{\boldsymbol{p}})]_{aa^{\prime}bb^{\prime}}[{\bf
w}_{\beta}({\boldsymbol{x}},{\boldsymbol{p}},t)]_{bb^{\prime}}$

for $1\leq a,a^{\prime}\leq A$. The $A\times A$ matrix
${\boldsymbol{\Sigma}}_{\alpha}({\boldsymbol{x}},{\boldsymbol{k}})$ describes
how the matrix ${\bf w}_{\alpha}$ in the direction ${\boldsymbol{k}}$ is
locally converted to all other modes $\beta$ and directions. Lastly, the Dirac
measure in the integral indicates that scattering occurs with possible mode
conversion whenever for some directions ${\boldsymbol{p}}$ and
${\boldsymbol{k}}$ one has
$\omega_{\beta}({\boldsymbol{x}},{\boldsymbol{p}})=\omega_{\alpha}({\boldsymbol{x}},{\boldsymbol{k}})$.
This includes the situation when no mode conversion occurs, $\alpha=\beta$,
and $\omega_{\alpha}$ has locally the same value for possibly two different
directions. For an isotropic, lossless background medium with
${\boldsymbol{K}}_{0}({\boldsymbol{x}})=\operatorname{diag}(\epsilon_{0}({\boldsymbol{x}}){\boldsymbol{I}},\mu_{0}({\boldsymbol{x}}){\boldsymbol{I}})$
for example, one has three modes of propagation each of multiplicity $2$ with
eigenvalues $\omega_{1}=0$,
$\omega_{2}({\boldsymbol{x}},{\boldsymbol{k}})=+c_{0}({\boldsymbol{x}})|{\boldsymbol{k}}|$,
and
$\omega_{2}({\boldsymbol{x}},{\boldsymbol{k}})=-c_{0}({\boldsymbol{x}})|{\boldsymbol{k}}|$
where $c_{0}=1/\sqrt{\epsilon_{0}\mu_{0}}$ [46]. Our results (18) and (22)
generalize [46] to general bianisotropic dielectric media, and [18, 19] to the
situation where these media are in addition heterogeneous and dissipative.

We now turn to the detailed analysis of the high-frequency limit
$\varepsilon\to 0$ of Eq. (10) and how these results are derived. For that
purpose, we use a classical pseudo-differential calculus, which also
contributes to significantly simplify these derivations.

##  3\. Transport of high-frequency electromagnetic waves in bianisotropic
media

In this section the transport equations (18) for bianisotropic dielectric
media are first derived. They describe the evolution of the high-frequency
electromagnetic energy density without random fluctuations of the optical
response. As illustrated in Sect. 2.2, the high-frequency limit
$\varepsilon\rightarrow 0$ in the previous setting shall be derived for
quadratic observables of the electromagnetic field ${\boldsymbol{u}}$. More
particularly, we introduce the spatio-temporal Wigner transform of that field
and its high-frequency limit, _i.e._ its Wigner measure, as in _e.g._ [2].
This is because the spatial and temporal scales in Maxwell’s system (10) play
symmetric roles, and their oscillations should be accounted for altogether.
The main objective here is to outline the (formal) mathematical tools used for
the derivation of the Wigner measure’s properties. They will prove useful in
the subsequent Sect. 4 where the effects of the random fluctuations of the
optical response will be considered. The analysis is derived from [23], where
first-order hyperbolic systems with constant and slowly varying coefficients
are addressed, [1], where arbitrary order hyperbolic systems with slowly
varying coefficients are addressed, and [2], where time-varying random media
and consistent pseudo-differential calculus and Wigner transforms with a
single quantization are used. In this way we first recall in Sect. 3.1 how
Wigner transforms are linked to semi-classical pseudo-differential operators
[37, 53]. Then writing (10) as a semi-classical operator applied to a rescaled
version of the electromagnetic field ${\boldsymbol{u}}$ in Sect. 3.2, the
dispersion properties of the Wigner measure (related to the Stokes parameters
of the electromagnetic waves) are derived in Sect. 3.3, and the transport
equations (18) are obtained Sect. 3.4.

###  3.1. Semi-classical operators and Wigner measure

From now on let us introduce the space-time variable
${\boldsymbol{s}}=({\boldsymbol{x}},t)\in{\mathcal{O}}\times\mathbb{R}$ and
its dual variable
${\boldsymbol{\xi}}=({\boldsymbol{k}},\omega)\in\mathbb{R}^{4}$ in the wave
vector-frequency Fourier domain. Let $\boldsymbol{L}$ be a smooth, compactly
supported real matrix-valued function of both the space-time variable
${\boldsymbol{s}}$ and impulse variable ${\boldsymbol{\xi}}$. For a vector
field ${\boldsymbol{u}}$ in $\smash{L^{2}(\mathbb{R}^{4})}$ endowed with the
scalar product
$\smash{\left({\boldsymbol{u}},{\boldsymbol{v}}\right)_{L^{2}}}=\smash{\int_{\mathbb{R}^{4}}{\boldsymbol{u}}({\boldsymbol{s}})\cdot\overline{{\boldsymbol{v}}}({\boldsymbol{s}})\,{\mathrm{d}}{\boldsymbol{s}}}$
where $\smash{\overline{{\boldsymbol{v}}}}$ stands for complex conjugation,
consider the (semiclassical) pseudo-differential operator:

(23) $\boldsymbol{L}^{\vartheta}({\boldsymbol{s}},\varepsilon{\bf
D}){\boldsymbol{u}}({\boldsymbol{s}})=\frac{1}{(2\pi)^{4}}\int_{\mathbb{R}^{4}\times\mathbb{R}^{4}}\operatorname{e}^{\mathrm{i}{\boldsymbol{\xi}}\cdot({\boldsymbol{s}}-{\boldsymbol{\tau}})}\boldsymbol{L}((1-\vartheta){\boldsymbol{s}}+\vartheta{\boldsymbol{\tau}},\varepsilon{\boldsymbol{\xi}}){\boldsymbol{u}}({\boldsymbol{\tau}})\,{\mathrm{d}}{\boldsymbol{\tau}}{\mathrm{d}}{\boldsymbol{\xi}}\,,$

for $\vartheta\in[0,1]$. This parameter defines the so-called quantization of
the operator $\boldsymbol{L}^{\vartheta}$. The case $\vartheta=0$ corresponds
to the standard quantization. It is simply denoted by
$\smash{\boldsymbol{L}({\boldsymbol{s}},\varepsilon{\bf D})}$ such that:

(24) $\boldsymbol{L}({\boldsymbol{s}},\varepsilon{\bf
D}){\boldsymbol{u}}({\boldsymbol{s}})=\frac{1}{(2\pi)^{4}}\int_{\mathbb{R}^{4}}\operatorname{e}^{\mathrm{i}{\boldsymbol{\xi}}\cdot{\boldsymbol{s}}}\boldsymbol{L}({\boldsymbol{s}},\varepsilon{\boldsymbol{\xi}})\widehat{{\boldsymbol{u}}}({\boldsymbol{\xi}})\,{\mathrm{d}}{\boldsymbol{\xi}}\,,$

where:

(25)
$\widehat{{\boldsymbol{u}}}({\boldsymbol{\xi}})=\int_{\mathbb{R}^{4}}\operatorname{e}^{-\mathrm{i}{\boldsymbol{\xi}}\cdot{\boldsymbol{s}}}{\boldsymbol{u}}({\boldsymbol{s}}){\mathrm{d}}{\boldsymbol{s}}$

stands for the space-time Fourier transform of
${\boldsymbol{u}}({\boldsymbol{s}})$. The case $\smash{\vartheta=1/2}$
corresponds to the Weyl quantization, which is usually denoted by
$\smash{\boldsymbol{L}^{W}({\boldsymbol{s}},\varepsilon{\bf D})}$.

Now let ${\boldsymbol{u}},{\boldsymbol{v}}$ be temperate distributions defined
on $\smash{\mathbb{R}^{4}}$ and let
${\boldsymbol{W}}_{\varepsilon}^{\vartheta}[{\boldsymbol{u}},{\boldsymbol{v}}]({\boldsymbol{s}},{\boldsymbol{\xi}})$
be their Wigner transform (9) (with $m=4$ here). Then one has the trace
formula [23]:

(26) $\left(\boldsymbol{L}^{\vartheta}({\boldsymbol{s}},\varepsilon{\bf
D}){\boldsymbol{u}},{\boldsymbol{v}}\right)_{L^{2}}=\int_{\mathbb{R}^{4}\times\mathbb{R}^{4}}\operatorname{Tr}\boldsymbol{L}({\boldsymbol{s}},{\boldsymbol{\xi}}){\boldsymbol{W}}_{\varepsilon}^{\vartheta}[{\boldsymbol{u}},{\boldsymbol{v}}]({\mathrm{d}}{\boldsymbol{s}},{\mathrm{d}}{\boldsymbol{\xi}})\,.$

For a sequence $\smash{({\boldsymbol{u}}_{\varepsilon})}$ uniformly bounded in
$\smash{L^{2}(\mathbb{R}^{4})}$, one can establish in particular that there
exists a positive measure ${\boldsymbol{W}}[{\boldsymbol{u}}_{\varepsilon}]$
such that, up to extracting a subsequence if need be (see _e.g._ [53, Theorem
5.2]):

(27) $\lim_{\varepsilon\rightarrow
0}\left(\boldsymbol{L}^{\vartheta}({\boldsymbol{s}},\varepsilon{\bf
D}){\boldsymbol{u}}_{\varepsilon},{\boldsymbol{u}}_{\varepsilon}\right)_{L^{2}}=\int_{\mathbb{R}^{4}\times\mathbb{R}^{4}}\operatorname{Tr}\boldsymbol{L}({\boldsymbol{s}},{\boldsymbol{\xi}}){\boldsymbol{W}}[{\boldsymbol{u}}_{\varepsilon}]({\mathrm{d}}{\boldsymbol{s}},{\mathrm{d}}{\boldsymbol{\xi}})\,,\quad\forall\boldsymbol{L}\,,$

independently of the quantization $\vartheta$. We thus recover the Wigner
measure ${\boldsymbol{W}}[{\boldsymbol{u}}_{\varepsilon}]$ of
$\smash{({\boldsymbol{u}}_{\varepsilon})}$ invoked in Sect. 2.2, which can
also be interpreted as the weak limit of its self Wigner transform
$\smash{{\boldsymbol{W}}_{\varepsilon}^{\vartheta}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{u}}_{\varepsilon}]:={\boldsymbol{W}}_{\varepsilon}^{\vartheta}[{\boldsymbol{u}}_{\varepsilon}]}$.
It describes the limit energy density of the sequence
$\smash{({\boldsymbol{u}}_{\varepsilon})}$ in the phase space
$\smash{\mathbb{R}^{4}_{\boldsymbol{s}}\times\mathbb{R}^{4}_{\boldsymbol{\xi}}}$.
The observable $\boldsymbol{L}({\boldsymbol{s}},{\boldsymbol{\xi}})$ is used
to select any quadratic observable or quantity of interest associated to this
energy: the kinetic energy, or the free energy, or the Poynting vector as in
Eq. (19) for example.

Next, we recall some formulas of pseudo-differential calculus which involve
Wigner transforms and semiclassical pseudo-differential operators. In [8] for
example, various formal rules of pseudo-differential calculus were given for
quantities like
$\smash{{\boldsymbol{W}}_{\varepsilon}^{\frac{1}{2}}[\boldsymbol{L}({\boldsymbol{x}},\varepsilon{\bf
D}_{\boldsymbol{x}}){\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}]}$,
mixing the quantization chosen for the Wigner transform (9)
($\vartheta=\frac{1}{2}$) and the standard quantization of the semiclassical
operator $\boldsymbol{L}({\boldsymbol{x}},\varepsilon{\bf
D}_{\boldsymbol{x}})$ ($\vartheta=0$). These formulas have been revisited in
[2] using the same standard quantization for both the Wigner transforms and
pseudo-differential operators, which seems more natural in view of the trace
formula (26). They also contribute to significantly simplify the calculations
below. They are recalled here, extending slightly the results of [2] for
scalar-valued Wigner transforms to matrix-valued Wigner transforms. The proofs
are straightforward applications of the proofs in [2, Appendix B] for the
scalar case. For the chosen standard quantization $\vartheta=0$, the Wigner
transform (9) reads:

(28)
${\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}]({\boldsymbol{s}},{\boldsymbol{\xi}}):=\frac{1}{(2\pi)^{4}}\int_{\mathbb{R}^{4}}\operatorname{e}^{\mathrm{i}{\boldsymbol{\xi}}\cdot{\boldsymbol{\tau}}}{\boldsymbol{u}}_{\varepsilon}({\boldsymbol{s}}-\varepsilon{\boldsymbol{\tau}}){{\boldsymbol{v}}}^{*}_{\varepsilon}({\boldsymbol{s}})\,{\mathrm{d}}{\boldsymbol{\tau}}\,,$

dropping the superscript $\vartheta$ from now on. The self Wigner transform
${\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{u}}_{\varepsilon}]$
is denoted by
${\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]$, as
implicitly done in Eq. (27). Let
$\boldsymbol{L}({\boldsymbol{s}},{\boldsymbol{\xi}})$ be a smooth matrix-
valued function defined on
$\smash{\mathbb{R}^{4}_{\boldsymbol{s}}\times\mathbb{R}_{\boldsymbol{\xi}}^{4}}$
of which derivatives increase only slowly,111That is, $\boldsymbol{L}$ satisfy
for some $m\geq 0$ [23, Prop. 1.8]:
$|\partial_{\boldsymbol{s}}^{\alpha}\partial_{\boldsymbol{\xi}}^{\beta}\boldsymbol{L}({\boldsymbol{s}},{\boldsymbol{\xi}})|\leq
C_{\alpha,\beta}(1+|{\boldsymbol{\xi}}|)^{m}\,,\quad\forall\alpha,\beta\in\mathbb{N}^{4}\,.$
and recall the notation of Eq. (24) for the operator
$\smash{\boldsymbol{L}({\boldsymbol{s}},\varepsilon{\bf D})}$ with
$\vartheta=0$. Then one has:

(29)
${\boldsymbol{W}}_{\varepsilon}[\boldsymbol{L}({\boldsymbol{s}},\varepsilon{\bf
D}){\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}]=\boldsymbol{L}\left({\boldsymbol{s}},{\boldsymbol{\xi}}\right){\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}]-\frac{\varepsilon}{\mathrm{i}}{\boldsymbol{\nabla}}_{\boldsymbol{s}}\boldsymbol{L}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{\xi}}{\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}]\\\
-\frac{\varepsilon}{\mathrm{i}}\left({\boldsymbol{\nabla}}_{\boldsymbol{s}}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{\xi}}\boldsymbol{L}\right){\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}]+\mathrm{O}(\varepsilon^{2})\,,$

and:

(30)
$\begin{split}{\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},\boldsymbol{L}({\boldsymbol{s}},\varepsilon{\bf
D}){\boldsymbol{v}}_{\varepsilon}]&={(\boldsymbol{L}({\boldsymbol{s}},{\boldsymbol{\xi}}-\varepsilon{\bf
D}){\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}])}^{*}\\\
&={\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}]{\boldsymbol{L}}^{*}({\boldsymbol{s}},{\boldsymbol{\xi}})-\frac{\varepsilon}{\mathrm{i}}{\boldsymbol{\nabla}}_{\boldsymbol{s}}{\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}]\cdot{\boldsymbol{\nabla}}_{\boldsymbol{\xi}}{\boldsymbol{L}}^{*}+\mathrm{O}(\varepsilon^{2})\,.\end{split}$

Here
${\boldsymbol{\nabla}}_{\boldsymbol{s}}{\boldsymbol{A}}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{\xi}}{\boldsymbol{B}}:={\boldsymbol{\nabla}}_{{\boldsymbol{x}}}{\boldsymbol{A}}\cdot{\boldsymbol{\nabla}}_{{\boldsymbol{k}}}{\boldsymbol{B}}+\partial_{t}{\boldsymbol{A}}\partial_{\omega}{\boldsymbol{B}}$,
and the differential operator $\smash{\varepsilon{\bf D}}$ within the
observable $\boldsymbol{L}$ acts on
$\smash{{\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}]}$
so that:

$\boldsymbol{L}({\boldsymbol{s}},{\boldsymbol{\xi}}-\varepsilon{\bf
D}){\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}]({\boldsymbol{s}},{\boldsymbol{\xi}})=\frac{1}{(2\pi)^{4}}\int_{\mathbb{R}^{4}}\operatorname{e}^{\mathrm{i}{\boldsymbol{\eta}}\cdot{\boldsymbol{s}}}\boldsymbol{L}({\boldsymbol{s}},{\boldsymbol{\xi}}-\varepsilon{\boldsymbol{\eta}})\widehat{{\boldsymbol{W}}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}]({\boldsymbol{\eta}},{\boldsymbol{\xi}}){\mathrm{d}}{\boldsymbol{\eta}}\,,$

where
$\smash{\widehat{{\boldsymbol{W}}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}]({\boldsymbol{\eta}},{\boldsymbol{\xi}})}$
is the Fourier transform (25) of
$\smash{{\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}]({\boldsymbol{s}},{\boldsymbol{\xi}})}$
with respect to the space-time variable ${\boldsymbol{s}}$.

###  3.2. Maxwell’s equations as a semiclassical operator

Here Eq. (10) is conveniently written as a pseudo-differential operator for
the derivation of the high-frequency regime $\varepsilon\to 0$. We consider
slowly varying materials characterized by the Hermitian optical response
${\boldsymbol{K}}_{0}({\boldsymbol{x}})$ which is independent of the small
parameter $\varepsilon$ and we ignore its random fluctuations for the time
being ($\sigma=0$). Premultiplying it by
$\smash{\frac{\varepsilon}{\mathrm{i}}}$, Maxwell’s system (4) reads:

(31) ${\boldsymbol{L}}_{\varepsilon}({\boldsymbol{s}},\varepsilon{\bf
D}_{\boldsymbol{s}}){\boldsymbol{u}}_{\varepsilon}={\bf
0}\,,\quad{\boldsymbol{s}}=({\boldsymbol{x}},t)\in{\mathcal{O}}\times\mathbb{R}_{+}^{*}\,,$

where $\smash{{\bf D}_{\boldsymbol{s}}}=\smash{({\bf
D}_{\boldsymbol{x}},\mathrm{D}_{t})}$ for $\smash{{\bf
D}_{\boldsymbol{x}}=\frac{1}{\mathrm{i}}{\boldsymbol{\nabla}}_{\boldsymbol{x}}}$,
$\smash{\mathrm{D}_{t}=\frac{1}{\mathrm{i}}\partial_{t}}$. The pseudo-
differential operator ${\boldsymbol{L}}_{\varepsilon}$ is given by:

(32)
${\boldsymbol{L}}_{\varepsilon}({\boldsymbol{s}},{\boldsymbol{\xi}})=\omega{\boldsymbol{I}}+{\boldsymbol{L}}_{0}({\boldsymbol{x}},{\boldsymbol{k}})+\frac{\varepsilon}{\mathrm{i}}{\boldsymbol{L}}_{1}({\boldsymbol{x}},\omega)\,,\\\
$

where the $6\times 6$ matrices
${\boldsymbol{L}}_{0}({\boldsymbol{x}},{\boldsymbol{k}})$ and
${\boldsymbol{L}}_{1}({\boldsymbol{x}},\omega)$ are given by Eq. (13), and the
symmetric Maxwell operator $\boldsymbol{M}({\boldsymbol{k}})$ is given by Eq.
(11). In (32) it is assumed that the rescaled susceptibility kernel
$\smash{{\boldsymbol{K}}_{d}}$ acts as:

$({\boldsymbol{K}}_{d}\star{\boldsymbol{u}}_{\varepsilon})({\boldsymbol{x}},t)=\int_{0}^{t}{\boldsymbol{K}}_{d}\left({\boldsymbol{x}},\frac{\tau}{\varepsilon}\right){\boldsymbol{u}}_{\varepsilon}({\boldsymbol{x}},t-\tau){\mathrm{d}}\tau\,.$

Also $\smash{{\boldsymbol{u}}_{\varepsilon}({\boldsymbol{x}},t)}$ stands for
the electromagnetic field satisfying Eq. (31) for highly oscillating initial
conditions characterized by the small scale $\varepsilon$. We derive its
Wigner measure and its evolution properties inside ${\mathcal{O}}$ in the
subsequent subsections, ignoring boundary effects for now. They can be
addressed as in [1] for example, at least for the hyperbolic-elliptic set.

######  Remark 1.

All derivations involving the Wigner transforms and Wigner measures of
$\smash{({\boldsymbol{u}}_{\varepsilon})}$ in the subsequent Sect. 3.3 and
Sect. 3.4 are formal. They can be made rigorous as follows. For all
${\boldsymbol{s}}\in{\mathcal{O}}\times\smash{\mathbb{R}^{*}_{+}}$, let
${\mathcal{U}}$ be an open set such that
$\overline{\mathcal{U}}\subset{\mathcal{O}}\times\smash{\mathbb{R}^{*}_{+}}$,
${\boldsymbol{s}}\in{\mathcal{U}}$, and let
$\varphi\in\smash{{\mathcal{C}}^{\infty}_{0}({\mathcal{O}}\times\mathbb{R}^{*}_{+})}$
with $\varphi\equiv 1$ on $\overline{\mathcal{U}}$. Multiplying Eq. (31) by
$\varphi$, we obtain an alternative equation where
$\smash{{\boldsymbol{u}}_{\varepsilon}}$ is replaced by
$\varphi\smash{{\boldsymbol{u}}_{\varepsilon}}$, with additional terms. The
Wigner measures related to these additional terms are concentrated outside
$\overline{\mathcal{U}}$ because of the semiclassical pseudo-differential
calculus [53, Theorem 4.24]. Besides, it can be shown that
${\boldsymbol{W}}[\varphi{\boldsymbol{u}}_{\varepsilon}]({\boldsymbol{s}},{\boldsymbol{\xi}})$
is actually independent of $\varphi$ choosen as at the beginning of this
remark. In the sequel we shall define the Wigner measure of
$\smash{({\boldsymbol{u}}_{\varepsilon})}$,
$\smash{{\boldsymbol{W}}[{\boldsymbol{u}}_{\varepsilon}]}({\boldsymbol{s}},{\boldsymbol{\xi}})$
for ${\boldsymbol{s}}\in{\mathcal{O}}\times\smash{\mathbb{R}^{*}_{+}}$,
${\boldsymbol{\xi}}\in\mathbb{R}^{4}\setminus\\{{\boldsymbol{k}}={\bf 0}\\}$,
as
$\smash{{\boldsymbol{W}}[\varphi{\boldsymbol{u}}_{\varepsilon}]}({\boldsymbol{s}},{\boldsymbol{\xi}})$,
with this choice of $\varphi$.

###  3.3. Dispersion properties

The foregoing pseudo-differential calculus and space-time Wigner transform of
Sect. 3.1 are now used with Maxwell’s system (31). Computing the space-time
Wigner transform (28) of ${\boldsymbol{u}}_{\varepsilon}$ from Eq. (31),
yields:

(33)
${\boldsymbol{W}}_{\varepsilon}\left[{\boldsymbol{L}}_{\varepsilon}({\boldsymbol{s}},\varepsilon{\bf
D}_{{\boldsymbol{s}}}){\boldsymbol{u}}_{\varepsilon},{\boldsymbol{u}}_{\varepsilon}\right]-{\boldsymbol{W}}_{\varepsilon}\left[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{L}}_{\varepsilon}({\boldsymbol{s}},\varepsilon{\bf
D}_{{\boldsymbol{s}}}){\boldsymbol{u}}_{\varepsilon}\right]={\bf 0}\,.$

However, invoking the rules (29) and (30) we get:

(34)
${\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{L}}_{\varepsilon}({\boldsymbol{s}},\varepsilon{\bf
D}_{{\boldsymbol{s}}}){\boldsymbol{u}}_{\varepsilon},{\boldsymbol{u}}_{\varepsilon}]=\left(\omega{\boldsymbol{I}}+{\boldsymbol{L}}_{0}+\frac{\varepsilon}{\mathrm{i}}{\boldsymbol{L}}_{1}\right){\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]-\frac{\varepsilon}{\mathrm{i}}{\boldsymbol{\nabla}}_{\boldsymbol{x}}{\boldsymbol{L}}_{0}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}{\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]\\\
-\frac{\varepsilon}{\mathrm{i}}({\boldsymbol{\nabla}}_{\boldsymbol{x}}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}{\boldsymbol{L}}_{0}){\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]+\mathrm{O}(\varepsilon^{2})\,,$

and:

(35)
${\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{L}}_{\varepsilon}({\boldsymbol{s}},\varepsilon{\bf
D}_{{\boldsymbol{s}}}){\boldsymbol{u}}_{\varepsilon}]={\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]\left(\omega{\boldsymbol{I}}+{{\boldsymbol{L}}}^{*}_{0}-\frac{\varepsilon}{\mathrm{i}}{{\boldsymbol{L}}}^{*}_{1}\right)-\frac{\varepsilon}{\mathrm{i}}\partial_{t}{\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]\\\
-\frac{\varepsilon}{\mathrm{i}}{\boldsymbol{\nabla}}_{\boldsymbol{x}}{\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}{{\boldsymbol{L}}}^{*}_{0}+\mathrm{O}(\varepsilon^{2})\,.$

Considering the leading-order term, one obtains as $\varepsilon\rightarrow 0$
(see also for instance [53, Theorem 5.3]):

(36)
${\boldsymbol{L}}_{0}({\boldsymbol{x}},{\boldsymbol{k}}){\boldsymbol{W}}({\boldsymbol{s}},{\boldsymbol{\xi}})-{\boldsymbol{W}}({\boldsymbol{s}},{\boldsymbol{\xi}}){{\boldsymbol{L}}}^{*}_{0}({\boldsymbol{x}},{\boldsymbol{k}})={\bf
0}$

for the Wigner measure
${\boldsymbol{W}}:={\boldsymbol{W}}[{\boldsymbol{u}}_{\varepsilon}]$
(independent of $\varepsilon$) of the sequence
$\smash{({\boldsymbol{u}}_{\varepsilon})}$ (dependent of $\varepsilon$) given
by Eq. (27) in phase space
$({\boldsymbol{s}},{\boldsymbol{\xi}})\in{\mathcal{O}}\times\smash{\mathbb{R}_{+}^{*}}\times\smash{\mathbb{R}^{4}_{\boldsymbol{\xi}}}\setminus\\{{\boldsymbol{k}}={\bf
0}\\}$. Again, the notation of [23, p. 330] is used and the left-hand side in
Eq. (36) above is clearly independent of $\varepsilon$. The Wigner measure is
a non-negative, Hermitian matrix. Besides, reminding Eq. (14) one deduces from
Eqs. (34)–(36):

(37)
${\boldsymbol{W}}=\sum_{\alpha}\delta(\omega+\omega_{\alpha}){\boldsymbol{W}}_{\alpha}\,,$

where
${\boldsymbol{W}}_{\alpha}=\boldsymbol{\Pi}_{\alpha}{\boldsymbol{W}}={\boldsymbol{W}}{\boldsymbol{\Pi}}^{*}_{\alpha}=\boldsymbol{\Pi}_{\alpha}{\boldsymbol{W}}{\boldsymbol{\Pi}}^{*}_{\alpha}$;
this is Eq. (16). We note that here we do not consider eigenvalues crossing,
when their multiplicities possibly change. We leave this case for future
works. It is considered for example in [20] for the Schrödinger equation.

######  Remark 2.

As noticed in Sect. 2.3, the non propagative mode $\smash{\omega_{0}}=0$ (for
$\alpha=\text{"0"}$) is always an eigenvalue of $\smash{{\boldsymbol{L}}_{0}}$
with multiplicity $2$, with the eigenvectors
$\smash{{\boldsymbol{a}}_{1}}=\smash{({\hat{\boldsymbol{k}}},{\bf 0})}$ and
$\smash{{\boldsymbol{a}}_{2}}=\smash{({\bf 0},{\hat{\boldsymbol{k}}})}$.
Recall that the optical response ${\boldsymbol{K}}_{0}$ is positive definite
and let $\boldsymbol{\mathcal{K}}$ be the $2\times 2$ matrix with elements
$\smash{{\mathcal{K}}_{jk}}=\smash{\left\langle{\boldsymbol{a}}_{j},{\boldsymbol{a}}_{k}\right\rangle}$,
$j,k=1,2$. One has
$\boldsymbol{\mathcal{K}}={\boldsymbol{Q}}^{*}{\boldsymbol{Q}}$ where:

${\boldsymbol{Q}}=\begin{bmatrix}\sqrt{\hat{\epsilon}_{0}}&\frac{\hat{\xi}_{0}}{\sqrt{\hat{\epsilon}_{0}}}\\\
0&\sqrt{\hat{\mu}_{0}-\frac{|\hat{\xi}_{0}|^{2}}{\hat{\epsilon}_{0}}}\end{bmatrix}\,,$

with
$\smash{\hat{\epsilon}_{0}({\boldsymbol{x}},{\boldsymbol{k}})}=\smash{{{\hat{\boldsymbol{k}}}}^{*}{\boldsymbol{\epsilon}}_{0}({\boldsymbol{x}}){\hat{\boldsymbol{k}}}}$,
$\smash{\hat{\mu}_{0}({\boldsymbol{x}},{\boldsymbol{k}})}=\smash{{{\hat{\boldsymbol{k}}}}^{*}{\boldsymbol{\mu}}_{0}({\boldsymbol{x}}){\hat{\boldsymbol{k}}}}$,
and
$\smash{\hat{\xi}_{0}({\boldsymbol{x}},{\boldsymbol{k}})}=\smash{{{\hat{\boldsymbol{k}}}}^{*}{\boldsymbol{\xi}}_{0}({\boldsymbol{x}}){\hat{\boldsymbol{k}}}}$.
Then the orthonormal eigenvectors with respect to the scalar product
$\left\langle\cdot,\cdot\right\rangle$ are
$\smash{{\boldsymbol{b}}_{0}^{j}}={\boldsymbol{Q}}^{-1}{\boldsymbol{a}}_{j}$
since
$\smash{{{\boldsymbol{a}}}^{*}_{j}{\boldsymbol{a}}_{k}}=\smash{\delta_{jk}}$,
yielding:

(38)
${\boldsymbol{b}}_{0}^{(1)}=\frac{1}{\sqrt{\hat{\epsilon}_{0}}}\begin{pmatrix}{\hat{\boldsymbol{k}}}\\\
{\bf
0}\end{pmatrix}\,,\quad{\boldsymbol{b}}_{0}^{(2)}=\frac{1}{\sqrt{\hat{\epsilon}_{0}\hat{\mu}_{0}-|\hat{\xi}_{0}|^{2}}}\begin{pmatrix}-\frac{\hat{\xi}_{0}}{\sqrt{\hat{\epsilon}_{0}}}{\hat{\boldsymbol{k}}}\\\
\sqrt{\hat{\epsilon}_{0}}\,{\hat{\boldsymbol{k}}}\end{pmatrix}\,.$

###  3.4. Evolution properties

Gathering the first order terms proportional to $\varepsilon$ in Eq. (33) one
obtains in the limit $\varepsilon\to 0$:

(39)
$\partial_{t}{\boldsymbol{W}}+{\boldsymbol{\nabla}}_{\boldsymbol{x}}{\boldsymbol{W}}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}{{\boldsymbol{L}}}^{*}_{0}-{\boldsymbol{\nabla}}_{\boldsymbol{x}}{\boldsymbol{L}}_{0}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}{\boldsymbol{W}}-({\boldsymbol{\nabla}}_{\boldsymbol{x}}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}{\boldsymbol{L}}_{0}){\boldsymbol{W}}+{\boldsymbol{L}}_{1}{\boldsymbol{W}}+{\boldsymbol{W}}{{\boldsymbol{L}}}^{*}_{1}={\bf
0}\,.$

Next we multiply Eq. (39) on the left by $\boldsymbol{\Pi}_{\alpha}$ and on
the right by ${\boldsymbol{\Pi}}^{*}_{\alpha}$. From the property of
projectors
$\boldsymbol{\Pi}_{\alpha}\boldsymbol{\Pi}_{\beta}=\delta_{\alpha\beta}\boldsymbol{\Pi}_{\alpha}$
we first derive:

$\begin{split}\boldsymbol{\Pi}_{\alpha}\partial_{x_{j}}{\boldsymbol{L}}_{0}&=\boldsymbol{\Pi}_{\alpha}\sum_{\beta}\left((\partial_{x_{j}}\omega_{\beta})\boldsymbol{\Pi}_{\beta}+\omega_{\beta}(\partial_{x_{j}}\boldsymbol{\Pi}_{\beta})\right)\\\
&=(\partial_{x_{j}}\omega_{\alpha})\boldsymbol{\Pi}_{\alpha}+(\partial_{x_{j}}\boldsymbol{\Pi}_{\alpha})(\omega_{\alpha}-{\boldsymbol{L}}_{0})\,.\end{split}$

The same holds for partial derivatives with respect to $k_{j}$. Besides:

$\begin{split}(\partial_{k_{j}}{\boldsymbol{W}}){\boldsymbol{\Pi}}^{*}_{\alpha}&=\partial_{k_{j}}({\boldsymbol{W}}{\boldsymbol{\Pi}}^{*}_{\alpha})-{\boldsymbol{W}}(\partial_{k_{j}}{\boldsymbol{\Pi}}^{*}_{\alpha})\\\
&=\partial_{k_{j}}{\boldsymbol{W}}_{\alpha}-{\boldsymbol{W}}{(\partial_{k_{j}}\boldsymbol{\Pi}_{\alpha})}^{*}\,\end{split}$

where again the same holds for partial derivatives with respect to $x_{j}$.
Therefore (with summation over the repeated index $j$):

$\begin{split}\boldsymbol{\Pi}_{\alpha}(\partial_{x_{j}}{\boldsymbol{L}}_{0})(\partial_{k_{j}}{\boldsymbol{W}}){\boldsymbol{\Pi}}^{*}_{\alpha}&=(\partial_{x_{j}}\omega_{\alpha})(\partial_{k_{j}}(\boldsymbol{\Pi}_{\alpha}{\boldsymbol{W}}_{\alpha})-(\partial_{k_{j}}\boldsymbol{\Pi}_{\alpha}){\boldsymbol{W}}_{\alpha}-{\boldsymbol{W}}_{\alpha}{(\partial_{k_{j}}\boldsymbol{\Pi}_{\alpha})}^{*})\\\
&\quad+(\partial_{x_{j}}\boldsymbol{\Pi}_{\alpha})(\omega_{\alpha}-{\boldsymbol{L}}_{0})(\partial_{k_{j}}{\boldsymbol{W}}_{\alpha}-{\boldsymbol{W}}{(\partial_{k_{j}}\boldsymbol{\Pi}_{\alpha})}^{*})\\\
&=(\partial_{x_{j}}\omega_{\alpha})(\partial_{k_{j}}{\boldsymbol{W}}_{\alpha}-(\partial_{k_{j}}\boldsymbol{\Pi}_{\alpha}){\boldsymbol{W}}_{\alpha}-{\boldsymbol{W}}_{\alpha}{(\partial_{k_{j}}\boldsymbol{\Pi}_{\alpha})}^{*})\\\
&\quad-(\partial_{x_{j}}\boldsymbol{\Pi}_{\alpha})((\partial_{k_{j}}(\omega_{\alpha}-{\boldsymbol{L}}_{0})){\boldsymbol{W}}_{\alpha}+(\omega_{\alpha}-{\boldsymbol{L}}_{0}){\boldsymbol{W}}{(\partial_{k_{j}}\boldsymbol{\Pi}_{\alpha})}^{*})\end{split}$

since $(\omega_{\alpha}-{\boldsymbol{L}}_{0}){\boldsymbol{W}}_{\alpha}={\bf
0}$. Because $\partial_{k_{j}}$ and $\partial_{x_{j}}$ play symmetric roles,
one also obtains:

$\begin{split}\boldsymbol{\Pi}_{\alpha}(\partial_{x_{j}}{\boldsymbol{W}})(\partial_{k_{j}}{{\boldsymbol{L}}}^{*}_{0}){\boldsymbol{\Pi}}^{*}_{\alpha}&={(\boldsymbol{\Pi}_{\alpha}(\partial_{k_{j}}{\boldsymbol{L}}_{0})(\partial_{x_{j}}{\boldsymbol{W}}){\boldsymbol{\Pi}}^{*}_{\alpha})}^{*}\\\
&=(\partial_{k_{j}}\omega_{\alpha})(\partial_{x_{j}}{\boldsymbol{W}}_{\alpha}-(\partial_{x_{j}}\boldsymbol{\Pi}_{\alpha}){\boldsymbol{W}}_{\alpha}-{\boldsymbol{W}}_{\alpha}{(\partial_{x_{j}}\boldsymbol{\Pi}_{\alpha})}^{*})\\\
&\quad-({\boldsymbol{W}}_{\alpha}\partial_{x_{j}}(\omega_{\alpha}-{{\boldsymbol{L}}}^{*}_{0})+(\partial_{x_{j}}\boldsymbol{\Pi}_{\alpha}){\boldsymbol{W}}(\omega_{\alpha}-{{\boldsymbol{L}}}^{*}_{0})){(\partial_{k_{j}}\boldsymbol{\Pi}_{\alpha})}^{*}\,.\end{split}$

Gathering these expressions, Eq. (39) yields:

(40)
$\partial_{t}{\boldsymbol{W}}_{\alpha}+\\{\omega_{\alpha},{\boldsymbol{W}}_{\alpha}\\}+({\boldsymbol{\nabla}}_{\boldsymbol{x}}\omega_{\alpha}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}\boldsymbol{\Pi}_{\alpha}-{\boldsymbol{\nabla}}_{\boldsymbol{x}}\boldsymbol{\Pi}_{\alpha}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}{\boldsymbol{L}}_{0}-\boldsymbol{\Pi}_{\alpha}{\boldsymbol{\nabla}}_{\boldsymbol{x}}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}{\boldsymbol{L}}_{0}){\boldsymbol{W}}_{\alpha}\\\
-{\boldsymbol{W}}_{\alpha}{({\boldsymbol{\nabla}}_{\boldsymbol{k}}\omega_{\alpha}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{x}}\boldsymbol{\Pi}_{\alpha}-{\boldsymbol{\nabla}}_{\boldsymbol{k}}\boldsymbol{\Pi}_{\alpha}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{x}}{\boldsymbol{L}}_{0})}^{*}+{\boldsymbol{L}}_{\alpha}{\boldsymbol{W}}_{\alpha}+{\boldsymbol{W}}_{\alpha}{{\boldsymbol{L}}}^{*}_{\alpha}={\bf
0}\,,$

where
$\\{\omega_{\alpha},{\boldsymbol{W}}_{\alpha}\\}={\boldsymbol{\nabla}}_{\boldsymbol{k}}\omega_{\alpha}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{x}}{\boldsymbol{W}}_{\alpha}-{\boldsymbol{\nabla}}_{\boldsymbol{x}}\omega_{\alpha}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}{\boldsymbol{W}}_{\alpha}$
and
$\smash{{\boldsymbol{L}}_{\alpha}}=\smash{\boldsymbol{\Pi}_{\alpha}{\boldsymbol{L}}_{1}}$.

Next consider the matrix
${\boldsymbol{M}}_{\alpha}=\partial_{x_{j}}\omega_{\alpha}\partial_{k_{j}}\boldsymbol{\Pi}_{\alpha}-\partial_{x_{j}}\boldsymbol{\Pi}_{\alpha}\partial_{k_{j}}{\boldsymbol{L}}_{0}-\boldsymbol{\Pi}_{\alpha}(\partial_{x_{j}}\partial_{k_{j}}{\boldsymbol{L}}_{0})$;
then:

$\begin{split}{\boldsymbol{M}}_{\alpha}&=\partial_{k_{j}}[(\partial_{x_{j}}\omega_{\alpha})\boldsymbol{\Pi}_{\alpha}]-(\partial_{x_{j}}\partial_{k_{j}}\omega_{\alpha})\boldsymbol{\Pi}_{\alpha}-\partial_{x_{j}}(\boldsymbol{\Pi}_{\alpha}\partial_{k_{j}}{\boldsymbol{L}}_{0})\\\
&=\partial_{k_{j}}[\cancel{\partial_{x_{j}}(\omega_{\alpha}\boldsymbol{\Pi}_{\alpha})}-\omega_{\alpha}\partial_{x_{j}}\boldsymbol{\Pi}_{\alpha}]-(\partial_{x_{j}}\partial_{k_{j}}\omega_{\alpha})\boldsymbol{\Pi}_{\alpha}-\partial_{x_{j}}[\cancel{\partial_{k_{j}}(\boldsymbol{\Pi}_{\alpha}{\boldsymbol{L}}_{0})}-(\partial_{k_{j}}\boldsymbol{\Pi}_{\alpha}){\boldsymbol{L}}_{0}]\\\
&=-(\partial_{k_{j}}\omega_{\alpha})(\partial_{x_{j}}\boldsymbol{\Pi}_{\alpha})-(\partial_{k_{j}}\partial_{x_{j}}\boldsymbol{\Pi}_{\alpha})(\omega_{\alpha}-{\boldsymbol{L}}_{0})-(\partial_{x_{j}}\partial_{k_{j}}\omega_{\alpha})\boldsymbol{\Pi}_{\alpha}+(\partial_{k_{j}}\boldsymbol{\Pi}_{\alpha})(\partial_{x_{j}}{\boldsymbol{L}}_{0})\,.\end{split}$

Letting (because again
$(\omega_{\alpha}-{\boldsymbol{L}}_{0}){\boldsymbol{W}}_{\alpha}={\bf 0}$):

(41)
${\boldsymbol{N}}_{\alpha}={\boldsymbol{\nabla}}_{\boldsymbol{k}}\boldsymbol{\Pi}_{\alpha}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{x}}{\boldsymbol{L}}_{0}-{\boldsymbol{\nabla}}_{\boldsymbol{k}}\omega_{\alpha}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{x}}\boldsymbol{\Pi}_{\alpha}-\frac{1}{2}({\boldsymbol{\nabla}}_{\boldsymbol{x}}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}\omega_{\alpha}){\boldsymbol{I}}\,,$

Eq. (40) finally reads:

(42)
$\partial_{t}{\boldsymbol{W}}_{\alpha}+\\{\omega_{\alpha},{\boldsymbol{W}}_{\alpha}\\}+({\boldsymbol{L}}_{\alpha}+{\boldsymbol{N}}_{\alpha}){\boldsymbol{W}}_{\alpha}+{\boldsymbol{W}}_{\alpha}{({\boldsymbol{L}}_{\alpha}+{\boldsymbol{N}}_{\alpha})}^{*}={\bf
0}\,.$

Alternatively, one may wish to expand the Wigner measure ${\boldsymbol{W}}$ of
Eq. (16) onto the eigenvectors ${\boldsymbol{b}}_{\alpha}$ as:

${\boldsymbol{W}}=\sum_{\alpha}\delta(\omega+\omega_{\alpha}){\boldsymbol{W}}_{\alpha}=\sum_{\alpha}\delta(\omega+\omega_{\alpha}){\boldsymbol{b}}_{\alpha}{\bf
w}_{\alpha}{{\boldsymbol{b}}}^{*}_{\alpha}\,,$

where ${\bf
w}_{\alpha}={{\boldsymbol{c}}}^{*}_{\alpha}{\boldsymbol{W}}_{\alpha}{\boldsymbol{c}}_{\alpha}$.
The evolution properties of the so-called specific intensity, or coherence
matrices ${\bf w}_{\alpha}$ are derived multiplying Eq. (42) by
${{\boldsymbol{c}}}^{*}_{\alpha}$ on the left-hand side, and by
${\boldsymbol{c}}_{\alpha}$ on the right-hand side. Firstly:

$\begin{split}{{\boldsymbol{c}}}^{*}_{\alpha}(\partial_{k_{j}}\omega_{\alpha})(\partial_{x_{j}}{\boldsymbol{W}}_{\alpha}){\boldsymbol{c}}_{\alpha}&=\partial_{k_{j}}\omega_{\alpha}[\partial_{x_{j}}{\bf
w}_{\alpha}-(\partial_{x_{j}}{{\boldsymbol{c}}}^{*}_{\alpha}){\boldsymbol{W}}_{\alpha}{\boldsymbol{c}}_{\alpha}-{{\boldsymbol{c}}}^{*}_{\alpha}{\boldsymbol{W}}_{\alpha}\partial_{x_{j}}{\boldsymbol{c}}_{\alpha}]\\\
&=\partial_{k_{j}}\omega_{\alpha}[\partial_{x_{j}}{\bf
w}_{\alpha}-(\partial_{x_{j}}{{\boldsymbol{c}}}^{*}_{\alpha}){\boldsymbol{b}}_{\alpha}{\bf
w}_{\alpha}-{\bf
w}_{\alpha}{{\boldsymbol{b}}}^{*}_{\alpha}\partial_{x_{j}}{\boldsymbol{c}}_{\alpha}]\\\
&=\partial_{k_{j}}\omega_{\alpha}[\partial_{x_{j}}{\bf
w}_{\alpha}+{{\boldsymbol{c}}}^{*}_{\alpha}(\partial_{x_{j}}{\boldsymbol{b}}_{\alpha}){\bf
w}_{\alpha}+{\bf
w}_{\alpha}(\partial_{x_{j}}{{\boldsymbol{b}}}^{*}_{\alpha}){\boldsymbol{c}}_{\alpha}]\end{split}$

owing to the orthonormality condition
${{\boldsymbol{c}}}^{*}_{\alpha}{\boldsymbol{b}}_{\beta}=\delta_{\alpha\beta}{\boldsymbol{I}}_{A}$
which yields
$(\partial_{x_{j}}{{\boldsymbol{c}}}^{*}_{\alpha}){\boldsymbol{b}}_{\beta}=-{{\boldsymbol{c}}}^{*}_{\alpha}(\partial_{x_{j}}{\boldsymbol{b}}_{\beta})$
(and a similar property with $\partial_{k_{j}}$). Likewise:

${{\boldsymbol{c}}}^{*}_{\alpha}(\partial_{x_{j}}\omega_{\alpha})(\partial_{k_{j}}{\boldsymbol{W}}_{\alpha}){\boldsymbol{c}}_{\alpha}=\partial_{x_{j}}\omega_{\alpha}[\partial_{k_{j}}{\bf
w}_{\alpha}+{{\boldsymbol{c}}}^{*}_{\alpha}(\partial_{k_{j}}{\boldsymbol{b}}_{\alpha}){\bf
w}_{\alpha}+{\bf
w}_{\alpha}(\partial_{k_{j}}{{\boldsymbol{b}}}^{*}_{\alpha}){\boldsymbol{c}}_{\alpha}]\,,$

and consequently:

${{\boldsymbol{c}}}^{*}_{\alpha}\\{\omega_{\alpha},{\boldsymbol{W}}_{\alpha}\\}{\boldsymbol{c}}_{\alpha}=\\{\omega_{\alpha},{\bf
w}_{\alpha}\\}+{{\boldsymbol{c}}}^{*}_{\alpha}\\{\omega_{\alpha},{\boldsymbol{b}}_{\alpha}\\}{\bf
w}_{\alpha}+{\bf
w}_{\alpha}\\{\omega_{\alpha},{{\boldsymbol{b}}}^{*}_{\alpha}\\}{\boldsymbol{c}}_{\alpha}\,.$

Secondly:

$\begin{split}{{\boldsymbol{c}}}^{*}_{\alpha}(\partial_{k_{j}}\boldsymbol{\Pi}_{\alpha})(\partial_{x_{j}}{\boldsymbol{L}}_{0}){\boldsymbol{W}}_{\alpha}{\boldsymbol{c}}_{\alpha}&=[\partial_{k_{j}}({{\boldsymbol{c}}}^{*}_{\alpha}\boldsymbol{\Pi}_{\alpha})-(\partial_{k_{j}}{{\boldsymbol{c}}}^{*}_{\alpha})\boldsymbol{\Pi}_{\alpha}](\partial_{x_{j}}{\boldsymbol{L}}_{0}){\boldsymbol{b}}_{\alpha}{\bf
w}_{\alpha}\\\
&=(\partial_{k_{j}}{{\boldsymbol{c}}}^{*}_{\alpha})({\boldsymbol{I}}-\boldsymbol{\Pi}_{\alpha})[\partial_{x_{j}}({\boldsymbol{L}}_{0}{\boldsymbol{b}}_{\alpha})-{\boldsymbol{L}}_{0}\partial_{x_{j}}{\boldsymbol{b}}_{\alpha}]{\bf
w}_{\alpha}\\\
&=(\partial_{k_{j}}{{\boldsymbol{c}}}^{*}_{\alpha})({\boldsymbol{I}}-\boldsymbol{\Pi}_{\alpha})[\partial_{x_{j}}(\omega_{\alpha}{\boldsymbol{b}}_{\alpha})-{\boldsymbol{L}}_{0}\partial_{x_{j}}{\boldsymbol{b}}_{\alpha}]{\bf
w}_{\alpha}\\\
&=(\partial_{k_{j}}{{\boldsymbol{c}}}^{*}_{\alpha})({\boldsymbol{I}}-\boldsymbol{\Pi}_{\alpha})(\omega_{\alpha}-{\boldsymbol{L}}_{0})(\partial_{x_{j}}{\boldsymbol{b}}_{\alpha}){\bf
w}_{\alpha}\\\
&=(\partial_{k_{j}}{{\boldsymbol{c}}}^{*}_{\alpha})(\omega_{\alpha}-{\boldsymbol{L}}_{0})(\partial_{x_{j}}{\boldsymbol{b}}_{\alpha}){\bf
w}_{\alpha}\\\
&=-{{\boldsymbol{c}}}^{*}_{\alpha}[\partial_{k_{j}}(\omega_{\alpha}-{\boldsymbol{L}}_{0})](\partial_{x_{j}}{\boldsymbol{b}}_{\alpha}){\bf
w}_{\alpha}\,;\end{split}$

and:

$\begin{split}{{\boldsymbol{c}}}^{*}_{\alpha}(\partial_{k_{j}}\omega_{\alpha})(\partial_{x_{j}}\boldsymbol{\Pi}_{\alpha}){\boldsymbol{W}}_{\alpha}{\boldsymbol{c}}_{\alpha}&=(\partial_{k_{j}}\omega_{\alpha})[\partial_{x_{j}}({{\boldsymbol{c}}}^{*}_{\alpha}\boldsymbol{\Pi}_{\alpha})-(\partial_{x_{j}}{{\boldsymbol{c}}}^{*}_{\alpha})\boldsymbol{\Pi}_{\alpha}]{\boldsymbol{b}}_{\alpha}{\bf
w}_{\alpha}\\\
&=(\partial_{k_{j}}\omega_{\alpha})(\partial_{x_{j}}{{\boldsymbol{c}}}^{*}_{\alpha})({\boldsymbol{I}}-\boldsymbol{\Pi}_{\alpha}){\boldsymbol{b}}_{\alpha}{\bf
w}_{\alpha}\\\ &={\bf 0}\,.\end{split}$

Thirdly:

$\begin{split}{{\boldsymbol{c}}}^{*}_{\alpha}\\{\omega_{\alpha},{\boldsymbol{b}}_{\alpha}\\}-{{\boldsymbol{c}}}^{*}_{\alpha}[\partial_{k_{j}}(\omega_{\alpha}-{\boldsymbol{L}}_{0})]\partial_{x_{j}}{\boldsymbol{b}}_{\alpha}&={{\boldsymbol{c}}}^{*}_{\alpha}[(\partial_{k_{j}}{\boldsymbol{L}}_{0})(\partial_{x_{j}}{\boldsymbol{b}}_{\alpha})-(\partial_{x_{j}}\omega_{\alpha})(\partial_{k_{j}}{\boldsymbol{b}}_{\alpha})]\,.\end{split}$

Letting:

(43)
$\begin{split}{\boldsymbol{n}}_{\alpha}&={{\boldsymbol{c}}}^{*}_{\alpha}(\\{\omega_{\alpha},{\boldsymbol{b}}_{\alpha}\\}+{\boldsymbol{N}}_{\alpha}{\boldsymbol{b}}_{\alpha})\\\
&={{\boldsymbol{c}}}^{*}_{\alpha}({\boldsymbol{\nabla}}_{\boldsymbol{k}}{\boldsymbol{L}}_{0}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{x}}{\boldsymbol{b}}_{\alpha}-{\boldsymbol{\nabla}}_{\boldsymbol{x}}\omega_{\alpha}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}{\boldsymbol{b}}_{\alpha})-\frac{1}{2}({\boldsymbol{\nabla}}_{\boldsymbol{x}}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}\omega_{\alpha}){\boldsymbol{I}}_{A}\,,\end{split}$

one can show that ${\boldsymbol{n}}_{\alpha}$ is skew-symmetric owing to Eq.
(13). Indeed:

$\begin{split}{{\boldsymbol{c}}}^{*}_{\alpha}(\partial_{k_{j}}{\boldsymbol{L}}_{0}){\boldsymbol{b}}_{\alpha}&={{\boldsymbol{c}}}^{*}_{\alpha}{\boldsymbol{K}}_{0}^{-1}({\boldsymbol{x}})(\partial_{k_{j}}\boldsymbol{M}({\boldsymbol{k}})){\boldsymbol{b}}_{\alpha}\\\
&={{\boldsymbol{b}}}^{*}_{\alpha}(\partial_{k_{j}}\boldsymbol{M}){\boldsymbol{b}}_{\alpha}\\\
&=(\partial_{k_{j}}\omega_{\alpha}){\boldsymbol{I}}_{A}\,.\end{split}$

Then
$(\partial_{x_{j}}{{\boldsymbol{b}}}^{*}_{\alpha})(\partial_{k_{j}}\boldsymbol{M}){\boldsymbol{b}}_{\alpha}+{{\boldsymbol{b}}}^{*}_{\alpha}(\partial_{k_{j}}\boldsymbol{M})(\partial_{x_{j}}{\boldsymbol{b}}_{\alpha})=(\partial_{x_{j}}\partial_{k_{j}}\omega_{\alpha}){\boldsymbol{I}}_{A}$,
or:

(44)
$\displaystyle{{\boldsymbol{c}}}^{*}_{\alpha}(\partial_{k_{j}}{\boldsymbol{L}}_{0})(\partial_{x_{j}}{\boldsymbol{b}}_{\alpha})$
$\displaystyle=(\partial_{x_{j}}\partial_{k_{j}}\omega_{\alpha}){\boldsymbol{I}}_{A}-{({{\boldsymbol{b}}}^{*}_{\alpha}(\partial_{k_{j}}\boldsymbol{M})(\partial_{x_{j}}{\boldsymbol{b}}_{\alpha}))}^{*}$
$\displaystyle=(\partial_{x_{j}}\partial_{k_{j}}\omega_{\alpha}){\boldsymbol{I}}_{A}-{({{\boldsymbol{c}}}^{*}_{\alpha}(\partial_{k_{j}}{\boldsymbol{L}}_{0})(\partial_{x_{j}}{\boldsymbol{b}}_{\alpha}))}^{*}\,;$

besides:

(45)
$\displaystyle(\partial_{x_{j}}\omega_{\alpha}){{\boldsymbol{c}}}^{*}_{\alpha}\partial_{k_{j}}{\boldsymbol{b}}_{\alpha}$
$\displaystyle=-(\partial_{x_{j}}\omega_{\alpha})(\partial_{k_{j}}{{\boldsymbol{c}}}^{*}_{\alpha}){\boldsymbol{b}}_{\alpha}$
$\displaystyle=-(\partial_{x_{j}}\omega_{\alpha})(\partial_{k_{j}}{{\boldsymbol{b}}}^{*}_{\alpha}){\boldsymbol{c}}_{\alpha}\,;$

thus one concludes with (44) and (45) that
${\boldsymbol{n}}_{\alpha}=-{{\boldsymbol{n}}}^{*}_{\alpha}$. Therefore,
setting
$\smash{{\boldsymbol{\ell}}_{\alpha}}=\smash{{{\boldsymbol{c}}}^{*}_{\alpha}{\boldsymbol{L}}_{1}{\boldsymbol{b}}_{\alpha}}$,
Eq. (42) for the specific intensities ${\bf w}_{\alpha}$ finally reads:

(46) $\partial_{t}{\bf w}_{\alpha}+\\{\omega_{\alpha},{\bf
w}_{\alpha}\\}+{\boldsymbol{\ell}}_{\alpha}{\bf w}_{\alpha}+{\bf
w}_{\alpha}{{\boldsymbol{\ell}}}^{*}_{\alpha}+{\boldsymbol{n}}_{\alpha}{\bf
w}_{\alpha}-{\bf w}_{\alpha}{\boldsymbol{n}}_{\alpha}={\bf 0}\,.$

One can check in particular that for the eigenspace (38) associated with the
null eigenvalue ($\alpha=\text{"0"}$), one has
${\boldsymbol{n}}_{0}={{\boldsymbol{c}}}^{*}_{0}{\boldsymbol{\nabla}}_{\boldsymbol{k}}{\boldsymbol{L}}_{0}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{x}}{\boldsymbol{b}}_{0}={{\boldsymbol{b}}}^{*}_{0}{\boldsymbol{\nabla}}_{\boldsymbol{k}}\boldsymbol{M}({\boldsymbol{k}})\cdot{\boldsymbol{\nabla}}_{\boldsymbol{x}}{\boldsymbol{b}}_{0}={\bf
0}$.

The geometrical structure of transport equations (46) is analyzed in _e.g._
[34]. Here we show how solving Eq. (46) amounts to track $\smash{{\bf
w}_{\alpha}}$ on its bicharacteristic rays in
${\mathcal{O}}\times\mathbb{R}^{3}$. Consider the bicharacteristic
(Hamiltonian) equations associated to
$\omega_{\alpha}({\boldsymbol{x}},{\boldsymbol{k}})$:

(47)
$\begin{split}\frac{{\mathrm{d}}{\boldsymbol{x}}_{\alpha}}{{\mathrm{d}}t}&={\boldsymbol{\nabla}}_{\boldsymbol{k}}\omega_{\alpha}({\boldsymbol{x}}_{\alpha}(t),{\boldsymbol{k}}_{\alpha}(t))\,,\\\
\frac{{\mathrm{d}}{\boldsymbol{k}}_{\alpha}}{{\mathrm{d}}t}&=-{\boldsymbol{\nabla}}_{\boldsymbol{x}}\omega_{\alpha}({\boldsymbol{x}}_{\alpha}(t),{\boldsymbol{k}}_{\alpha}(t))\,,\end{split}$

with initial conditions
$\smash{{\boldsymbol{x}}_{\alpha}(0)}=\smash{{\boldsymbol{x}}_{0}}$ and
$\smash{{\boldsymbol{k}}_{\alpha}(0)}=\smash{{\boldsymbol{k}}_{0}}$ lying in
the support of the Wigner measure $\smash{{\boldsymbol{W}}_{I}}$ of the
initial conditions $\smash{{\boldsymbol{u}}_{\varepsilon}(\cdot,t=0)}$. Also
introduce
$\smash{{\boldsymbol{\Omega}}_{\alpha}}={\boldsymbol{\ell}}_{\alpha}+{\boldsymbol{n}}_{\alpha}$
and the $A\times A$ matrix $\smash{{\boldsymbol{R}}_{\alpha}}$ which is such
that:

(48)
$\frac{{\mathrm{d}}{\boldsymbol{R}}_{\alpha}}{{\mathrm{d}}t}=-{\boldsymbol{\Omega}}_{\alpha}({\boldsymbol{x}}_{\alpha}(t),{\boldsymbol{k}}_{\alpha}(t)){\boldsymbol{R}}_{\alpha}(t)\,,\quad{\boldsymbol{R}}_{\alpha}(0)={\boldsymbol{I}}_{A}\,.$

Then:

$\begin{split}{\boldsymbol{\Omega}}_{\alpha}{\bf w}_{\alpha}+{\bf
w}_{\alpha}{{\boldsymbol{\Omega}}}^{*}_{\alpha}&=-\frac{{\mathrm{d}}{\boldsymbol{R}}_{\alpha}}{{\mathrm{d}}t}{\boldsymbol{R}}^{-1}_{\alpha}{\bf
w}_{\alpha}-{\bf
w}_{\alpha}({{\boldsymbol{R}}}^{*}_{\alpha})^{-1}\frac{{\mathrm{d}}{{\boldsymbol{R}}}^{*}_{\alpha}}{{\mathrm{d}}t}\\\
&=-\frac{{\mathrm{d}}{\boldsymbol{R}}_{\alpha}}{{\mathrm{d}}t}\tilde{{\bf
w}}_{\alpha}{{\boldsymbol{R}}}^{*}_{\alpha}-{\boldsymbol{R}}_{\alpha}\tilde{{\bf
w}}_{\alpha}\frac{{\mathrm{d}}{{\boldsymbol{R}}}^{*}_{\alpha}}{{\mathrm{d}}t}\end{split}$

if one lets $\tilde{{\bf w}}_{\alpha}={\boldsymbol{R}}^{-1}_{\alpha}{\bf
w}_{\alpha}({{\boldsymbol{R}}}^{*}_{\alpha})^{-1}$, and owing to Eq. (47) the
transport equation (18) reads:

$\frac{{\mathrm{d}}{\bf
w}_{\alpha}}{{\mathrm{d}}t}-\frac{{\mathrm{d}}{\boldsymbol{R}}_{\alpha}}{{\mathrm{d}}t}\tilde{{\bf
w}}_{\alpha}{{\boldsymbol{R}}}^{*}_{\alpha}-{\boldsymbol{R}}_{\alpha}\tilde{{\bf
w}}_{\alpha}\frac{{\mathrm{d}}{{\boldsymbol{R}}}^{*}_{\alpha}}{{\mathrm{d}}t}={\bf
0}\,,$

or:

(49) $\frac{{\mathrm{d}}}{{\mathrm{d}}t}\tilde{{\bf
w}}_{\alpha}({\boldsymbol{x}}_{\alpha}(t),{\boldsymbol{k}}_{\alpha}(t),t)={\bf
0}\,.$

Consequently:

(50) ${\bf
w}_{\alpha}({\boldsymbol{x}}_{\alpha}(t),{\boldsymbol{k}}_{\alpha}(t),t)={\boldsymbol{R}}_{\alpha}(t){\bf
w}^{I}_{\alpha}({\boldsymbol{x}}_{0},{\boldsymbol{k}}_{0}){{\boldsymbol{R}}}^{*}_{\alpha}(t)$

where $\smash{{\bf
w}^{I}_{\alpha}}=\smash{{{\boldsymbol{c}}}^{*}_{\alpha}{\boldsymbol{W}}_{I}{\boldsymbol{c}}_{\alpha}}$,
and one comes down to solving Eqs. (47)–(48) in order to track $\smash{{\bf
w}_{\alpha}}$ on its bicharacteristics
$t\mapsto\smash{({\boldsymbol{x}}_{\alpha}(t),{\boldsymbol{k}}_{\alpha}(t))}$.

######  Example 2.

For the Lorentz model with damping (6) one has three eigenvalues
$\omega_{0}=0$ ($\alpha=\text{"0"}$),
$\omega_{+}({\boldsymbol{k}})=+c_{0}|{\boldsymbol{k}}|$ ($\alpha=\text{"+"}$),
and $\omega_{-}({\boldsymbol{k}})=-c_{0}|{\boldsymbol{k}}|$
($\alpha=\text{"$-$"}$), each of multiplicity two ($A=2$),
$\smash{c_{0}}=\smash{1/\sqrt{\epsilon_{0}\mu_{0}}}$, with associated
eigenvectors [46]:

${\boldsymbol{b}}_{0}^{(1)}({\boldsymbol{k}})=\frac{1}{\sqrt{\epsilon_{0}}}\begin{pmatrix}{\hat{\boldsymbol{k}}}\\\
{\bf
0}\end{pmatrix}\,,\quad{\boldsymbol{b}}_{0}^{(2)}({\boldsymbol{k}})=\frac{1}{\sqrt{\mu_{0}}}\begin{pmatrix}{\bf
0}\\\ {\hat{\boldsymbol{k}}}\end{pmatrix}\,,$

and for $j=1,2$:

${\boldsymbol{b}}_{\pm}^{(j)}({\boldsymbol{k}})=\begin{pmatrix}\displaystyle\frac{\hat{\boldsymbol{e}}_{j}}{\sqrt{2\epsilon_{0}}}\\\
\pm\displaystyle\frac{{\hat{\boldsymbol{k}}}\times\hat{\boldsymbol{e}}_{j}}{\sqrt{2\mu_{0}}}\end{pmatrix}\,,\quad\hat{\boldsymbol{e}}_{1}\perp\hat{\boldsymbol{e}}_{2}\in{\hat{\boldsymbol{k}}}^{\perp}\,,$

such that
$({\hat{\boldsymbol{k}}},\hat{\boldsymbol{e}}_{1}({\boldsymbol{k}}),\hat{\boldsymbol{e}}_{2}({\boldsymbol{k}}))$
form a right-handed orthonormal system. Besides
$\smash{{\boldsymbol{\ell}}_{0}}\smash{{\bf w}_{0}}={\bf 0}$ since
$\omega=-\smash{\omega_{0}}=0$ on the support of $\smash{{\bf w}_{0}}$, and:

${\boldsymbol{\ell}}_{\pm}=\frac{\mathrm{i}\omega\widehat{\epsilon}_{1}(\omega)}{2}\begin{bmatrix}1&0\\\
0&1\end{bmatrix}\,,$

such that for the propagating modes
$\smash{\omega_{\pm}}({\boldsymbol{k}})=\pm\smash{c_{0}}|{\boldsymbol{k}}|$:

(51) $\begin{split}{\boldsymbol{\ell}}_{\pm}{\bf w}_{\pm}+{\bf
w}_{\pm}{{\boldsymbol{\ell}}}^{*}_{\pm}&=\Re\mathfrak{e}\\{\mathrm{i}\omega\widehat{\epsilon}_{1}(\omega)\\}{\bf
w}_{\pm}\\\
&=\frac{\omega_{p}^{2}c_{0}^{2}|{\boldsymbol{k}}|^{2}\Gamma}{(\omega_{0}^{2}-c_{0}^{2}|{\boldsymbol{k}}|^{2})^{2}+c_{0}^{2}|{\boldsymbol{k}}|^{2}\Gamma^{2}}{\bf
w}_{\pm}\\\ &=\tilde{\Gamma}(c_{0}|{\boldsymbol{k}}|){\bf
w}_{\pm}\,.\end{split}$

Also $\smash{{\boldsymbol{n}}_{\pm}}={\bf 0}$ for an homogeneous material.
Therefore the solution of Eq. (46) reads:

(52) ${\bf
w}_{\pm}({\boldsymbol{x}},{\boldsymbol{k}},t)=\operatorname{e}^{-\tilde{\Gamma}(c_{0}\left|{\boldsymbol{k}}\right|)t}{\bf
w}_{I}({\boldsymbol{x}}\mp c_{0}{\hat{\boldsymbol{k}}}t,{\boldsymbol{k}})$

where ${\bf w}_{I}({\boldsymbol{x}},{\boldsymbol{k}})$,
$({\boldsymbol{x}},{\boldsymbol{k}})\in{\mathcal{O}}\times\mathbb{R}^{3}$, are
the initial (incident) specific intensities supported in some subset of the
physical domain ${\mathcal{O}}$. Note that similar results were obtained in
[5, 42] for memoryless viscoelastic media. Recently a semi-linear radiative
transfer model has been derived in media with non linear absorption [33].

##  4\. Radiative transfer

In this section the radiative transfer equations (22) for a bianistropic
dielectric medium are now derived. They describe the evolution of the high-
frequency electromagnetic energy density accounting for random fluctuations of
the optical response. We resort to the same formal pseudo-differential
calculus as recalled in Sect. 3.1. However additional rules are needed to
account for the dependency of the optical response on both the slow scale of
variation of the background medium, and the fast scale of variation of the
random fluctuations. In Sect. 4.1 below we start by rewriting Maxwell’s system
as a pseudo-differential operator when these random fluctuations are
considered. In Sect. 4.2 we give the additional rules of pseudo-differential
calculus for this situation, and we finally obtain in Sect. 4.3 the radiative
transfer equations by a multiscale analysis.

###  4.1. Maxwell’s equations with random perturbations

In view of the scattering regime of Eqs. (20)–(21), the rescaled optical
response ${\boldsymbol{K}}({\boldsymbol{x}})$ in Eq. (10) is written :

${\boldsymbol{K}}({\boldsymbol{x}})={\boldsymbol{K}}_{0}({\boldsymbol{x}})\left[{\boldsymbol{I}}+\sqrt{\varepsilon}{\boldsymbol{V}}\left(\frac{{\boldsymbol{x}}}{\varepsilon}\right)\right]\,.$

The power spectral density matrix $\widehat{{\boldsymbol{R}}}$ of the
homogeneous random field
$({\boldsymbol{V}}({\boldsymbol{y}}),\,{\boldsymbol{y}}\in\mathbb{R}^{3})$ is
given by Eq. (12) and satisfies:

(53)
${\mathbb{E}}\left\\{\smash{\widehat{{\boldsymbol{V}}}({\boldsymbol{p}})\otimes\widehat{{\boldsymbol{V}}}({\boldsymbol{q}})}\right\\}=(2\pi)^{6}\widehat{{\boldsymbol{R}}}({\boldsymbol{p}})\delta({\boldsymbol{p}}+{\boldsymbol{q}})\,,$

for the Fourier transform $\widehat{{\boldsymbol{V}}}$ of ${\boldsymbol{V}}$;
see Eq. (25). Besides, homogeneity (spatial stationarity) also implies that:

(54)
$\widehat{R}_{jklm}({\boldsymbol{p}})=\widehat{R}_{lmjk}(-{\boldsymbol{p}})\,,\quad
1\leq j,k,l,m\leq 6\,.$

We remind that in order to preserve the Hermiticity of the actual optical
response ${\boldsymbol{K}}$, it is assumed that
${\boldsymbol{K}}_{0}({\boldsymbol{x}}){\boldsymbol{V}}({\boldsymbol{y}})={{\boldsymbol{V}}}^{*}({\boldsymbol{y}}){\boldsymbol{K}}_{0}({\boldsymbol{x}})$
for all ${\boldsymbol{x}}\in{\mathcal{O}}$,
${\boldsymbol{y}}=\frac{{\boldsymbol{x}}}{\varepsilon}$. The pseudo-
differential operator ${\boldsymbol{L}}_{\varepsilon}$ of Eq. (32) is then
modified to:

(55)
$\mathcal{L}_{\varepsilon}\left({\boldsymbol{s}},\frac{{\boldsymbol{x}}}{\varepsilon},{\boldsymbol{\xi}}\right)={\boldsymbol{L}}_{\varepsilon}({\boldsymbol{s}},{\boldsymbol{\xi}})+\sqrt{\varepsilon}{\boldsymbol{L}}_{\text{0.5}}\left(\frac{{\boldsymbol{x}}}{\varepsilon},\omega\right)\,,$

where again ${\boldsymbol{s}}=({\boldsymbol{x}},t)$,
${\boldsymbol{\xi}}=({\boldsymbol{k}},\omega)$, and
$\smash{{\boldsymbol{L}}_{\text{0.5}}({\boldsymbol{y}},\omega)=\omega{\boldsymbol{V}}({\boldsymbol{y}})}$.
This additional term arises from the random fluctuations of the optical
response. Since it depends on the fast variable
${\boldsymbol{y}}=\frac{{\boldsymbol{x}}}{\varepsilon}$, the rules (29) and
(30) have to be adapted.

###  4.2. Rules of pseudo-differential calculus with oscillating coefficients

The rules (29) and (30) hold true if the Wigner transform
$\smash{{\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}]}$
does not depend on $\smash{\frac{{\boldsymbol{s}}}{\varepsilon}}$. If however
it depends on this fast variable they have to be modified accordingly. Let
${\boldsymbol{\tau}}\mapsto{\boldsymbol{A}}({\boldsymbol{\tau}})$ be a smooth
matrix-valued function. Then we have:

(56)
$\begin{split}{\boldsymbol{W}}_{\varepsilon}\left[{\boldsymbol{A}}\left(\frac{{\boldsymbol{s}}}{\varepsilon}\right){\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}\right]&=\frac{1}{(2\pi)^{4}}\int_{\mathbb{R}^{4}}\operatorname{e}^{\mathrm{i}\frac{{\boldsymbol{s}}}{\varepsilon}\cdot{\boldsymbol{\eta}}}\widehat{{\boldsymbol{A}}}({\boldsymbol{\eta}}){\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}]\left({\boldsymbol{s}},{\boldsymbol{\xi}}-{\boldsymbol{\eta}}\right){\mathrm{d}}{\boldsymbol{\eta}}\,,\\\
{\boldsymbol{W}}_{\varepsilon}\left[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{A}}\left(\frac{{\boldsymbol{s}}}{\varepsilon}\right){\boldsymbol{v}}_{\varepsilon}\right]&={\boldsymbol{W}}_{\varepsilon}\left[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}\right]{\boldsymbol{A}\left(\frac{{\boldsymbol{s}}}{\varepsilon}\right)}^{*}\,,\end{split}$

where $\widehat{{\boldsymbol{A}}}({\boldsymbol{\eta}})$ is the Fourier
transform (25) of ${\boldsymbol{A}}({\boldsymbol{\tau}})$. Applying these
formulas for a matrix-valued observable
$\smash{{\boldsymbol{A}}(\frac{{\boldsymbol{s}}}{\varepsilon})\boldsymbol{L}({\boldsymbol{s}},{\boldsymbol{\xi}})}$
yields:

(57)
${\boldsymbol{W}}_{\varepsilon}\left[{\boldsymbol{A}}\left(\frac{{\boldsymbol{s}}}{\varepsilon}\right)\boldsymbol{L}({\boldsymbol{s}},\varepsilon{\bf
D}){\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}\right]\\\
=\int\frac{\operatorname{e}^{\mathrm{i}\frac{{\boldsymbol{s}}}{\varepsilon}\cdot{\boldsymbol{\eta}}}{\mathrm{d}}{\boldsymbol{\eta}}}{(2\pi)^{4}}\widehat{{\boldsymbol{A}}}({\boldsymbol{\eta}})\boldsymbol{L}({\boldsymbol{s}},{\boldsymbol{\xi}}-{\boldsymbol{\eta}})W_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}]\left({\boldsymbol{s}},{\boldsymbol{\xi}}-{\boldsymbol{\eta}}\right)+\mathrm{O}(\varepsilon)\,,$

and:

(58)
${\boldsymbol{W}}_{\varepsilon}\left[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{A}}\left(\frac{{\boldsymbol{s}}}{\varepsilon}\right)\boldsymbol{L}({\boldsymbol{s}},\varepsilon{\bf
D}){\boldsymbol{v}}_{\varepsilon}\right]={\left({\boldsymbol{A}}\left(\frac{{\boldsymbol{s}}}{\varepsilon}\right)\boldsymbol{L}({\boldsymbol{s}},{\boldsymbol{\xi}}-\varepsilon{\bf
D}){\boldsymbol{W}}_{\varepsilon}\left[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{v}}_{\varepsilon}\right]\right)}^{*}\,,$

owing to Eq. (29) and Eq. (30). The proofs of these formulas are given in [2,
Appendix C] for scalar fields. Their extensions to vector fields
${\boldsymbol{u}}_{\varepsilon}$ and ${\boldsymbol{v}}_{\varepsilon}$ are
straightforward. Consequently one has:

(59)
$\begin{split}{\boldsymbol{W}}_{\varepsilon}\left[{\boldsymbol{L}}_{\text{0.5}}\left(\frac{{\boldsymbol{x}}}{\varepsilon},\varepsilon\mathrm{D}_{t}\right){\boldsymbol{u}}_{\varepsilon},{\boldsymbol{u}}_{\varepsilon}\right]&=\frac{\omega}{(2\pi)^{3}}\int_{\mathbb{R}^{3}}\operatorname{e}^{\mathrm{i}\frac{{\boldsymbol{x}}}{\varepsilon}\cdot{\boldsymbol{p}}}\widehat{{\boldsymbol{V}}}({\boldsymbol{p}}){\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]({\boldsymbol{s}},{\boldsymbol{k}}-{\boldsymbol{p}},\omega){\mathrm{d}}{\boldsymbol{p}}\,,\\\
{\boldsymbol{W}}_{\varepsilon}\left[{\boldsymbol{u}}_{\varepsilon},{\boldsymbol{L}}_{\text{0.5}}\left(\frac{{\boldsymbol{x}}}{\varepsilon},\varepsilon\mathrm{D}_{t}\right){\boldsymbol{u}}_{\varepsilon}\right]&=(\omega-\varepsilon\mathrm{D}_{t}){\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]({\boldsymbol{s}},{\boldsymbol{\xi}}){{\boldsymbol{V}}\left(\frac{{\boldsymbol{x}}}{\varepsilon}\right)}^{*}\,.\end{split}$

Considering again the rules (29) and (30) one obtains that the Wigner
transform ${\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]$
satisfies now:

(60)
${\boldsymbol{L}}_{0}({\boldsymbol{x}},{\boldsymbol{k}}){\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]-{({\boldsymbol{L}}_{0}({\boldsymbol{x}},{\boldsymbol{k}}-\varepsilon{\bf
D}_{\boldsymbol{x}}){\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}])}^{*}\\\
+\sqrt{\varepsilon}\left(\frac{\omega}{(2\pi)^{3}}\int_{\mathbb{R}^{3}}\operatorname{e}^{\mathrm{i}\frac{{\boldsymbol{x}}}{\varepsilon}\cdot{\boldsymbol{p}}}\widehat{{\boldsymbol{V}}}({\boldsymbol{p}}){\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]({\boldsymbol{s}},{\boldsymbol{k}}-{\boldsymbol{p}},\omega){\mathrm{d}}{\boldsymbol{p}}-\omega{\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]{{\boldsymbol{V}}\left(\frac{{\boldsymbol{x}}}{\varepsilon}\right)}^{*}\right)\\\
+\frac{\varepsilon}{\mathrm{i}}\left(\partial_{t}{\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]-{\boldsymbol{\nabla}}_{\boldsymbol{x}}{\boldsymbol{L}}_{0}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}{\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]-({\boldsymbol{\nabla}}_{\boldsymbol{x}}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}{\boldsymbol{L}}_{0}){\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]+{\boldsymbol{\nabla}}_{\boldsymbol{x}}{\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}{{\boldsymbol{L}}}^{*}_{0}\right)\\\
+\frac{\varepsilon}{\mathrm{i}}{\boldsymbol{L}}_{1}{\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]+\frac{\varepsilon}{\mathrm{i}}{\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]{{\boldsymbol{L}}}^{*}_{1}=\mathrm{O}(\varepsilon^{\frac{3}{2}})\,.$

###  4.3. Multiscale expansion

Because of the two-scale dependency of $\mathcal{L}_{\varepsilon}$ in terms of
the slow variable ${\boldsymbol{x}}$ and the fast variable
${\boldsymbol{y}}=\frac{{\boldsymbol{x}}}{\varepsilon}$, we expand
${\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]$ as:

(61)
${\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]\left({\boldsymbol{s}},{\boldsymbol{\xi}}\right)={\boldsymbol{W}}_{0}({\boldsymbol{s}},{\boldsymbol{\xi}})+\sqrt{\varepsilon}{\boldsymbol{W}}_{\text{0.5}}({\boldsymbol{s}},{\boldsymbol{y}},{\boldsymbol{\xi}})+\varepsilon{\boldsymbol{W}}_{1}({\boldsymbol{s}},{\boldsymbol{y}},{\boldsymbol{\xi}})+\mathrm{o}(\varepsilon)\,,$

and replace $\varepsilon{\bf D}_{\boldsymbol{x}}$ by $\varepsilon{\bf
D}_{\boldsymbol{x}}+{\bf D}_{\boldsymbol{y}}$ in Eq. (60). The $0$–th order
equation (36) is recovered for ${\boldsymbol{W}}_{0}$ which is independent of
the fast variable ${\boldsymbol{y}}$, so that it is given by Eq. (16) alike:

(62)
${\boldsymbol{W}}_{0}({\boldsymbol{s}},{\boldsymbol{k}},\omega)=\sum_{\alpha}\delta(\omega+\omega_{\alpha}({\boldsymbol{x}},{\boldsymbol{k}})){\boldsymbol{W}}_{\alpha}({\boldsymbol{s}},{\boldsymbol{k}})\,.$

The $\mathrm{O}(\varepsilon^{\frac{1}{2}})$ half-order terms yield:

(63)
${\boldsymbol{L}}_{0}({\boldsymbol{x}},{\boldsymbol{k}}){\boldsymbol{W}}_{\text{0.5}}({\boldsymbol{s}},{\boldsymbol{y}},{\boldsymbol{\xi}})-{({\boldsymbol{L}}_{0}({\boldsymbol{x}},{\boldsymbol{k}}-{\bf
D}_{\boldsymbol{y}}){\boldsymbol{W}}_{\text{0.5}})}^{*}({\boldsymbol{s}},{\boldsymbol{y}},{\boldsymbol{\xi}})\\\
+\frac{\omega}{(2\pi)^{3}}\int_{\mathbb{R}^{3}}\operatorname{e}^{\mathrm{i}{\boldsymbol{y}}\cdot{\boldsymbol{p}}}\widehat{{\boldsymbol{V}}}({\boldsymbol{p}}){\boldsymbol{W}}_{0}({\boldsymbol{s}},{\boldsymbol{k}}-{\boldsymbol{p}},\omega){\mathrm{d}}{\boldsymbol{p}}-\omega{\boldsymbol{W}}_{0}({\boldsymbol{s}},{\boldsymbol{k}},\omega){{\boldsymbol{V}}({\boldsymbol{y}})}^{*}={\bf
0}\,.$

Taking the Fourier transform with respect to ${\boldsymbol{y}}$ of the above
equation yields:

(64)
${\boldsymbol{L}}_{0}({\boldsymbol{x}},{\boldsymbol{k}})\widehat{{\boldsymbol{W}}}_{\text{0.5}}({\boldsymbol{s}},{\boldsymbol{p}},{\boldsymbol{k}},\omega)-\widehat{{\boldsymbol{W}}}_{\text{0.5}}({\boldsymbol{s}},{\boldsymbol{p}},{\boldsymbol{k}},\omega){{\boldsymbol{L}}}^{*}_{0}({\boldsymbol{x}},{\boldsymbol{k}}-{\boldsymbol{p}})\\\
+\omega\widehat{{\boldsymbol{V}}}({\boldsymbol{p}}){\boldsymbol{W}}_{0}({\boldsymbol{s}},{\boldsymbol{k}}-{\boldsymbol{p}},\omega)-\omega{\boldsymbol{W}}_{0}({\boldsymbol{s}},{\boldsymbol{k}},\omega){\widehat{{\boldsymbol{V}}}(-{\boldsymbol{p}})}^{*}={\bf
0}\,.$

The corrector $\widehat{{\boldsymbol{W}}}_{\text{0.5}}$ is then expanded as
(we remove the ${\boldsymbol{s}}$ and $\omega$ dependences from now on for
clarity purposes):

$\widehat{{\boldsymbol{W}}}_{\text{0.5}}({\boldsymbol{p}},{\boldsymbol{k}})=\sum_{\alpha,\beta}{\boldsymbol{b}}_{\alpha}({\boldsymbol{k}})\widehat{{\bf
w}}_{\alpha\beta}({\boldsymbol{p}},{\boldsymbol{k}}){{\boldsymbol{b}}}^{*}_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}})\,,$

where from Eq. (64) one deduces that:

$\widehat{{\bf
w}}_{\alpha\beta}({\boldsymbol{p}},{\boldsymbol{k}})={{\boldsymbol{c}}}^{*}_{\alpha}({\boldsymbol{k}})\frac{\omega_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}})\widehat{{\boldsymbol{V}}}({\boldsymbol{p}}){\boldsymbol{W}}_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}}){\boldsymbol{W}}_{\alpha}({\boldsymbol{k}}){\widehat{{\boldsymbol{V}}}(-{\boldsymbol{p}})}^{*}}{\omega_{\alpha}({\boldsymbol{k}})-\omega_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}})-\mathrm{i}\theta}{\boldsymbol{c}}_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}})\,.$

Here it is implicitly understood from Eq. (62) that
$\omega=-\omega_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}})$ on the support of
${\boldsymbol{W}}_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}})$ and that
$\omega=-\omega_{\alpha}({\boldsymbol{k}})$ on the support of
${\boldsymbol{W}}_{\alpha}({\boldsymbol{k}})$, while these terms vanish away
from these supports. Also the regularization parameter $\theta\in\mathbb{R}$
is introduced to evade the case
$\omega_{\alpha}({\boldsymbol{k}})=\omega_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}})$
for the time being. It will eventually be sent to $0$ at the end of the
derivation. Therefore one has:222This result should be compared with [18, Eq.
(35)] or [19, Eq. (35)] where no coupling between modes is considered.

(65)
$\widehat{{\boldsymbol{W}}}_{\text{0.5}}({\boldsymbol{p}},{\boldsymbol{k}})=\\\
\sum_{\alpha,\beta}\frac{\omega_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}})\boldsymbol{\Pi}_{\alpha}({\boldsymbol{k}})\widehat{{\boldsymbol{V}}}({\boldsymbol{p}}){\boldsymbol{W}}_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}}){\boldsymbol{W}}_{\alpha}({\boldsymbol{k}}){\widehat{{\boldsymbol{V}}}(-{\boldsymbol{p}})}^{*}{\boldsymbol{\Pi}_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}})}^{*}}{\omega_{\alpha}({\boldsymbol{k}})-\omega_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}})-\mathrm{i}\theta}\,.$

At last, the $\mathrm{O}(\varepsilon)$ terms yield:

$\partial_{t}{\boldsymbol{W}}_{0}-{\boldsymbol{\nabla}}_{\boldsymbol{x}}{\boldsymbol{L}}_{0}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}{\boldsymbol{W}}_{0}-({\boldsymbol{\nabla}}_{\boldsymbol{x}}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}{\boldsymbol{L}}_{0}){\boldsymbol{W}}_{0}+{\boldsymbol{\nabla}}_{\boldsymbol{x}}{\boldsymbol{W}}_{0}\cdot{\boldsymbol{\nabla}}_{\boldsymbol{k}}{{\boldsymbol{L}}}^{*}_{0}+{\boldsymbol{L}}_{1}{\boldsymbol{W}}_{0}+{\boldsymbol{W}}_{0}{{\boldsymbol{L}}}^{*}_{1}\\\
+\frac{\mathrm{i}\omega}{(2\pi)^{3}}\int_{\mathbb{R}^{3}}\operatorname{e}^{\mathrm{i}{\boldsymbol{y}}\cdot{\boldsymbol{p}}}\widehat{{\boldsymbol{V}}}({\boldsymbol{p}}){\boldsymbol{W}}_{\text{0.5}}({\boldsymbol{y}},{\boldsymbol{k}}-{\boldsymbol{p}}){\mathrm{d}}{\boldsymbol{p}}-\mathrm{i}\omega{\boldsymbol{W}}_{\text{0.5}}({\boldsymbol{y}},{\boldsymbol{k}}){{\boldsymbol{V}}({\boldsymbol{y}})}^{*}\\\
+\mathrm{i}{\boldsymbol{L}}_{0}({\boldsymbol{k}}){\boldsymbol{W}}_{1}+{(\mathrm{i}{\boldsymbol{L}}_{0}({\boldsymbol{k}}-{\bf
D}_{\boldsymbol{y}}){\boldsymbol{W}}_{1})}^{*}={\bf 0}\,.$

We plug Eq. (65) into the above equation, multiply it on the left by
$\smash{{{\boldsymbol{c}}}^{*}_{\alpha}({\boldsymbol{k}})}$ and on the right
by $\smash{{\boldsymbol{c}}_{\alpha}({\boldsymbol{k}})}$, and finally take the
average ${\mathbb{E}}\left\\{\cdot\right\\}$ with respect to the random
perturbations ${\boldsymbol{V}}$ with due consideration of Eq. (53). The terms
in ${\boldsymbol{W}}_{1}$ vanish once they are averaged since
$\smash{{\mathbb{E}}\left\\{{\boldsymbol{W}}_{1}\right\\}}={\bf 0}$ for
$\smash{{\mathbb{E}}\left\\{{\boldsymbol{W}}_{\varepsilon}[{\boldsymbol{u}}_{\varepsilon}]\right\\}}=\smash{{\mathbb{E}}\left\\{{\boldsymbol{W}}_{0}\right\\}}$
by construction, and one gets:

(66) $\partial_{t}{\mathbb{E}}\left\\{{\bf
w}_{\alpha}\right\\}+\\{\omega_{\alpha},{\mathbb{E}}\left\\{{\bf
w}_{\alpha}\right\\}\\}+{\boldsymbol{\ell}}_{\alpha}{\mathbb{E}}\left\\{{\bf
w}_{\alpha}\right\\}+{\mathbb{E}}\left\\{{\bf
w}_{\alpha}\right\\}{{\boldsymbol{\ell}}}^{*}_{\alpha}\\\
+{\boldsymbol{n}}_{\alpha}{\mathbb{E}}\left\\{{\bf
w}_{\alpha}\right\\}-{\mathbb{E}}\left\\{{\bf
w}_{\alpha}\right\\}{\boldsymbol{n}}_{\alpha}={\mathbb{E}}\left\\{{\mathcal{I}}_{1}\right\\}-{\mathbb{E}}\left\\{{\mathcal{I}}_{2}\right\\}\,,$

where:

(67)
$\begin{split}{\mathcal{I}}_{1}&=\mathrm{i}\omega{{\boldsymbol{c}}}^{*}_{\alpha}({\boldsymbol{k}}){\boldsymbol{W}}_{\text{0.5}}({\boldsymbol{y}},{\boldsymbol{k}}){{\boldsymbol{V}}}^{*}({\boldsymbol{y}}){\boldsymbol{c}}_{\alpha}({\boldsymbol{k}})\\\
&=\sum_{\beta}\iint\operatorname{e}^{\mathrm{i}{\boldsymbol{y}}\cdot({\boldsymbol{p}}-{\boldsymbol{q}})}\frac{\omega_{\beta}^{2}({\boldsymbol{k}}-{\boldsymbol{p}}){\bf
s}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{p}},{\boldsymbol{k}}-{\boldsymbol{p}}){\bf
w}_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}}){\bf
s}_{\beta\alpha}({\boldsymbol{k}}-{\boldsymbol{p}},-{\boldsymbol{q}},{\boldsymbol{k}})}{\mathrm{i}(\omega_{\alpha}({\boldsymbol{k}})-\omega_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}}))+\theta}\frac{{\mathrm{d}}{\boldsymbol{p}}{\mathrm{d}}{\boldsymbol{q}}}{(2\pi)^{6}}\\\
&\quad-\omega_{\alpha}^{2}({\boldsymbol{k}}){\bf
w}_{\alpha}({\boldsymbol{k}}){\left(\sum_{\beta}\iint\operatorname{e}^{-\mathrm{i}{\boldsymbol{y}}\cdot({\boldsymbol{p}}-{\boldsymbol{q}})}\frac{{\bf
s}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{q}},{\boldsymbol{k}}-{\boldsymbol{p}}){\bf
s}_{\beta\alpha}({\boldsymbol{k}}-{\boldsymbol{p}},-{\boldsymbol{p}},{\boldsymbol{k}})}{\mathrm{i}(\omega_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}}))+\theta}\frac{{\mathrm{d}}{\boldsymbol{p}}{\mathrm{d}}{\boldsymbol{q}}}{(2\pi)^{6}}\right)}^{*}\,;\end{split}$

and:

(68)
${\scriptsize\begin{split}{\mathcal{I}}_{2}&=\mathrm{i}\omega{{\boldsymbol{c}}}^{*}_{\alpha}({\boldsymbol{k}})\left(\iint\operatorname{e}^{\mathrm{i}{\boldsymbol{y}}\cdot({\boldsymbol{p}}+{\boldsymbol{q}})}\widehat{{\boldsymbol{V}}}({\boldsymbol{p}})\widehat{{\boldsymbol{W}}}_{\text{0.5}}({\boldsymbol{q}},{\boldsymbol{k}}-{\boldsymbol{p}})\frac{{\mathrm{d}}{\boldsymbol{p}}{\mathrm{d}}{\boldsymbol{q}}}{(2\pi)^{6}}\right){\boldsymbol{c}}_{\alpha}({\boldsymbol{k}})\\\
&=\sum_{\gamma,\beta}\left(\iint\operatorname{e}^{\mathrm{i}{\boldsymbol{y}}\cdot({\boldsymbol{p}}+{\boldsymbol{q}})}\frac{\omega_{\beta}^{2}({\boldsymbol{k}}-{\boldsymbol{p}}-{\boldsymbol{q}}){\bf
s}_{\alpha\gamma}({\boldsymbol{k}},{\boldsymbol{p}},{\boldsymbol{k}}-{\boldsymbol{p}}){\bf
s}_{\gamma\beta}({\boldsymbol{k}}-{\boldsymbol{p}},{\boldsymbol{q}},{\boldsymbol{k}}-{\boldsymbol{p}}-{\boldsymbol{q}}){\bf
w}_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}}-{\boldsymbol{q}}){{\boldsymbol{b}}}^{*}_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}}-{\boldsymbol{q}})}{\mathrm{i}(\omega_{\gamma}({\boldsymbol{k}}-{\boldsymbol{p}})-\omega_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}}-{\boldsymbol{q}}))+\theta}\frac{{\mathrm{d}}{\boldsymbol{p}}{\mathrm{d}}{\boldsymbol{q}}}{(2\pi)^{6}}\right){\boldsymbol{c}}_{\alpha}({\boldsymbol{k}})\\\
&\quad-\sum_{\gamma,\beta}\left(\iint\operatorname{e}^{\mathrm{i}{\boldsymbol{y}}\cdot({\boldsymbol{p}}+{\boldsymbol{q}})}\frac{\omega_{\gamma}^{2}({\boldsymbol{k}}-{\boldsymbol{p}}){\bf
s}_{\alpha\gamma}({\boldsymbol{k}},{\boldsymbol{p}},{\boldsymbol{k}}-{\boldsymbol{p}}){\bf
w}_{\gamma}({\boldsymbol{k}}-{\boldsymbol{p}}){\bf
s}_{\gamma\beta}({\boldsymbol{k}}-{\boldsymbol{p}},{\boldsymbol{q}},{\boldsymbol{k}}-{\boldsymbol{p}}-{\boldsymbol{q}}){{\boldsymbol{b}}}^{*}_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}}-{\boldsymbol{q}})}{\mathrm{i}(\omega_{\gamma}({\boldsymbol{k}}-{\boldsymbol{p}})-\omega_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}}-{\boldsymbol{q}}))+\theta}\frac{{\mathrm{d}}{\boldsymbol{p}}{\mathrm{d}}{\boldsymbol{q}}}{(2\pi)^{6}}\right){\boldsymbol{c}}_{\alpha}({\boldsymbol{k}})\,.\end{split}}$

In the above relationships we have introduced:

${\bf
s}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{p}},{\boldsymbol{q}}):={{\boldsymbol{c}}}^{*}_{\alpha}({\boldsymbol{k}})\widehat{{\boldsymbol{V}}}({\boldsymbol{p}}){\boldsymbol{b}}_{\beta}({\boldsymbol{q}})={{\bf
s}_{\beta\alpha}({\boldsymbol{q}},-{\boldsymbol{p}},{\boldsymbol{k}})}^{*}\,,$

where the second equality stems from
${\boldsymbol{K}}_{0}{\boldsymbol{V}}={{\boldsymbol{V}}}^{*}{\boldsymbol{K}}_{0}$
by the Hermiticity of the optical response. By the change of variable
$-{\boldsymbol{q}}\to{\boldsymbol{q}}$, the use of Eq. (53), and then the
change of variable ${\boldsymbol{k}}-{\boldsymbol{p}}\to{\boldsymbol{p}}$, the
average of the first sum in Eq. (67) reads:

$\begin{split}{\mathbb{E}}\left\\{{\mathcal{I}}_{11}\right\\}&=\sum_{\beta}\iint\operatorname{e}^{\mathrm{i}{\boldsymbol{y}}\cdot({\boldsymbol{p}}+{\boldsymbol{q}})}\frac{\omega_{\beta}^{2}({\boldsymbol{k}}-{\boldsymbol{p}}){\mathbb{E}}\left\\{{\bf
s}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{p}},{\boldsymbol{k}}-{\boldsymbol{p}}){\bf
w}_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}}){\bf
s}_{\beta\alpha}({\boldsymbol{k}}-{\boldsymbol{p}},{\boldsymbol{q}},{\boldsymbol{k}})\right\\}}{\mathrm{i}(\omega_{\alpha}({\boldsymbol{k}})-\omega_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}}))+\theta}\frac{{\mathrm{d}}{\boldsymbol{p}}{\mathrm{d}}{\boldsymbol{q}}}{(2\pi)^{6}}\\\
&=\sum_{\beta}\int\frac{\omega_{\beta}^{2}({\boldsymbol{k}}-{\boldsymbol{p}})\widehat{{\boldsymbol{R}}}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{k}}-{\boldsymbol{p}}):{\mathbb{E}}\left\\{{\bf
w}_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}})\right\\}}{\mathrm{i}(\omega_{\alpha}({\boldsymbol{k}})-\omega_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}}))+\theta}{\mathrm{d}}{\boldsymbol{p}}\\\
&=\sum_{\beta}\int\frac{\omega_{\beta}^{2}({\boldsymbol{p}})\widehat{{\boldsymbol{R}}}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{p}}):{\mathbb{E}}\left\\{{\bf
w}_{\beta}({\boldsymbol{p}})\right\\}}{\mathrm{i}(\omega_{\alpha}({\boldsymbol{k}})-\omega_{\beta}({\boldsymbol{p}}))+\theta}{\mathrm{d}}{\boldsymbol{p}}\\\
\end{split}$

where the linear operator $\smash{\widehat{{\boldsymbol{R}}}_{\alpha\beta}}$
is:

$[\widehat{{\boldsymbol{R}}}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{p}})]_{aa^{\prime}bb^{\prime}}:=\overline{{\boldsymbol{c}}_{\alpha}^{a}({\boldsymbol{k}})}\otimes{\boldsymbol{b}}_{\beta}^{b}({\boldsymbol{p}}):\widehat{{\boldsymbol{R}}}({\boldsymbol{k}}-{\boldsymbol{p}}):\overline{{\boldsymbol{c}}_{\beta}^{b^{\prime}}({\boldsymbol{p}})}\otimes{\boldsymbol{b}}_{\alpha}^{a^{\prime}}({\boldsymbol{k}})\,,$

for $1\leq a,a^{\prime}\leq A$ and $1\leq b,b^{\prime}\leq B$, $A$ and $B$
being the orders of multiplicity of the modes $\alpha$ and $\beta$,
respectively.333This result should be compared with [18, Eq. (42)] or [19, Eq.
(42)] where no coupling between modes occurs and this operator rather reads:
$[\widehat{{\boldsymbol{R}}}_{\alpha}({\boldsymbol{k}},{\boldsymbol{p}})]_{aa^{\prime}bb^{\prime}}:=\overline{{\boldsymbol{c}}_{\alpha}^{a}({\boldsymbol{k}})}\otimes{\boldsymbol{b}}_{\alpha}^{b}({\boldsymbol{p}}):\widehat{{\boldsymbol{R}}}({\boldsymbol{k}}-{\boldsymbol{p}}):{\boldsymbol{c}}_{\alpha}^{a^{\prime}}({\boldsymbol{k}})\otimes\overline{{\boldsymbol{b}}_{\alpha}^{b^{\prime}}({\boldsymbol{p}})}\,.$
Here we have used a mixing assumption by which expectations over products of
$\smash{\widehat{{\boldsymbol{V}}}}$ and over $\smash{{\bf w}_{\beta}}$ are
independent and get decoupled, because both quantities vary on different
scales. Now using the changes of variable
$-{\boldsymbol{p}}\to{\boldsymbol{p}}$, and then
${\boldsymbol{k}}+{\boldsymbol{p}}\to{\boldsymbol{p}}$, the average of the
second sum in Eq. (67) reads:

(69)
$\begin{split}{\mathbb{E}}\left\\{{\mathcal{I}}_{12}\right\\}&=\omega_{\alpha}^{2}({\boldsymbol{k}})\sum_{\beta}\iint\operatorname{e}^{\mathrm{i}{\boldsymbol{y}}\cdot({\boldsymbol{q}}-{\boldsymbol{p}})}\frac{{\mathbb{E}}\left\\{{\bf
s}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{q}},{\boldsymbol{k}}-{\boldsymbol{p}}){\bf
s}_{\beta\alpha}({\boldsymbol{k}}-{\boldsymbol{p}},-{\boldsymbol{p}},{\boldsymbol{k}})\right\\}}{\mathrm{i}(\omega_{\beta}({\boldsymbol{k}}-{\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}}))+\theta}\frac{{\mathrm{d}}{\boldsymbol{p}}{\mathrm{d}}{\boldsymbol{q}}}{(2\pi)^{6}}\\\
&=\omega_{\alpha}^{2}({\boldsymbol{k}})\sum_{\beta}\iint\operatorname{e}^{\mathrm{i}{\boldsymbol{y}}\cdot({\boldsymbol{q}}+{\boldsymbol{p}})}\frac{{\mathbb{E}}\left\\{{\bf
s}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{q}},{\boldsymbol{k}}+{\boldsymbol{p}}){\bf
s}_{\beta\alpha}({\boldsymbol{k}}+{\boldsymbol{p}},{\boldsymbol{p}},{\boldsymbol{k}})\right\\}}{\mathrm{i}(\omega_{\beta}({\boldsymbol{k}}+{\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}}))+\theta}\frac{{\mathrm{d}}{\boldsymbol{p}}{\mathrm{d}}{\boldsymbol{q}}}{(2\pi)^{6}}\\\
&=\sum_{\beta}\int\frac{\omega_{\alpha}^{2}({\boldsymbol{k}})\widehat{{\boldsymbol{R}}}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{k}}+{\boldsymbol{p}}):{\boldsymbol{I}}_{B}}{\mathrm{i}(\omega_{\beta}({\boldsymbol{k}}+{\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}}))+\theta}{\mathrm{d}}{\boldsymbol{p}}\\\
&=\sum_{\beta}\int\frac{\omega_{\alpha}^{2}({\boldsymbol{k}})\widehat{{\boldsymbol{R}}}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{p}}):{\boldsymbol{I}}_{B}}{\mathrm{i}(\omega_{\beta}({\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}}))+\theta}{\mathrm{d}}{\boldsymbol{p}}\,.\end{split}$

Besides, taking the average of Eq. (68), one has with Eq. (53):

$\begin{split}{\mathbb{E}}\left\\{{\mathcal{I}}_{2}\right\\}&=\left(\sum_{\gamma}\int\frac{\omega_{\alpha}^{2}({\boldsymbol{k}})\widehat{{\boldsymbol{R}}}_{\alpha\gamma}({\boldsymbol{k}},{\boldsymbol{k}}-{\boldsymbol{p}}):{\boldsymbol{I}}_{G}}{\mathrm{i}(\omega_{\gamma}({\boldsymbol{k}}-{\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}}))+\theta}{\mathrm{d}}{\boldsymbol{p}}\right){\mathbb{E}}\left\\{{\bf
w}_{\alpha}({\boldsymbol{k}})\right\\}\\\
&\quad\quad-\sum_{\gamma}\int\frac{\omega_{\gamma}^{2}({\boldsymbol{k}}-{\boldsymbol{p}})\widehat{{\boldsymbol{R}}}_{\alpha\gamma}({\boldsymbol{k}},{\boldsymbol{k}}-{\boldsymbol{p}}):{\mathbb{E}}\left\\{{\bf
w}_{\gamma}({\boldsymbol{k}}-{\boldsymbol{p}})\right\\}}{\mathrm{i}(\omega_{\gamma}({\boldsymbol{k}}-{\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}}))+\theta}{\mathrm{d}}{\boldsymbol{p}}\\\
&=\left(\sum_{\gamma}\int\frac{\omega_{\alpha}^{2}({\boldsymbol{k}})\widehat{{\boldsymbol{R}}}_{\alpha\gamma}({\boldsymbol{k}},{\boldsymbol{p}}):{\boldsymbol{I}}_{G}}{\mathrm{i}(\omega_{\gamma}({\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}}))+\theta}{\mathrm{d}}{\boldsymbol{p}}\right){\mathbb{E}}\left\\{{\bf
w}_{\alpha}({\boldsymbol{k}})\right\\}\\\
&\quad\quad-\sum_{\gamma}\int\frac{\omega_{\gamma}^{2}({\boldsymbol{p}})\widehat{{\boldsymbol{R}}}_{\alpha\gamma}({\boldsymbol{k}},{\boldsymbol{p}}):{\mathbb{E}}\left\\{{\bf
w}_{\gamma}({\boldsymbol{p}})\right\\}}{\mathrm{i}(\omega_{\gamma}({\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}}))+\theta}{\mathrm{d}}{\boldsymbol{p}}\,,\end{split}$

where $G$ stands for the order of multiplicity of the mode $\gamma$.

The last step is to let $\theta\to 0$, which is done thank to the (formal)
identity
$\smash{\frac{1}{\mathrm{i}x+\theta}}\to\smash{\frac{1}{\mathrm{i}x}}+\hat{\theta}\pi\delta(x)$,
for $\smash{\hat{\theta}}$ being the sign of $\theta$. This yields:

${\mathbb{E}}\left\\{{\mathcal{I}}_{1}\right\\}-{\mathbb{E}}\left\\{{\mathcal{I}}_{2}\right\\}=\sum_{\beta}\int
2\pi\hat{\theta}\delta(\omega_{\alpha}({\boldsymbol{k}})-\omega_{\beta}({\boldsymbol{p}}))\omega_{\beta}^{2}({\boldsymbol{p}})\widehat{{\boldsymbol{R}}}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{p}}):{\mathbb{E}}\left\\{{\bf
w}_{\beta}({\boldsymbol{p}})\right\\}{\mathrm{d}}{\boldsymbol{p}}\\\
-{\mathbb{E}}\left\\{{\bf
w}_{\alpha}({\boldsymbol{k}})\right\\}{\left[\sum_{\beta}\int\left(\frac{1}{\mathrm{i}(\omega_{\beta}({\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}}))}+\pi\hat{\theta}\delta(\omega_{\beta}({\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}}))\right)\omega_{\alpha}^{2}({\boldsymbol{k}})\widehat{{\boldsymbol{R}}}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{p}}):{\boldsymbol{I}}_{B}{\mathrm{d}}{\boldsymbol{p}}\right]}^{*}\\\
-\left[\sum_{\beta}\int\left(\frac{1}{\mathrm{i}(\omega_{\beta}({\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}}))}+\pi\hat{\theta}\delta(\omega_{\beta}({\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}}))\right)\omega_{\alpha}^{2}({\boldsymbol{k}})\widehat{{\boldsymbol{R}}}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{p}}):{\boldsymbol{I}}_{B}{\mathrm{d}}{\boldsymbol{p}}\right]{\mathbb{E}}\left\\{{\bf
w}_{\alpha}({\boldsymbol{k}})\right\\}\,.$

Consider $\hat{\theta}=1$ (for causality purposes) and let:

(70)
${\boldsymbol{\sigma}}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{p}})=2\pi\omega_{\alpha}({\boldsymbol{k}})\omega_{\beta}({\boldsymbol{p}})\widehat{{\boldsymbol{R}}}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{p}})\,,$

and:

(71)
${\boldsymbol{\Sigma}}_{\alpha}({\boldsymbol{k}})=\frac{1}{2}\sum_{\beta}\int\delta(\omega_{\beta}({\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}})){\boldsymbol{\sigma}}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{p}}):{\boldsymbol{I}}_{B}{\mathrm{d}}{\boldsymbol{p}}-\frac{\mathrm{i}}{2\pi}\sum_{\beta}\oint\frac{{\boldsymbol{\sigma}}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{p}}):{\boldsymbol{I}}_{B}}{\omega_{\beta}({\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}})}{\mathrm{d}}{\boldsymbol{p}}\,,$

where the second integral holds as a Cauchy principal value; then Eq. (66)
finally reads as the radiative transfer equations:

(72) $\partial_{t}{\mathbb{E}}\left\\{{\bf
w}_{\alpha}\right\\}+\\{\omega_{\alpha},{\mathbb{E}}\left\\{{\bf
w}_{\alpha}\right\\}\\}+{\boldsymbol{\ell}}_{\alpha}{\mathbb{E}}\left\\{{\bf
w}_{\alpha}\right\\}+{\mathbb{E}}\left\\{{\bf
w}_{\alpha}\right\\}{{\boldsymbol{\ell}}}^{*}_{\alpha}+{\boldsymbol{n}}_{\alpha}{\mathbb{E}}\left\\{{\bf
w}_{\alpha}\right\\}-{\mathbb{E}}\left\\{{\bf
w}_{\alpha}\right\\}{\boldsymbol{n}}_{\alpha}\\\
=\sum_{\beta}\int\delta(\omega_{\beta}({\boldsymbol{p}})-\omega_{\alpha}({\boldsymbol{k}})){\boldsymbol{\sigma}}_{\alpha\beta}({\boldsymbol{k}},{\boldsymbol{p}}):{\mathbb{E}}\left\\{{\bf
w}_{\beta}({\boldsymbol{p}})\right\\}{\mathrm{d}}{\boldsymbol{p}}\\\
-{\boldsymbol{\Sigma}}_{\alpha}({\boldsymbol{k}}){\mathbb{E}}\left\\{{\bf
w}_{\alpha}({\boldsymbol{k}})\right\\}-{\mathbb{E}}\left\\{{\bf
w}_{\alpha}({\boldsymbol{k}})\right\\}{{\boldsymbol{\Sigma}}_{\alpha}({\boldsymbol{k}})}^{*}\,,$

possibly coupling all modes of propagation while keeping the frequency
constant: $\omega_{\beta}({\boldsymbol{p}})=\omega_{\alpha}({\boldsymbol{k}})$
in the scattering processes. This result extends [46, Eq. (4.32)] for
isotropic dielectric media to the case of general bianisotropic dielectric
media.444The sole difference lies in the sign of the imaginary part of
${\boldsymbol{\Sigma}}_{\alpha}$ and the $1/\pi$ factor, to be compared with
[46, Eq. (4.34)]. It revisits the derivation in [18, 19] by considering the
time-dependent case with damping effects, and possible coupling between the
modes. Moreover, it uses Wigner transforms and pseudo-differential calculus in
the same standard quantization rather than mixing it with the Weyl
quantization. Matrix-valued radiative transfer equations of the form (72) can
be solved numerically by the Monte-Carlo method, which is based on their
probabilistic representation in terms of a jump Markov process as outlined in
[40]. It extends ray tracing to scattering media.

######  Remark 3.

In [12] the derivation of a radiative transfer equation of the form of Eq.
(72) is proved rigorously for the scalar wave equation with a random speed of
sound. Here it is shown that the average Wigner transform of the wave field
converges (weakly) to the solution of a linear Boltzmann equation similar to
(72), provided that the initial conditions satisfy some _adhoc_ boundedness
and tightness conditions and that the random fluctuations fulfill some
regularity assumptions as well. A self-averaging property is also
demonstrated, by which the deviations of the Wigner transform from its average
are shown to vanish in the high-frequency limit. This means that Eq. (72)
would also hold for ${\bf w}_{\alpha}$ and not only ${\mathbb{E}}\left\\{{\bf
w}_{\alpha}\right\\}$, a result that is implicitly acknowledged in the
statement of Eq. (22). The proofs in [12] are however quite involved and very
technical so that their extension to the present case is much beyond the scope
of this research. Radiative transfer equations can also be rigorously derived
for some particular configurations as in [16, 36].

##  5\. Summary and outlook

We have derived a system of coupled radiative transfer equations to describe
the propagation of high-frequency electromagnetic waves in randomly
fluctuating bianisotropic dielectric media including dispersive and
dissipative effects. The fluctuations are weak and vary at the same length
scales as the typical wavelength. We have used a Wigner functional approach
and its interpretation in terms of semiclassical pseudo-differential operators
to derive the transfer equations from a multiscale analysis. Radiative
transfer equations have a geometrical interpretation in terms of
bicharacteristic rays that is well suited for their numerical integration by
ray methods and ray tracing solvers. Our results extend this approach to
heterogeneous bianisotropic media or plasmas beyond the classical applications
of ray tracing in homogeneous or piecewise homogeneous isotropic media. In
future works one aims to derive the diffusion limit of the radiative transfer
model, which describes the evolution of the energy density in physical space
solely after the waves have travelled several scattering mean free paths and
have lost the memory of their angular distribution at earlier times. The issue
of deriving boundary conditions adapted to the Wigner distribution from the
boundary conditions applied to the electromagnetic fields is also worth
considering, for applications in radar imaging or remote sensing with scenes
exhibiting sharp interfaces for example.

## Acknowledgement

We thank Thomas Lepetit at ONERA for valuable discussions.

## References

  * [1] J.-L. Akian. Space-time semiclassical measures for three-dimensional elastodynamics: boundary conditions for the hyperbolic set. Asymp. Anal., 78 (2012), pp. 37–83.
  * [2] J.-L. Akian and É. Savin. Kinetic modeling of multiple scattering of acoustic waves in randomly heterogeneous flows. SIAM Multiscale Model. Simul., 19 (2021), pp. 1394–1424.
  * [3] E. Akkermans and G. Montambaux. Mesoscopic Physics of Electrons and Photons. Cambridge University Press, Cambridge, 2007.
  * [4] R. Alonso and L. Borcea. Electromagnetic wave propagation in random waveguides. SIAM Multiscale Model. Simul., 13 (2021), pp. 847–889.
  * [5] A. Atallah-Baraket and C. Fermanian Kammerer. High frequency analysis of families of solutions to the equation of viscoelasticity of Kelvin-Voigt. J. Hyperbolic Differ. Equ., 1 (2004), pp. 789–812.
  * [6] S. Auer, S. Hinz, and R. Bamler. Ray-tracing simulation techniques for understanding high-resolution SAR images. IEEE Trans. Geosci. Remote Sens., 48 (2010), pp. 1445–1456.
  * [7] V. M. Babic and V. S. Buldyrev. Short-Wavelength Diffraction Theory: Asymptotic Methods. Springer, Berlin, 1991.
  * [8] G. Bal. Kinetics of scalar wave fields in random media. Wave Motion, 43 (2005), pp. 132–157.
  * [9] G. Bal, T. Komorowski, and L. V. Ryzhik. Kinetic limits for waves in a random medium. Kinet. Relat. Models, 3 (2010), pp. 529–644.
  * [10] I. Baydoun, É. Savin, R. Cottereau, D. Clouteau, and J. Guilleminot. Kinetic modeling of multiple scattering of elastic waves in heterogeneous anisotropic media. Wave Motion, 51 (2014), pp. 1325–1348.
  * [11] L. Borcea and J. Garnier. Polarization effects for electromagnetic wave propagation in random media. Wave Motion, 63 (2016), pp. 179–208.
  * [12] M. J. Butz. Kinetic Limit for Wave Propagation in a Continuous, Weakly Random Medium. Self-averaging and Convergence to a Linear Boltzmann Equation. PhD Thesis, Technische Universität München, August 2015.
  * [13] R. Carminati and J. C. Schotland. Principles of Scattering and Transport of Light. Cambridge University Press, Cambridge, 2021.
  * [14] M. Cassier, P. Joly, and M. Kachanovska. Mathematical models for dispersive electromagnetic waves: An overview. Comput. Math. Appl., 74 (2017), pp. 2792–2830.
  * [15] S. Chandrasekhar. Radiative Transfer. Dover, New York NY, 1960.
  * [16] L. Erdös and H.-T. Yau. Linear Boltzmann equation as the weak coupling limit of a random Schrödinger equation. Comm. Pure Appl. Math., LIII (2000), pp. 667–735.
  * [17] L. C. Evans. Partial Differential Equations, 2nd Ed. American Mathematical Society, Providence RI, 2010.
  * [18] A. Fannjiang. Two-frequency radiative transfer: Maxwell equations in random dielectrics. J. Opt. Soc. Am. A, 24 (2007), pp. 3680–3690.
  * [19] A. Fannjiang. Mutual coherence of polarized light in disordered media: two-frequency method extended. J. Phys. A: Math. Theor., 40 (2007), pp. 13667–13683.
  * [20] C. Fermanian-Kammerer and P. Gérard. Mesures semi-classiques et croisement de modes (Semiclassical measures and mode crossing). Bull. Soc. Math. France, 130 (2002), pp. 123–168.
  * [21] S. Filippas and G. E. Makrakis. Semiclassical Wigner function and geometrical optics. SIAM Multiscale Model. Simul., 1 (2003), pp. 674–710.
  * [22] J. Garnier and K. Sølna. Transmission and reflection of electromagnetic waves in randomly layered media. Commun. Math. Sci., 13 (2015), pp. 707–728.
  * [23] P. Gérard, P. A. Markowich, N. J. Mauser, and F. Poupaud. Homogenization limits and Wigner transforms. Comm. Pure Appl. Math., L (1997), pp. 323–379.
  * [24] A. S. Glassner (Ed.). An Introduction to Ray Tracing. Academic Press, San Diego CA, 1989.
  * [25] A. Goldsmith. Wireless Communications. Cambridge University Press, Cambridge, 2005.
  * [26] C. Gomez. Radiative transport limit for the random Schrödinger equation with long-range correlations. J. Math. Pures Appl., 98 (2012), pp. 295–327.
  * [27] M. Guo and X.-P. Wang. Transport equations for a general class of evolution equations with random perturbations. J. Math. Phys., 40 (1999), pp. 4828–4858.
  * [28] E. Haines and T. Akenine-Möller (Eds.). Ray Tracing Gems: High-Quality and Real-Time Rendering with DXR and Other APIs. Apress, Berkeley CA, 2019.
  * [29] A. D. Ioannidis, G. Kristensson, and I. G. Stratis. On the well-posedness of the Maxwell system for linear bianisotropic media. SIAM J. Math. Anal., 44 (2012), pp. 2459–2473.
  * [30] A. Ishimaru. Wave Propagation and Scattering in Random Media, vol. I & II. Academic Press, New York NY, 1978.
  * [31] A. Karlsson and G. Kristensson. Constitutive relations, dissipation, and reciprocity for the Maxwell equations in the time domain. J. Electromagn. Waves Appl., 6 (1992), pp. 537–551.
  * [32] Y. A. Kravtsov and Y. I. Orlov. Geometrical Optics of Inhomogeneous Media. Springer, Berlin, 1990.
  * [33] J. Kraisler, W. Li, K. Ren, J. C. Schotland, and Y. Zhong. Transport models for wave propagation in scattering media with nonlinear absorption. SIAM J. Appl. Math., 83 (2023), pp. .
  * [34] C. Lessig and A. L. Castro. The geometry of radiative transfer. In Geometry, Mechanics, and Dynamics (D. Chang, D. Holm, G. Patrick, and T. Ratiu, eds.), Fields Institute Communications vol. 27, Springer, New York NY (2015), pp. 239–255.
  * [35] P.-L. Lions and T. Paul. Sur les mesures de Wigner (On Wigner measures). Rev. Mat. Iberoamericana, 9 (1993), pp. 553–618.
  * [36] J. Lukkarinen and H. Spohn. Kinetic limit for wave propagation in a random medium. Arch. Rat. Mech. Anal., 183 (2007), pp. 93–162.
  * [37] A. Martinez. An Introduction to Semiclassical and Microlocal Analysis. Springer, Berlin, 2002.
  * [38] M. I. Mishchenko. Electromagnetic Scattering by Particles and Particle Groups. An Introduction. Cambridge University Press, Cambridge, 2014.
  * [39] G. C. Papanicolaou and R. Burridge. Transport equations for the Stokes parameters from Maxwell’s equations in a random medium. J. Math. Phys., 16 (1975), pp. 2074–2085.
  * [40] G. C. Papanicolaou, G. Bal, and L. V. Ryzhik. Probabilistic theory of transport processes with polarization. SIAM J. Appl. Math., 60 (2000), pp 1639–1666.
  * [41] S. G. Parker, J. Bigler, A. Dietrich, H. Friedrich, J. Hoberock, D. Luebke, D. McAllister, M. McGuire, K. Morley, A. Robison, and M. Stich. OptiX: a general purpose ray tracing engine. ACM Trans. Graphics, 29 (2010), Article no. 66.
  * [42] M. Petrini. Behaviour of the energy density associated to a Kelvin-Voight model in viscoelasticity. Asymp. Anal., 34 (2003), pp. 261–273.
  * [43] J. M. Powell and J. Vanneste. Transport equations for waves in randomly perturbed Hamiltonian systems, with application to Rossby waves. Wave Motion, 42 (2005), pp. 289–308.
  * [44] A. Raman and S. Fan. Photonic band structure of dispersive metamaterials formulated as a Hermitian eigenvalue problem. Phys. Rev. Lett., 104 (2010), 087401.
  * [45] M. C. W. van Rossum and Th. M. Nieuwenhuizen. Multiple scattering of classical waves: microscopy, mesoscopy, and diffusion. Rev. Mod. Phys., 71 (1999), pp. 313–371.
  * [46] L. V. Ryzhik, G. C. Papanicolaou, and J. B. Keller. Transport equations for elastic and other waves in random media. Wave Motion, 24 (1996), pp. 327–370.
  * [47] P. Sheng. Introduction to Wave Scattering, Localization and Mesoscopic Phenomena, 2nd Ed. Springer, Berlin, 2006.
  * [48] M. I. Skolnik (Ed.). Radar Handbook, 3rd Ed. McGraw-Hill, New York NY, 2008.
  * [49] C. Sparber, P. A. Markowich, and N. J. Mauser. Wigner functions versus WKB-methods in multivalued geometrical optics. Asymp. Anal., 33 (2003), pp. 153–187.
  * [50] E. R. Tracy, A. J. Brizard, A. S. Richardson, and A. N. Kaufman. Ray Tracing and Beyond. Cambridge University Press, Cambridge, 2014.
  * [51] L. Tsang, J. A. Kong, and K.-H. Ding. Scattering of Electromagnetic Waves: Theories and Applications. John Wiley & Sons, New York NY, 2000.
  * [52] V. G. Veselago The electrodynamics of substances with simultaneously negative values of $\epsilon$ and $\mu$. Sov. Phys. Usp., 10 (1968), pp. 509–514.
  * [53] M. Zworski. Semiclassical Analysis. American Mathematical Society, Providence RI, 2012.