<figure><img src="image_8.png" /><figcaption>Fig. 4.1: Realizations of Gaussian random fields (for $d = 2$) with correlation function (4.3). The left picture represents a field with short-range correlations, while the right picture depicts a field with long-range correlations with $\alpha = 0.999$.</figcaption></figure>

media to describe the energy density propagation [12, 45, 180]. In such a context, the velocity field $c(\mathbf{x})$ has the form

$$c^2(\mathbf{x}) = c_0^2 \left( 1 + \sqrt{\varepsilon} V \left( \frac{\mathbf{x}}{\varepsilon} \right) \right) \quad \mathbf{x} \in \mathbb{R}^3,$$

where $V$ is a mean-zero random field modeling fluctuations around the background velocity $c_0$. A singular scattering kernel can then be obtained by considering a correlation function of the form

$$\mathbb{E}[V(\mathbf{x})V(\mathbf{y})] = R(\mathbf{x} - \mathbf{y}) = \int \frac{a(|\mathbf{p}|)}{|\mathbf{p}|^{2+\alpha}} e^{i\mathbf{p}\cdot(\mathbf{x}-\mathbf{y})} d\mathbf{p}. \quad (4.3)$$

However, such a formula holds only for $\alpha < 1$. The Fourier transform of $R$ being singular at $\mathbf{p} = 0$ (unbounded), the correlation function $R$ cannot be integrable, and more precisely we have $R(\mathbf{x}) \sim 1/|\mathbf{x}|^{1-\alpha}$. This situation corresponds to a random media with long-range correlations where a realization is provided in Fig. 4.1. Because of the long-range correlations one can observe significantly larger statistical patterns than for short-range correlations. In the context of neutronics and light scattering by biological tissues, the following Gegenbauer scattering kernel $\rho_G$ and Henyey-Greenstein (HG) kernel $\rho_{HG}$ are commonly used in the peaked-forward regime [114, 123, 178]:

$$F_G(s) := \frac{\alpha g (1 + g^2 - 2g s)^{-1-\alpha/2}}{2\pi((1-g)^{-\alpha} - (1+g)^{-\alpha})}, \quad F_{HG}(s) := \frac{1}{4\pi} \frac{1-g^2}{(1+g^2 - 2g s)^{3/2}}. \quad (4.4)$$

The parameter $g \in (-1, 1)$ is called the anisotropy factor. The case $g = 0$ corresponds to isotropic energy transfer over the unit sphere, $g < 0$ to dominant transfer in the backward direction, and $g > 0$ to forward energy transfer. The peaked forward regime is obtained in the limit $g \to 1$, for which

$$\frac{1}{(1-g)^\alpha} F_G(\hat{\mathbf{k}} \cdot \hat{\mathbf{p}}) \underset{g \to 1}{\sim} \frac{\alpha}{2\pi(2 - 2\hat{\mathbf{k}} \cdot \hat{\mathbf{p}})^{1+\alpha/2}} = \frac{\alpha}{2\pi|\hat{\mathbf{k}} - \hat{\mathbf{p}}|^{2+\alpha}}. \quad (4.5)$$

A typical realization of the corresponding random field in 2D, as $g \to 1$, is depicted in Fig. 4.1. Another example is provided by the so-called Kolmogorov power spectrum for standard atmospheric turbulence:

$$\Phi(|\mathbf{k}|) \propto \frac{a(|\mathbf{k}|)}{|\mathbf{k}|^{11/3}},$$

for $|\mathbf{k}|$ in the inertial range of turbulence, which corresponds to the case $\alpha = 5/3$.

In this chapter, the particular role played by the singularity of the scattering operator $Q$ is discussed at both the theoretical and numerical levels. The first part of this chapter presents regularization effects in all variables produced by such a singular kernel, and the diffusion limit of the RTE. The second part presents Monte-Carlo (MC) methods allowing to handle efficiently the singularity of the scattering operator.

This chapter covers the results of [8, 9, 13, 15] of the bibliography list on page 6.
## 4.1 Regularizing effects

In the spirit of [3, 40, 141] for the Boltzmann equation, our main motivation was to derive hypoelliptic type estimates that characterize the transfer of regularity from the momentum variable $\hat{\mathbf{k}}$ to the spatial variable $\mathbf{x}$ for solutions of (4.1). Similar problems are considered in the Euclidean case where $\mathbf{k} \in \mathbb{R}^d$ [4, 40, 92, 141] using various approaches. These results have been extended to RTEs where the momentum is confined to the unit sphere. In the case of the Boltzmann equation, the non-linearity in the collision operator makes the problem much harder than in the present linear situation [5, 6, 102, 103], but the non-integrability of the scattering kernel leads to comparable phenomena. The operator $\mathcal{L}$ roughly acts as a fractional Laplace-Beltrami operator on the unit sphere. A standard energy estimate then yields some Sobolev regularity in $\hat{\mathbf{k}}$, and the main question is to figure how this regularity is propagated to the spatial variable $\mathbf{x}$. Bootstrapping the estimates, the solution to (4.1) turns to be $\mathcal{L}^\infty$ in all its variables for any time $t > 0$, which is equivalent to the hypoellipticity of the operator $\partial_t + \hat{\mathbf{k}} \cdot \nabla_{\mathbf{x}} - \mathcal{Q}$. Compared to the Euclidean case, where $\mathbf{k} \in \mathbb{R}^d$, the fact that $\hat{\mathbf{k}} \in \mathbb{S}^{d-1}$ brings serious additional technical difficulties. To prove proper hypoelliptic estimates we follow the approach of [40] and build the estimates by regularization of the momentum variable using averaging lemmas on the sphere. These latter lemmas were established in [41–44] in various configurations, and needed to be extended to the case where the scattering operator in the transport equation involves fractional derivatives on the unit sphere. The hypoelliptic estimate we use to prove the global regularity for (4.1) is the following.

**Theorem 4.1.1** Assume $g, h \in L^2(\mathbb{R} \times \mathbb{R}^d \times \mathbb{S}^{d-1})$, and let $u \in L^2(\mathbb{R} \times \mathbb{R}^d \times \mathbb{S}^{d-1})$ satisfies the transport equation

$$ \partial_t u + \hat{\mathbf{k}} \cdot \nabla_{\mathbf{x}} u = (-\Delta_{\mathbb{S}^{d-1}})^\beta h + g \quad (4.6) $$

in the distribution sense, where $\Delta_{\mathbb{S}^{d-1}}$ is the Laplace-Beltrami operator on $\mathbb{S}^{d-1}$, and $\beta \ge 0$. For some $\theta > 0$, suppose in addition that

$$ (-\Delta_{\mathbb{S}^{d-1}})^{\theta/2} u \in L^2(\mathbb{R} \times \mathbb{R}^d \times \mathbb{S}^{d-1}). $$

Then, for

$$ \Upsilon := \frac{\theta}{2(1 + 2\beta) + \theta}, $$

we have $\partial_{t, \mathbf{x}}^\Upsilon u \in L^2(\mathbb{R} \times \mathbb{R}^d \times \mathbb{S}^{d-1})$ with the estimate

$$ \|\partial_{t, \mathbf{x}}^\Upsilon u\|_{L^2} \le C \left( \|(-\Delta_{\mathbb{S}^{d-1}})^{\theta/2} u\|_{L^2} + \|u\|_{L^2} + \|g\|_{L^2} + \|h\|_{L^2} \right), $$

for some constant $C > 0$, and where the $L^2$ norm holds on all the variables.

The above fractional derivatives are defined in the Fourier domain using the following convention:

$$ \hat{f}(\xi) = \mathcal{F} f(\xi) = \int e^{-i\mathbf{x} \cdot \xi} f(\mathbf{x}) d\mathbf{x}, \quad \mathcal{F}^{-1} f(\mathbf{x}) = \frac{1}{(2\pi)^d} \int e^{i\mathbf{x} \cdot \xi} \hat{f}(\xi) d\xi, $$

and introduce the fractional derivative as,

$$ \partial_{\mathbf{x}_j}^\Upsilon f(\mathbf{x}) = \mathcal{F}^{-1}[(i\xi_j)^\Upsilon \hat{f}(\xi)](\mathbf{x}) \quad \text{for } \Upsilon \in (0, 1), $$

with a similar definition for the fractional derivative involving the time variable and both variables. We denote by $\partial_{t, \mathbf{x}}^\Upsilon f$ any of the fractional derivatives w.r.t. $t$ and $\mathbf{x}_j$, and $\partial_{t, \mathbf{x}}^\Upsilon f = \partial_{t, \mathbf{x}}^{\Upsilon - [\Upsilon]} \partial_{t, \mathbf{x}}^{[\Upsilon]} f$ when $\Upsilon \ge 1$, where $[\Upsilon]$ is the integer part of $\Upsilon$.

From this result, replacing formally $(-\Delta_{\mathbb{S}^{d-1}})^\beta h$ in (4.6) by $Qu$, and using bootstrap arguments, we obtain that $u$ has derivatives of any order in all variables. The regularity result is the following.

**Theorem 4.1.2** The operator $\partial_t + \hat{\mathbf{k}} \cdot \nabla_{\mathbf{x}} - \mathcal{Q}$ is hypoelliptic. Namely, for any $u_0 \in L^2(\mathbb{R}^d \times \mathbb{S}^{d-1})$, (4.1) admits a unique classical solution $u$ that satisfies

$$ u \in \mathcal{C}^\infty((0, \infty) \times \mathbb{R}^d \times \mathbb{S}^{d-1}). $$
## 4.2 Diffusion limit

The diffusion limit described in this section is obtained through a probabilistic interpretation of the solution to (4.1). This probabilistic representation is also the basis of the numerical Monte-Carlo method proposed in the following section. According to [75, Thm. 7.2.1 pp. 302], after a proper extension of the operator $Q$, there exists a Markov process $K$ on $S^{d-1}$ with infinitesimal generator $Q$. Setting

$$D := (X, K), \quad \text{with} \quad X_t := X_0 - \int_0^t K(s)ds \quad t \ge 0,$$

defines a Markov process with infinitesimal generator

$$\mathcal{L}f(\mathbf{x}, \hat{\mathbf{k}}) = -\hat{\mathbf{k}} \cdot \nabla_{\mathbf{x}} f(\mathbf{x}, \hat{\mathbf{k}}) + Qf(\mathbf{x}, \hat{\mathbf{k}}),$$

from which we can write

$$u(t, \mathbf{x}, \hat{\mathbf{k}}) = \mathbb{E}_{\mathbf{x},\hat{\mathbf{k}}} [u_0(D(t))] := \mathbb{E}[u_0(D(t)) | D(0) = (\mathbf{x}, \hat{\mathbf{k}})].$$

The diffusion regime takes place in the presence of strong scattering effects, which can be represented using the following changes of variables:

$$t \to \frac{t}{\varepsilon^2} \quad \text{and} \quad \mathbf{x} \to \frac{\mathbf{x}}{\varepsilon},$$

and considering

$$u^{\varepsilon}(t, \mathbf{x}, \hat{\mathbf{k}}) = u\left(\frac{t}{\varepsilon^2}, \frac{\mathbf{x}}{\varepsilon}, \hat{\mathbf{k}}\right),$$

with $u^{\varepsilon}(t = 0, \mathbf{x}, \mathbf{k}) = u_0(\mathbf{x}, \mathbf{k})$. This scaling corresponds to a long-time and space scaling providing at the macroscopic scale strong scattering effects. The transport equation (4.1) then becomes

$$\varepsilon^2 \partial_t u^{\varepsilon} + \varepsilon \hat{\mathbf{k}} \cdot \nabla_{\mathbf{x}} u^{\varepsilon} = Q u^{\varepsilon}, \quad (4.7)$$

and the associated probabilistic representation is obtained by considering the rescaled Markov process

$$D^{\varepsilon}(t) := \left( X^{\varepsilon}(0) - \varepsilon \int_{0}^{t/\varepsilon^2} K(s)ds, K(t/\varepsilon^2) \right),$$

with generator

$$\mathcal{L}_{\varepsilon}f = -\frac{1}{\varepsilon}\hat{\mathbf{k}} \cdot \nabla_{\mathbf{x}}f + \frac{1}{\varepsilon^2}Qf.$$

The solution to (4.7) then reads

$$u^{\varepsilon}(t, \mathbf{x}, \hat{\mathbf{k}}) = \mathbb{E}_{\mathbf{x},\hat{\mathbf{k}}} [u_0(D^{\varepsilon}(t))],$$

and its asymptotic behavior in $\varepsilon$ is given by the following result.

**Theorem 4.2.1** Suppose there exists $\eta > 0$ such that $a \ge \eta$ almost everywhere, and let

$$Y_{\mathbf{y}}^{\varepsilon}(t) = \mathbf{y} - \varepsilon \int_{0}^{t/\varepsilon^2} K(s)ds.$$

Then, the process $Y_{\mathbf{y}}^{\varepsilon}$ converges in law in $C^0([0, \infty), \mathbb{R}^d)$, as $\varepsilon \to 0$, to a diffusion process $Y_{\mathbf{y}}$, starting at $\mathbf{y}$, with generator

$$\mathcal{L}_0 = \nabla_{\mathbf{x}} \cdot A \nabla_{\mathbf{x}}.$$

The positive-definite diffusion matrix $A$ is given by

$$A_{jl} = \frac{1}{\mathcal{C}} \int_{\mathbb{S}^{d-1}} \tilde{\sigma}(d\hat{\mathbf{k}}) \hat{k}_j \hat{k}_l \quad \text{with} \quad \mathcal{C} = \sigma(\mathbb{S}^{d-2}) \int_{-1}^{1} F(s)(1-s^2)^{(d-3)/2} (1-s)ds.$$
Above, $\tilde{\sigma}$ is the uniform measure on $\mathbb{S}^{d-1}$ (i.e. $\tilde{\sigma} = \sigma/\sigma(\mathbb{S}^{d-1})$). Moreover, for any $u_0 \in L^2(\mathbb{R}^d \times \mathbb{S}^{d-1})$ and $t > 0$, $u^\varepsilon(t)$ converges weakly in $L^2(\mathbb{R}^d \times \mathbb{S}^{d-1})$, as $\varepsilon \to 0$, to the unique solution to

$$ \partial_t \tilde{u} = \nabla_x \cdot (A \nabla_x \tilde{u}), \quad \text{with} \quad \tilde{u}(t=0, \mathbf{x}) = \int_{\mathbb{S}^{d-1}} u_0(\mathbf{x}, \mathbf{p}) d\tilde{\sigma}(\mathbf{p}), $$

and where the function $\tilde{u}$ reads

$$ \tilde{u}(t, \mathbf{x}) = \mathbb{E}[\tilde{u}(t=0, Y_{\mathbf{x}}(t))]. $$

The main ingredient of the proof is a spectral gap estimate showing that the Markov process $Y^\varepsilon$ is ergodic, and allows the use of standard approximation-diffusion theorems to obtain the convergence. The hypothesis that $a$ is strictly positive is crucial to obtain this spectral gap. Moreover, the constant $\mathcal{C}$ is non-zero and finite, and so does $A$, despite the long-range correlations and the singular behavior of $F$.

## 4.3 A Monte-Carlo methods for radiative transfer with singular kernels in 2D propagation media

MC methods have been used for a long time to solve numerically (4.1) (see [140,192]). These methods have several advantages compared to classical finite element or finite volume methods. They can handle the high-dimensionality in a simple manner, since no mesh is required and the quantities of interest only need be discretized at a detector where measurements take place. They are very flexible in terms of geometry, translate boundary conditions onto the underlying stochastic process, and are easily applicable even when the coefficients in the equation depend on the variables $(\mathbf{x}, \hat{\mathbf{k}})$. MC methods are also very easy to implement and can be parallelized in a straightforward manner, which is a strong feature with today's technology. Note that discretization-based methods would have to handle carefully the singularity of the cross-section. Without singular scattering operators, the standard MC methods rely on the mean free time between two scattering events. However, as shown by (4.2), when considering a singular scattering operator this mean free time is zero. This section illustrates how to overcome efficiently this difficulty.

With a singular scattering operator $Q$, the underlying Markov process $K$ of the MC method no longer relies on a Poisson process to describe its jump times as in the case of a nonnull mean free time. In our context, $K$ corresponds to a general jump process with infinite jump intensity on $\mathbb{S}^1$. The simulation of such Lévy processes in the context of transport equations seems much less studied than its compound Poisson process counterpart. We focus in this section on the two-dimensional case with momentum on the circle $\mathbb{S}^1$, and will address the three-dimensional case in the next section. Note that the choice of $\mathbb{S}^1$ is arbitrary, any circle of given radius would be handle similarly.

The key for efficient MC methods resides in an appropriate simulation method of the process $K$ representing the scattering mechanism w.r.t. the momentum variable. To this end, we adapt a method of the probabilistic literature developed by Asmussen-Rosiński [9], and Cohen-Rosiński [53], that we refer to as the ACR method. The first ingredient is to introduce a cut-off in the scattering operator to remove the singularity. The resulting operator is then the one of a compound Poisson process with finite (but large) intensity. The main issue with this cut-off approach is that a small value of the cut-off parameter is required for good accuracy, leading to a high intensity and therefore to many scattering events increasing the computational time. This issue is fixed by the second ingredient of this method, which consists in adding a correction term describing asymptotically the behavior of the operator around the singularity. This term models frequent jumps of small amplitude, while the regularized generator describes more rare jumps with larger amplitudes. In the work of Asmussen-Rosiński, and Cohen-Rosiński, the correction takes the form of a Brownian motion on the real line, while in the present context it corresponds to a Brownian motion on the circle (the unit sphere in the 3D case).

For comparison, and to illustrate the efficiency of the ACR method for the simulation of RTEs, an alternative MC method that does not involve a regularization of the generator is also introduced. This alternative method is based on the so-called Alpha-Stable processes, and is referred to as the AS method.
### 4.3.1 The ACR method

To introduce the ACR method, the operator $\mathcal{Q}$ needs to be appropriately parametrized. The natural choice on $\mathbb{S}^1$ are the polar coordinates and we introduce the surjective map

$$P: \phi \in \mathbb{R} \mapsto (\cos(\phi), \sin(\phi)) \in \mathbb{S}^1.$$

Setting $\hat{\mathbf{k}} = (\cos(\theta), \sin(\theta))$ and $\hat{\mathbf{p}} = (\cos(\theta'), \sin(\theta'))$, so that $\hat{\mathbf{k}} \cdot \hat{\mathbf{p}} = \cos(\theta' - \theta)$, the scattering kernel $\mathcal{Q}$ can be recast as

$$ \begin{aligned} \mathcal{Q}(f)(k) &= \frac{1}{2^{(1+\alpha)/2}} \int_{-\pi}^{\pi} \frac{a(\theta' - \theta)(\tilde{f}(\theta') - \tilde{f}(\theta))}{(1 - \cos(\theta' - \theta))^{(1+\alpha)/2}} d\theta' \\ &= \frac{1}{2^{(1+\alpha)/2}} \int_{-\pi}^{\pi} \frac{a(\theta')(\tilde{f}(\theta' + \theta) - \tilde{f}(\theta))}{(1 - \cos(\theta'))^{(1+\alpha)/2}} d\theta' \\ &= \tilde{\mathcal{Q}}(\tilde{f})(\theta), \quad \theta \in [-\pi, \pi], \end{aligned} $$

where $\tilde{f} = f \circ P$ is a $2\pi$-periodic function. Note that the $2\pi$-periodic function $a$ is symmetric with respect to 0 and is nonnegative. The operator $\tilde{\mathcal{Q}}$ can be seen as the generator of a Lévy process on $\mathbb{R}$ starting from 0, that we denote by $(\Theta(t))_{t \ge 0}$. According to our parametrization,

$$K = P(\Theta + \theta)$$

in law, when $K(0) = \hat{\mathbf{k}}$. For $\varepsilon \in (0, 1)$, the ACR method consists in splitting $\tilde{\mathcal{Q}}$ into

$$ \begin{aligned} \tilde{\mathcal{Q}}(g)(\theta) &= \mathcal{Q}_{<}^{\varepsilon}(g)(\theta) + \mathcal{Q}_{>}^{\varepsilon}(g)(\theta) \\ &:= \frac{1}{2^{(1+\alpha)/2}} \int_{I_{\varepsilon}^{C}} \frac{a(\theta')(g(\theta' + \theta) - g(\theta))}{(1 - \cos(\theta'))^{(1+\alpha)/2}} d\theta' + \frac{1}{2^{(1+\alpha)/2}} \int_{I_{\varepsilon}} \frac{a(\theta')(g(\theta' + \theta) - g(\theta))}{(1 - \cos(\theta'))^{(1+\alpha)/2}} d\theta' \end{aligned} $$

where $I_{\varepsilon}$ is defined by

$$I_{\varepsilon} := \{\theta' \in [-\pi, \pi], \ |\tan(\theta'/4)| > \varepsilon/4\},$$

and $I_{\varepsilon}^C$ its complementary in $[-\pi, \pi]$. Because of the regularization, the second part of the generator, namely $\mathcal{Q}_{>}^{\varepsilon}(g)$, is now the generator of a classical compound Poisson process $\theta_{\varepsilon}$ with jump intensity

$$\mu_{\varepsilon} = \frac{\Pi_{\varepsilon}^{0}}{2^{(1+\alpha)/2}}, \quad \text{where} \quad \Pi_{\varepsilon}^{0} = \int_{I_{\varepsilon}} \frac{a(\theta)d\theta}{(1 - \cos(\theta))^{(1+\alpha)/2}},$$

and jump distribution

$$\Pi_{\varepsilon}(d\theta) = \frac{\mathbf{1}_{I_{\varepsilon}}(\theta)a(\theta)}{\Pi_{\varepsilon}^{0}(1 - \cos(\theta))^{(1+\alpha)/2}} d\theta.$$

The expected number of jumps for the compound Poisson process in any time interval $[0, T]$ is given by $\mu_{\varepsilon}T$, but the term $\mu_{\varepsilon}$ blows up as $\varepsilon \to 0$,

$$\mu_{\varepsilon} \propto \frac{1}{\alpha \varepsilon^{\alpha}}.$$

In other words the process jumps more and more as $\varepsilon$ becomes small, this rate blows up even more rapidly when $\alpha$ gets close to 2, and hence increases the simulation cost. An efficient simulation technique for the jump distribution is therefore required. Good performances to approximately simulate numerically the jump distribution is achieved using the stochastic collocation method [104].

The first term of the generator $\mathcal{Q}_{<}^{\varepsilon}$ models frequent jumps of small amplitude and can be asymptotically reduced to the generator of a Brownian motion on the circle. The following lemma describes the relation between $\mathcal{Q}_{<}^{\varepsilon}$ and the Laplace operator.
**Lemma 4.3.1** Suppose that $a \in C^2([-\pi, \pi])$ and that $g$ is a smooth $2\pi$-periodic bounded function with bounded derivatives. Then, for any $\phi \in [-\pi, \pi]$,

$$Q_{<}^{\epsilon}(g)(\phi) = \frac{\sigma_{\epsilon}}{2} \frac{\partial^2 g(\phi)}{\partial \phi^2} + \mathcal{R}_{\epsilon}[g](\phi),$$

where

$$\sigma_{\epsilon} = \frac{2a(0)\epsilon^{2-\alpha}}{2-\alpha} \quad \text{and} \quad \|\mathcal{R}_{\epsilon}[g]\|_{\infty} \le \frac{\epsilon^{4-\alpha}}{4-\alpha} \|a\|_{C^2([-\pi, \pi])} \|g\|_{C^4([-\pi, \pi])}.$$

The leading term in $Q_{<}^{\epsilon}$ is the generator of $\phi_{\epsilon} = (\sigma_{\epsilon}W_t \mod 2\pi)_t$, where $W$ is a Brownian motion. Brownian trajectories can easily be generated, but note that $\sigma_{\epsilon}$ is proportional to $\epsilon^{2-\alpha}/(2-\alpha)$, which becomes larger as $\alpha$ increases. We therefore expect the inclusion of the Brownian term in the simulation to be critical in order to obtain good accuracy and computational cost.

## Error analysis

The main idea of the ACR method is that the Lévy process $K$ on $\mathbb{S}^1$, starting at $\hat{\mathbf{k}}$ with generator $Q$, is approximated by the process

$$K^{\epsilon} = (\cos(\phi_{\epsilon} + \theta_{\epsilon} + \theta), \sin(\phi_{\epsilon} + \theta_{\epsilon} + \theta)),$$

which is the one simulated numerically in practice. Here $\theta_{\epsilon}$ accounts for the regularized compound Poisson process and $\phi_{\epsilon}$ for the Brownian correction. Note that these two processes are independent. Introducing

$$X_{\mathbf{x}}^{\epsilon}(t) = \mathbf{x} - \int_{0}^{t} K^{\epsilon}(s)ds,$$

the function

$$u_{\epsilon}(t, \mathbf{x}, \hat{\mathbf{k}}) = \mathbb{E}[u_0(X_{\mathbf{x}}^{\epsilon}(t), K^{\epsilon}(t)) | K^{\epsilon}(0) = \hat{\mathbf{k}}] \quad (4.8)$$

is then the solution of the following transport equation, with the parametrization $\hat{\mathbf{k}} = (\cos(\theta), \sin(\theta))$:

$$\partial_t u_{\epsilon} + \hat{\mathbf{k}} \cdot \nabla_{\mathbf{x}} u_{\epsilon} = Q_{>}^{\epsilon}(u_{\epsilon}) + \frac{D_{\epsilon}}{2} \partial_{\theta}^{2} u_{\epsilon}, \quad u_{\epsilon}(t=0) = u_0.$$

The following proposition quantifies the approximation error, in term of the cut-off parameter, when using this model to approximate the original RTE.

**Proposition 4.3.1** Suppose that $a \in C^2([-\pi, \pi])$. Then, for all $t \in (0, T)$,

$$\|u_{\epsilon}(t, \cdot, \cdot) - u(t, \cdot, \cdot)\|_{L^2(\mathbb{R}^2 \times \mathbb{S}^1)} \le \frac{\epsilon^{4-\alpha}}{4-\alpha} \|a\|_{C^2([-\pi, \pi])} \int_0^t \|u(s)\|_{L^2(\mathbb{R}^2, C^4(\mathbb{S}^1))} ds.$$

Recall that $u$ satisfies the regularity requirement thanks to the hypoelliptic property of the RTE.

## 4.3.2 The AS method

This section introduces an alternative method to compare against the ACR method. The main idea is to exactly simulate the process $K$, without using the small jumps - larger jumps approximation of the generator. The method is applicable provided an additional assumption on the function $a$. In this section $a(\theta)$ is assumed to admit a global minimizer at $\theta = 0$, or equivalently $A(s) = a(\theta)$ (for $s = \hat{\mathbf{k}} \cdot \hat{\mathbf{p}} = \cos(\theta)$) admits a global minimizer at $s = 1$. This assumption allows us to decompose $Q$ in terms of two infinitesimal generators

$$\begin{aligned} Q(f)(\hat{\mathbf{k}}) &= \frac{A(1)}{2^{(1+\alpha)/2}} \int_{\mathbb{S}^1} \frac{f(\hat{\mathbf{p}}) - f(\hat{\mathbf{k}})}{(1 - \cos(\hat{\mathbf{k}} \cdot \hat{\mathbf{p}}))^{(1+\alpha)/2}} \sigma(d\hat{\mathbf{p}}) \\ &+ \frac{1}{2^{(1+\alpha)/2}} \int_{\mathbb{S}^1} \frac{A(\hat{\mathbf{k}} \cdot \hat{\mathbf{p}}) - A(1)}{(1 - \cos(\hat{\mathbf{k}} \cdot \hat{\mathbf{p}}))^{(1+\alpha)/2}} (f(\hat{\mathbf{p}}) - f(\hat{\mathbf{k}})) \sigma(d\hat{\mathbf{p}}) \\ &=: Q^1(f)(\hat{\mathbf{k}}) + Q^2(f)(\hat{\mathbf{k}}). \end{aligned}$$
The second term $Q^2$ is the generator of a compound Poisson process on $\mathbb{S}^1$ with a finite jump intensity, which does not pose any computational issue since now the average length of the jumps is large compared to those generated by $Q^\varepsilon$. The main difficulty is to simulate numerically random trajectories for a Lévy process $K^1$ with self-adjoint generator $Q^1$. Such trajectories are simulated numerically through piecewise constant trajectories on a fine time grid, and each transitions are simulated using the transition density function

$$p(h, \hat{\mathbf{k}}, \hat{\mathbf{q}}) = \mathbb{P}(K^1(h) = \hat{\mathbf{k}} | K^1(0) = \hat{\mathbf{q}}),$$

where $h$ corresponds to the time stepsize. This function can be derived explicitly by solving the forward Kolmogorov equation

$$\partial_t p(\hat{\mathbf{k}}) = Q^1(p)(\hat{\mathbf{k}}), \quad \text{with} \quad p(t=0, \hat{\mathbf{k}}, \hat{\mathbf{q}}) = \delta(\hat{\mathbf{k}} - \hat{\mathbf{q}}),$$

for $\hat{\mathbf{k}}, \hat{\mathbf{q}} \in \mathbb{S}^1$, and where $\delta$ is the Dirac measure. With the parametrization $\hat{\mathbf{k}} \cdot \hat{\mathbf{q}} = \cos(\theta)$, the operator $Q$ can be diagonalized using Fourier series, and we obtain the following exact expression for $p$ (written as a function of $\theta$),

$$p(t, \theta) = \frac{1}{2\pi} + \frac{1}{\pi} \sum_{\ell=1}^{\infty} e^{t\lambda_{\ell}} \cos(\ell\theta),$$

where (with $a_0 = a(0) = A(1)$)

$$\lambda_{\ell} = \frac{a_0 \pi^{\frac{1}{2}} \Gamma(-\frac{\alpha}{2})}{2^{\alpha} \Gamma(\frac{1+\alpha}{2})} \left( \frac{\Gamma(\ell + \frac{1+\alpha}{2})}{\Gamma(\ell + \frac{1-\alpha}{2})} - \frac{\Gamma(\frac{1+\alpha}{2})}{\Gamma(\frac{1-\alpha}{2})} \right),$$

according to [182] with a slight adjustment of the constant prefactors, and $\Gamma$ the usual Gamma function. It does not seem possible to directly simulate the increments using $p$. To address this issue, we should use an approximation of it that enables fast simulations. Once we have this approximation, we can once again apply the stochastic collocation method [104] to provide accurate simulations to $p$. The operator $Q^1$ should not be far from a fractional Laplace-Beltrami operator on $\mathbb{S}^1$. Therefore, we introduce

$$p_0(h, \theta) := \frac{1}{2\pi} + \frac{1}{\pi} \sum_{\ell=1}^{\infty} e^{-hD_{\alpha}\ell^{\alpha}} \cos(\ell\theta)$$

as a rough approximation to $p(h, \cdot)$, which corresponds to the distribution to $(hD_\alpha)^{1/\alpha} X \mod 2\pi$, where $X$ stands for a random variable with an $\alpha$-stable distribution, which has characteristic function

$$\mathbb{E}[e^{ivX}] = e^{-|v|^\alpha}.$$

Following the strategy of [208], the $\alpha$-stable distribution can be easily exactly simulated numerically from the uniform distribution over $(0, 1)$.

It is worth mentioning that the density $p_0$ is already an excellent approximation of $p$, as shown in Fig. 4.2. For both cases illustrated in this figure, the $L^2$ relative error is about 0.3%. Even though these small errors translate into small errors in the simulation of the RTE, it is not possible to obtain an arbitrary precision by using $p_0$ instead of $p$. For arbitrary accuracy, a random variable drawn according to $p_0$ needs to be corrected to be seen as distributed according to $p$. This can be achieved with the use of the efficient stochastic collocation method [104].

### 4.3.3 Numerical Illustrations

To perform relevant numerical simulations, we need first to define a typical time scale. When the scattering kernel is integrable this time scale is given by the mean free time (4.2). In our non-integrable case, this mean free time is zero and we need a different quantity. A natural one is the inverse of $-\lambda_1$, that is the inverse of the first non zero eigenvalue of $-Q$. When $t \gg t_S := (-\lambda_1)^{-1}$, the distribution of the angles is almost uniform over $\mathbb{S}^1$. Fig. 4.3 represents the characteristic time $t_S$ as a function of $\alpha$ for $a \equiv a_0 = 0.1$.
<figure><img src="image_9.png" /><figcaption>Fig. 4.2: Comparison between $p$ (blue) and $p_0$ (green), $a_0 = 0.1$ and $h = 0.03$, for $\alpha = 0.8$ (left) and $\alpha = 1.6$ (right).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 4.3: Characteristic time $t_S$ as a function of $\alpha$ for $a(s) = a_0 = 0.1$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 4.4: Comparison between the reference solution and the numerical solutions including or not the Brownian correction. The non corrected solution is completely off when $\alpha = 1.9$, while there is a good accuracy when the Brownian part is included. Here, $h = 0.03$ for the time stepsize, $\epsilon = a_0 = 0.1$, $t = 3$, $N_B = 80$, $z_1 = 2.7$, $z_2 = 3$, and $10^7$ trajectories have been used for the MC method.</figcaption></figure>

## Impact of the Brownian correction in the ACR method

In view of Lem. 4.3.1 this correction is expected to play a crucial role for large $\alpha$. To highlight the role of this correction, we compare the numerical simulations against a reference solution. In our context no explicit analytical formula for the solution of the RTE is available, so that we only rely on a semi-analytical expression corresponding to

$$J_i(t) = \int_{z_1}^{z_2} dx_1 \int_{\mathbb{R}} dx_2 \int_{\theta_i}^{\theta_{i+1}} d\theta \, u(t, x_1, x_2, \theta)$$

with $\mathbf{x} = (x_1, x_2)$, and $\hat{\mathbf{k}} = (\cos(\theta), \sin(\theta))$. Here, $u$ is the solution to the RTE with a particular initial condition

$$u_0(\mathbf{x}, \theta) = \frac{1}{\sqrt{2\pi c^2}} e^{-\frac{|\mathbf{x}|^2}{2c^2}} \frac{1}{2\pi} (1 + \cos(\theta)),$$

with $c = 0.1$ for the numerical evaluation. The time is here fixed at $t = 3$, which is roughly the characteristic time $t_S$ depicted in Fig. 4.3. When $\alpha = 1.9$, the observation time is about $3t_S$. The solution $u$ is integrated over several bins $[\theta_i, \theta_{i+1}]$, $i = 1, \dots, N_B$, $\theta_{i+1} = \theta_i + \Delta\theta$. In Fig. 4.4, we represent the reference solution $(J_i(t))_{i=1,\dots,N_B}$ and the numerical solution to the RTE $u_{\epsilon,h}$, including or not the Brownian correction. When $\alpha = 1.9$ the solution is completely off when the correction is not included, while the relative error is about 3% when it is included. When $\alpha = 1.1$, the relative errors are 6% without the Brownian correction, and 0.3% with it, and 1.4% and 0.4% respectively for $\alpha = 0.3$. In all cases, this shows that the Brownian correction is crucial in order to obtain a good accuracy in the simulations.

## Comparison ACR/AS methods

The ACR and AS methods are compared to each other in terms of accuracy versus computational cost for $\alpha = 0.7, 1.3, 1.8$. For the evaluation of the accuracy, the outcomes of both methods are compared to our reference solution $(J_i(t))_{i=1,\dots,N_B}$ with $N_B = 20$. All the other parameters remain as before. The errors are averaged over 20 simulations of the MC methods to obtain more stable estimates. The results are depicted in Fig. 4.5. The cost is expressed in units of 0.17s, which is the computational time to achieve an error of about 10% for a code written in the Julia language, and ran in parallel on a 32 cores Intel Xeon E5-2697A at 2.60GHz.

In the case $\alpha = 0.7$, the cost to achieve an error of 0.2% is about 10 times higher for the AS method than for the ACR method. This difference is reduced for lower accuracies. Similar observations hold for $\alpha = 1.3$ and $\alpha = 1.8$, with the gap shrinking as $\alpha$ increases. This can be explained on the one hand by the fact that the ACR method is more sensitive to the parameter $\alpha$ than the AS method. In fact, there are more and more scattering events as $\alpha$ grows, and on the other hand the jumps are smaller
<figure><img src="image_12.png" /><figcaption>Fig. 4.5: Comparison of the AS and ACR methods in terms of cost versus accuracy for $\alpha = 0.7, 1.3, 1.8$. The ACR method is consistently more efficient.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 4.6: Realizations of Gaussian random fields. In the left picture, we have $\lambda = 1_{\{|x_1|<15\}}$, and $\alpha = 0.1 \cdot 1_{\{|x|\le 10\}} + 1 \cdot 1_{\{10<|x|\}}$. In the right picture, we have $\lambda = 1_{\{x_2\in(0,20)\}}$ and $\alpha(x_2) = 5/3 \cdot 1_{\{x_2\le 2\}} + 0.5 \cdot 1_{\{2<x_2\le 8\}} + 1.9 \cdot 1_{\{8<x_2\}}$.</figcaption></figure>

for large $\alpha$, which reduces the error when evaluating the integral for the position variable $X(t)$ in the AS method. Nevertheless, the cost of a 0.5% error when $\alpha = 1.8$ is more than twice higher for the AS method, and the cost significantly increases when going from an accuracy of about 10% to about 0.1%.

## 4.4 A Monte-Carlo method for 3D radiative transfer equations with multifractional singular kernels

The radiative transfer model (4.1) can be extended to more general situations by letting the properties of the scattering operator depends on the position $x$:

$$ \Phi(\mathbf{x}, |\hat{\mathbf{p}} - \hat{\mathbf{k}}|) := \frac{\lambda(\mathbf{x})a(|\hat{\mathbf{p}} - \hat{\mathbf{k}}|)}{|\hat{\mathbf{p}} - \hat{\mathbf{k}}|^{2+\alpha(\mathbf{x})}} = \frac{\lambda(\mathbf{x})}{2^{1+\alpha(\mathbf{x})/2}} \rho(\mathbf{x}, \hat{\mathbf{k}} \cdot \hat{\mathbf{p}}), $$

with

$$ \rho(\mathbf{x}, s) := \frac{a(\sqrt{2(1-s)})}{(1-s)^{1+\alpha(\mathbf{x})/2}} \quad s \in [-1, 1]. \qquad (4.9) $$

Above, $\alpha: \mathbb{R}^3 \to (0,2)$ accounts for variations of scattering properties across the ambient space, and $\lambda \ge 0$ is a function modeling the support of the scattering region. Regions where $\lambda(\mathbf{x}) = 0$ are homogeneous and $u$ undergoes a free transport. In Fig. 4.6, we provide examples of 2D random fields for which $\alpha$ and $\lambda$ are not constant. In these two pictures, we highlight the roles of $\lambda$ and $\alpha$: $\lambda$ characterizes scattering regions, and $\alpha$ defines the correlation structure. In the inner circle of
the left picture we have $\alpha \equiv 0.1$, which tends to create shorter correlation patterns than outside the circle where $\alpha \equiv 1$. In the right picture, we have a three-layer model for $\alpha$ in which the inner band exhibits smaller statistical patterns than the outer ones. This type of model is used for modeling non-Kolmogorov atmospheric turbulences.

From a numerical perspective, the 3D context outline in this section is more involved compared to the 2D setting, where the small-jumps part can simply be approximated by a Brownian motion mod $2\pi$. Here, the singularity is approximated using a Laplace-Beltrami operator (with respect to the momentum variable) on the unit sphere $S^2$, so that simulating paths of a jump-diffusion process over the unit sphere is necessary to handle the small-jumps part. Remember that disregarding small jumps altogether in order to use standard MC methods results in substantial errors, and mitigating these errors comes with a notable increase in computational cost.

### 4.4.1 Approximation and probabilistic representation

In order to adapt the ACR method, we introduce the following small region over which the singularity of the kernel $\rho$ (in (4.9)) is not integrable:

$$S_{<}^{\epsilon} = S_{<}^{\epsilon}(\hat{\mathbf{k}}) := \{\hat{\mathbf{p}} \in S^2 : 1 - \hat{\mathbf{p}} \cdot \hat{\mathbf{k}} < \epsilon\} \quad \text{for } \epsilon \in (0, 1). \quad (4.10)$$

The ACR method is based on decomposing the scattering operator $Q$ into two parts:

$$\begin{aligned} \mathcal{L}Q(\mathbf{x}, \hat{\mathbf{k}}) &= +Q_{<}^{\epsilon}f(\mathbf{x}, \hat{\mathbf{k}}) + Q_{>}^{\epsilon}f(\mathbf{x}, \hat{\mathbf{k}}) \\ &:= \frac{\lambda(\mathbf{x})}{2^{1+\alpha(\mathbf{x})/2}} \left( \int_{S_{<}^{\epsilon}} + \int_{S_{>}^{\epsilon}} \right) \rho(\mathbf{x}, \hat{\mathbf{p}} \cdot \hat{\mathbf{k}}) (f(\mathbf{x}, \hat{\mathbf{p}}) - f(\mathbf{x}, \hat{\mathbf{k}})) \sigma(d\hat{\mathbf{p}}), \end{aligned}$$

where $S_{>}^{\epsilon} = (S_{<}^{\epsilon})^c$ is the complementary set of the region (4.10) over the unit sphere. The term $Q_{>}^{\epsilon}$ of the scattering operator, with no singularity, is the infinitesimal generator of a standard jump Markov process with large-jumps. Regarding $Q_{<}^{\epsilon}$, with the singularity, the following result justifies the approximation of this operator by a Laplace-Beltrami operator $\Delta_{S^2}$ over the unit sphere $S^2$. In what follows, we use the notations $r'_{\epsilon} = \sqrt{1 - (1 - \epsilon)^2/(2 - \epsilon)}$ and $0 \le \alpha_m \le \alpha(\mathbf{x}) \le \alpha_M < 2$.

**Proposition 4.4.1** Let $u$ be the solution to (4.1) with kernel (4.9), and $u_{\epsilon}$ be the solution of

$$\left\{ \begin{array}{l} \partial_t u_{\epsilon} + \hat{\mathbf{k}} \cdot \nabla_{\mathbf{x}} u_{\epsilon} = \sigma_{\epsilon}^2(\mathbf{x}) \Delta_{S^2} u_{\epsilon} + \frac{\lambda(\mathbf{x})}{2^{1+\alpha(\mathbf{x})/2}} \int_{S_{>}^{\epsilon}} \rho(\mathbf{x}, \hat{\mathbf{p}} \cdot \hat{\mathbf{k}}) (u_{\epsilon}(\hat{\mathbf{p}}) - u_{\epsilon}(\hat{\mathbf{k}})) \sigma(d\hat{\mathbf{p}}), \\ u_{\epsilon}(0, \mathbf{x}, \hat{\mathbf{k}}) = u_0(\mathbf{x}, \hat{\mathbf{k}}), \end{array} \right. \quad (4.11)$$

for $(t, \mathbf{x}, \hat{\mathbf{k}}) \in (0, \infty) \times \mathbb{R}^3 \times S^2$, where

$$\sigma_{\epsilon}^{2}(\mathbf{x}) := \frac{2^{1-\alpha(\mathbf{x})}a(0)\pi\lambda(\mathbf{x})}{2-\alpha(\mathbf{x})} r'_{\epsilon}^{2-\alpha(\mathbf{x})}. \quad (4.12)$$

Assuming $a'(0) = 0$, for any $T > 0$, we have

$$\sup_{t \in [0, T]} \|u(t, \cdot, \cdot) - u_{\epsilon}(t, \cdot, \cdot)\|_{L^2(\mathbb{R}^3 \times S^2)} \le \varepsilon^{2-(\alpha_M/2)} \sqrt{2T} E(u), \quad (4.13)$$

where $E(u)$ is a constant depending on the 2nd- and 4th-order derivatives of $u$.

As we know from Sect. 4.1, when $\lambda$ and $\alpha$ are constants, $u$ is infinitely differentiable in all variables for any square integrable initial conditions. However, when $\lambda$ and $\alpha$ are infinitely differentiable with bounded derivatives at all orders (from above and below), this result remains valid so $u$ is still assumed to be smooth.

Based on the approximation result (4.13), we then devise a MC method for (4.11) instead of (4.1). The advantage in using (4.11) lies in that the angular diffusion term $\sigma_{\epsilon}^2(\mathbf{x})\Delta_{S^2}$ is the generator of a Markov process that can be easily simulated. Indeed, for $W$ a standard 3D Brownian motion on $\mathbb{R}^3$, and $\times$ being the cross product in $\mathbb{R}^3$, it is shown in [205] that the process $B$ solving the SDE

$$dB = B \times dW - Bdt, \quad B(0) \in S^2,$$
has infinitesimal generator $\frac{1}{2}\Delta_{\mathbb{S}^2}$, and corresponds to a Brownian motion over $\mathbb{S}^2$. A simple adaptation gives the desired diffusion coefficient. Note also that $\sigma_{\varepsilon}(\mathbf{x})$ increases as $\alpha(\mathbf{x})$ gets close to 2, and diffusion on the sphere eventually becomes the dominant dynamics.

In this section the MC method is built by interpreting (4.11) as the forward Kolmogorov equation of an appropriate Markov process, and as a consequence focusing for simplicity on the forward MC method (see [140] for terminology). The Markov process we consider for this approach is defined by

$$D_{\varepsilon}(t) := \sum_{n \ge 0} \mathbf{1}_{[T_n, T_{n+1})}(t) \psi_n^{Z_n}(t-T_n) \quad t \ge 0, \quad (4.14)$$

where for $z = (\mathbf{x}, \hat{\mathbf{k}})$:

1. The flow $\psi_n^z = (X_n^x, K_n^{\hat{\mathbf{k}}})$ is the unique strong solution to the SDE

$$\left\{
\begin{array}{l}
dX_n^x(t) = K_n^{\hat{\mathbf{k}}}(t) dt \\
dK_n^{\hat{\mathbf{k}}}(t) = \sqrt{2} \sigma_{\varepsilon}(X_n^x(t)) K_n^{\hat{\mathbf{k}}}(t) \times dW_n(t) - 2 \sigma_{\varepsilon}^2(X_n^x(t)) K_n^{\hat{\mathbf{k}}}(t) dt,
\end{array}
\right.
\quad (4.15)$$

where $\psi_n^z(0) = z$, $(W_n)_n$ is a sequence of independent standard Brownian motions on $\mathbb{R}^3$, and $\sigma_{\varepsilon}$ is defined by (4.12).

2. The jump times $(T_n)_n$ are distributed according to

$$\mathbb{P}(T_{n+1} - T_n > t | D_{\varepsilon}(T_n) = z, (\psi_n^z(s))_{s \in [0, t]}) = \exp \left( - \int_0^t \Lambda_{\varepsilon}(\psi_n^z(s)) ds \right) \quad n \ge 0,$$

with $T_0 = 0$, and for $\rho$ given by (4.9),

$$\Lambda_{\varepsilon}(z) := \frac{\lambda(\mathbf{x})}{2^{1+\alpha(\mathbf{x})/2}} \int_{S_{\varepsilon}^{z}} \rho(\mathbf{x}, \hat{\mathbf{p}} \cdot \hat{\mathbf{k}}) \sigma(d\hat{\mathbf{p}}).$$

3. The jumps $(Z_n)_n$ describe a Markov chain with transition probability

$$\mathbb{P}(Z_{n+1} \in dy \otimes \sigma(d\hat{\mathbf{p}}) | Z_n, T_{n+1} - T_n) = \Pi_{\varepsilon}(z_{n+1}, dz),$$

where $z_{n+1} := \psi_n^{Z_n}(T_{n+1} - T_n)$, and

$$\Pi_{\varepsilon}(z, dz) := \pi_{\varepsilon}(z, \hat{\mathbf{p}}) \sigma(d\hat{\mathbf{p}}) \delta_{\mathbf{x}}(dy),$$

with probability density function

$$\pi_{\varepsilon}(z, \hat{\mathbf{p}}) := \frac{\rho(\mathbf{x}, \hat{\mathbf{p}} \cdot \hat{\mathbf{k}})}{\int_{S_{\varepsilon}^{z}} \rho(\mathbf{x}, \hat{\mathbf{p}}' \cdot \hat{\mathbf{k}}) d\sigma(\hat{\mathbf{p}}')} \mathbf{1}_{S_{\varepsilon}^{z}}(\hat{\mathbf{p}}).$$

The above Dirac mass $\delta_{\mathbf{x}}(dy) := \delta(\mathbf{x} - \mathbf{y})dy$ translates the fact that the jumps only hold w.r.t. the $\hat{\mathbf{k}}$-variable.

From this construction we have the following probabilistic representation for the solution to (4.11).

**Proposition 4.4.2** The Markov process $D_{\varepsilon}$ defined in (4.14) has for infinitesimal generator

$$\mathcal{A}_{\varepsilon} g(z) := \hat{\mathbf{k}} \cdot \nabla_{\mathbf{x}} g(z) + \sigma_{\varepsilon}^2(\mathbf{x}) \Delta_{\mathbb{S}^2} g(z) + \Lambda_{\varepsilon}(z) \int_{\mathbb{S}^2} \pi_{\varepsilon}(z, \hat{\mathbf{p}}) (g(\mathbf{x}, \hat{\mathbf{p}}) - g(\mathbf{x}, \hat{\mathbf{k}})) \sigma(d\hat{\mathbf{p}}),$$

and we have

$$\mathbb{P}_{\mu_0}(D_{\varepsilon}(t) \in dx \otimes \sigma(d\hat{\mathbf{k}})) = \frac{1}{\bar{u}_0} u_{\varepsilon}(t, \mathbf{x}, \hat{\mathbf{k}}) dx \sigma(d\hat{\mathbf{k}}), \quad (4.16)$$

where

$$\mu_0(dx, d\hat{\mathbf{k}}) := \mathbb{P}(D_{\varepsilon}(0) \in dx \otimes \sigma(d\hat{\mathbf{k}})) = \frac{u_0(\mathbf{x}, \hat{\mathbf{k}})}{\bar{u}_0} dx \sigma(d\hat{\mathbf{k}}), \quad \text{with} \quad \bar{u}_0 := \int_{\mathbb{R}^3 \times \mathbb{S}^2} u_0(\mathbf{x}, \hat{\mathbf{k}}) dx \sigma(d\hat{\mathbf{k}}).$$
The terminology forward, for this MC method, comes from the fact that the particles are emitted at random points at time $t = 0$, through $\mu_0$ and the initial condition $u_0$, and propagate toward the observation positions $z = (\mathbf{x}, \hat{\mathbf{k}})$. This approach is used when we want to simulate the propagation over the entire space, or at least over regions large compared to the support of the initial condition. For the 2D case of Sect. 4.3, the probabilistic representation (4.8) corresponds to the backward MC method, where the particles are emitted at the observation point $z = (\mathbf{x}, \hat{\mathbf{k}})$ and are evaluated at time $t$ through the initial condition $u_0$ of the RTE. This approach is convenient when the initial condition has a support larger than the region over which we want to evaluate the solution of the RTE.

Let us illustrate two aspects of the representation (4.16). In order to obtain an estimation of $u_\epsilon(t, \mathbf{x}, \hat{\mathbf{k}})$ at point $z = (\mathbf{x}, \hat{\mathbf{k}})$, we evaluate the probability

$$ \frac{\bar{u}_0}{|\mathbf{B}(z,r)|} \mathbb{P}_{\mu_0}(D_\epsilon(t) \in \mathbf{B}(z,r)) \approx u_\epsilon(t, \mathbf{x}, \hat{\mathbf{k}}), $$

where $\mathbf{B}(z,r) \subset \mathbb{R}^3 \times \mathbb{S}^2$ stands for the ball centered at $z = (\mathbf{x}, \hat{\mathbf{k}})$ with radius $r \ll 1$. If we are only interested in the energy density at point $\mathbf{x}$ for instance, we estimate

$$ \frac{\bar{u}_0}{|\mathbf{B}(\mathbf{x}, r)|} \mathbb{P}_{\mu_0}(D_{\epsilon}(t) \in \mathbf{B}(\mathbf{x}, r) \times \mathbb{S}^2) = \frac{\bar{u}_0}{|\mathbf{B}(\mathbf{x}, r)|} \mathbb{P}_{\mu_0}(D_{1,\epsilon}(t) \in \mathbf{B}(\mathbf{x}, r)) \approx \int_{\mathbb{S}^2} u_{\epsilon}(t, \mathbf{x}, \hat{\mathbf{k}}) \sigma(d\hat{\mathbf{k}}), $$

where $\mathbf{B}(\mathbf{x}, r) \subset \mathbb{R}^3$ stands for the ball centered at $\mathbf{x}$ with radius $r \ll 1$, and $D_{1,\epsilon}$ is the $\mathbf{x}$-component of $D_\epsilon$.

### 4.4.2 The Monte-Carlo method

Based on the previous probabilistic representation (4.16), solving (4.11) requires the generation of random paths of the stochastic process $D_\epsilon$. For any measurable bounded function $f$, the convergence of the estimator

$$ \mu_N(t, f) := \frac{1}{N} \sum_{j=1}^{N} f(D_{\epsilon}^{j}(t)) \underset{N \to \infty}{\longrightarrow} \int f(\mathbf{x}, \hat{\mathbf{k}}) u_{\epsilon}(t, \mathbf{x}, \hat{\mathbf{k}}) d\mathbf{x} \sigma(d\hat{\mathbf{k}}) \quad \mathbb{P}_{\mu_0} \text{ - almost surely,} $$

is guaranteed by the strong law of large numbers. Above $(D_{\epsilon}^j)_j$ is a sample of $D_\epsilon$. The process $D_\epsilon$ being inhomogeneous, i.e. $\Lambda_\epsilon$ and $\Pi_\epsilon$ both depend on $z = (\mathbf{x}, \hat{\mathbf{k}})$, we make use of the so-called thinning method, also referred to as the fictitious shocks method [140], to simulate the jump part. It is based on a acceptation/rejection step and consists in simulating at first more jumps (or shocks) than necessary. In a second step, some of the jumps are rejected according to an appropriate probability distribution in order to recover the original dynamics.

The diffusion part between two jumps satisfies the linear SDE (4.15), and can be simulated using the following Euler-Maruyama type scheme

$$ (S_{n,m}) : \left\{ \begin{array}{l} X_{n,m+1} = X_{n,m} + h_{n,m} \hat{K}_{n,m} \\ K_{n,m+1} = \hat{K}_{n,m} - 2h_{n,m} \sigma_{\epsilon}^{2}(X_{n,m}) \hat{K}_{n,m} + \sqrt{2h_{n,m}} \sigma_{\epsilon}(X_{n,m}) \hat{K}_{n,m} \times W_{n,m} \\ \hat{K}_{n,m+1} = \frac{K_{n,m+1}}{|K_{n,m+1}|}, \end{array} \right. $$

where the $(W_{n,m})_{m,n}$ are i.i.d. mean-zero Gaussian random vectors with identity covariance matrix, and $h_{n,m}$ some appropriate stepsizes.

The theoretical justification of our MC method is supported by the following convergence result.

**Theorem 4.4.1** Consider

$$ \mu_{N,h,\epsilon}(t, f) = \frac{1}{N} \sum_{j=1}^{N} f(D_{\hat{h},\epsilon}^{j}(t)), \quad \text{and} \quad \mu(t, f) = \int f(\mathbf{x}, \hat{\mathbf{k}}) u(t, \mathbf{x}, \hat{\mathbf{k}}) d\mathbf{x} \sigma(d\hat{\mathbf{k}}), $$

where $(D_{\hat{h},\epsilon}^j)_j$ is a sample of the numerical approximation to $D_\epsilon$. For any $T > 0$, $\eta > 0$ and any smooth bounded function $f$ on $\mathbb{R}^3 \times \mathbb{S}^2$, we have

$$ \limsup_{N \to \infty} \mathbb{P} \left( |\mu_{N,h,\epsilon}(T, f) - \mu(T, f)| > \frac{\eta \Sigma_{h,\epsilon}}{\sqrt{N}} + \epsilon 2^{-(\alpha_M/2)} F_0 + hF_1 \right) \le \text{erfc}(\eta/\sqrt{2}), \quad (4.17) $$
where

$$ \Sigma_{h,\varepsilon} = \sqrt{\mathrm{Var}(f(D_{h,\varepsilon}(T)))} \leq \sup |f|, $$

and $F_0$ and $F_1$ are explicit constants independent of $\varepsilon$ and the stepsize $h$.

In (4.17), there is three terms that quantify the approximation error of our estimator $\mu_{N,h,\varepsilon}(t, f)$: one of order $\varepsilon^{2-(\alpha_M/2)}$ due to the approximation of $u$ by $u_\varepsilon$ (the smaller the $\alpha_M$, i.e. the less singular the kernel is, the smaller the error), one of order $h$ due to the numerical approximation of the diffusion over the unit sphere, and one due to the MC approximation with the standard $1/\sqrt{N}$ convergence rate. Note that the discretization error of the diffusion process is only of order $h$ and not of order the standard $\sqrt{h}$. The reason is that we are only interested in the convergence of MC estimators, allowing us to consider this discretization error in the weak sense [200]. Note that a weak second-order Runge-Kutta method can be considered to provide an error in $h^2$ instead of $h$ for the Euler scheme [59], and a modifications of the SDE (4.15) can also be considered to provide weak higher-order schemes [1].

### 4.4.3 Numerical illustrations

We first illustrate the role of the correction and the role of $\alpha$. We also present two situations, related to optical tomography and propagation through a non-Kolmogorov turbulent atmosphere, where $\alpha$ is no longer a constant. Furthermore, the stepsize $h$ is chosen proportionally to the shortest mean free time,

$$ \bar{\Lambda}_{\varepsilon}^{-1} := \frac{\alpha_m \varepsilon^{\alpha_M/2}}{2\pi \sup \lambda \sup a}, $$

and small enough so that the approximation error w.r.t. $\varepsilon$ in dominant. $N$ is also chosen large enough for the MC error to be small compared to the approximation error in $\varepsilon$. The numerical simulations below are performed using the Julia programming language (v1.6.5) on a NVIDIA Quadro RTX 6000 GPU driven by a 24 Intel Xeon Silver 2.20GHz CPUs station.

#### Role of the correction

To highlight the role of the correction provided by the diffusion on the unit sphere, the MC simulations, with and without the correction, are compared to a semi-analytic reference solution in the same spirit as in Sect. 4.3.

<figure><img src="image_14.png" /><figcaption>Fig. 4.7: Illustration of the relative error and running time of the MC method with and without a diffusive correction. The reference time in the right picture is the one of the corrected method with $\varepsilon = 0.1$. The observation time $T = 3t_S$ ($t_S$ computed for $\alpha = 1$) and $N = 6 \times 10^8$.</figcaption></figure>

In Fig. 4.7, we represent the relative error for various sizes of the cutoff $\varepsilon$. The left picture illustrates the evolution of the relative error for various $\varepsilon$. The blue curve corresponds to the corrected MC method with $\varepsilon = 0.1$ (with still a fairly large stepsize $h = 0.5/\bar{\Lambda}_\varepsilon$) providing at most a relative error slightly larger than 1%. The other curves correspond to the noncorrected MC method for several values of $\varepsilon$. The corrected MC method consistently yields a better accuracy compared to the noncorrected versions, and even in weakly singular cases where $\alpha$ is less than one. A very small value of $\varepsilon$ (orange and green curves) is necessary to match the accuracy of the corrected method. The right
picture illustrates the evolution of the relative running time of the noncorrected method w.r.t. the corrected one. For values of $\alpha$ less than 0.7 (weakly singular kernels), corrected and noncorrected methods have similar computational times for comparable accuracy, while in the case of singular kernels, with $\alpha \ge 1$, the noncorrected method yields a much larger cost and a much lower accuracy.

<figure><img src="image_15.png" /><figcaption>Fig. 4.8: Illustration of the relative error and running time of the corrected MC method.</figcaption></figure>

In Fig. 4.8, we illustrate the precision and running time sensitivity of the (corrected) MC method w.r.t. the stepsize $h = h_0/\bar{\Lambda}_\epsilon$. As expected, we obtain a better precision for smaller step sizes but at the price of a longer running time. These effects are amplified as $\alpha$ increases due to the increasing strength of the diffusion correction. In what follows, we select $h_0 = 0.3$ since this value yields a relative error less than 1% for a wide range of $\alpha$'s while not changing significantly the running time.

<figure><img src="image_16.png" /><figcaption>Fig. 4.9: Illustration of the relative error and running time of the corrected MC method.</figcaption></figure>

In Fig. 4.9, we depict the precision and running time sensitivity w.r.t. the cutoff parameter $\epsilon$, and observe the same phenomena as in the case of the stepsize $h$. With our choice of $h$, the parameter $\epsilon$ defines not only the accuracy of the diffusion correction, but also the average number of jumps, and as a consequence the running time increases as $\epsilon$ decreases as in the case of the noncorrected MC method.

## The role of $\alpha$

In this section, we highlight the effects of the kernel singularity on the energy density. We consider here a constant $\alpha$'s, with $a \equiv 0.002$. Our setting is depicted in Fig. 4.10. The spatial variable $\mathbf{x}$ is decomposed into a main propagation axis $x_3$ and a transverse plane $\mathbf{x}_{\perp}$, i.e. $\mathbf{x} = (\mathbf{x}_{\perp}, x_3) \in \mathbb{R}^2 \times \mathbb{R}$. The same notation holds for the direction variable $\hat{\mathbf{k}} = (\hat{\mathbf{k}}_{\perp}, \hat{k}_3) \in \mathbb{S}^2$. We choose an initial condition for the RTE of the form

$$u_0(\mathbf{x}, \hat{\mathbf{k}}) = \delta(\mathbf{x})\delta(\hat{\mathbf{k}} - \hat{k}_0), \quad \text{with} \quad \hat{k}_0 = (0, 0, 1),$$

modeling a source located at $\mathbf{x} = 0$, embedded in the random medium, and emitting in the forward $x_3$-direction. We set a function $\lambda$ of the form $\lambda(\mathbf{x}) = 1_{(-5,40)}(x_3)$, that defines a scattering layer between $x_3 = -5$ and $x_3 = 40$. In such a configuration, both transmitted and reflected quantities at
<figure><img src="image_17.png" /><figcaption>Fig. 4.10: Illustration of the numerical setting.</figcaption></figure>

$x_3 = 40$ and $x_3 = -5$ are of interest. With our particular choice of $\hat{k}_0$, what is observed at $x_3 = -5$ is purely due to backscattering.

The observables, that are estimated by the MC method, are the following: the (time-integrated) transverse reflected and transmitted energy respectively defined by

$$F_{tr}^{T}(\mathbf{x}_{\perp}) := \int_{0}^{T} dt \int_{S^{2}} \sigma(d\hat{\mathbf{k}}) u(t, \mathbf{x}_{\perp}, x_{3} = 40, \hat{\mathbf{k}})$$

and

$$F_{ref}^{T}(\mathbf{x}_{\perp}) := \int_{0}^{T} dt \int_{S^{2}} \sigma(d\hat{\mathbf{k}}) u(t, \mathbf{x}_{\perp}, x_{3} = -5, \hat{\mathbf{k}}),$$

and the reflected and transmitted time evolution of the integrated energy exiting the slab by

$$\mathcal{F}_{tr}(t) := \int_{\mathbb{R}^2} d\mathbf{x}_{\perp} \int_{S^2} \sigma(d\hat{\mathbf{k}}) u(t, \mathbf{x}_{\perp}, x_3 = 40, \hat{\mathbf{k}})$$

and

$$\mathcal{F}_{ref}(t) := \int_{\mathbb{R}^2} d\mathbf{x}_{\perp} \int_{S^2} \sigma(d\hat{\mathbf{k}}) u(t, \mathbf{x}_{\perp}, x_3 = -5, \hat{\mathbf{k}}).$$

In what follows, the MC estimations of the observables are obtained using $N = 1 \times 10^9$ particles and a diffusion stepsize $h = 0.3/\Lambda_{\epsilon}$. We set $\epsilon = 0.01$ for the transmitted quantities, and $\epsilon = 0.1$ for the reflected ones. In the transmission case, and when $\epsilon$ is too large, the mean free time for the jumps is large as well and it is possible that particles escape the slab without undergoing any jumps, leading to inaccurate results. The choice $\epsilon = 0.01$ provides good results to simulate transmitted quantities. Nonetheless, a larger value of $\epsilon$ still provides a good accuracy for the reflected quantities since the particles exiting the slab at $x_3 = -5$ have necessarily undergone some number of jumps and diffusion steps.

The running times for the time-integrated transmitted ($\epsilon = 0.01$) and reflected ($\epsilon = 0.1$) signals for different values of $\alpha$ are the following:

<table><thead><tr><th>running time (s)</th><th>α = 0.3</th><th>α = 0.7</th><th>α = 0.1</th><th>α = 1.3</th><th>α = 1.5</th></tr></thead><tbody><tr><td>ε = 0.01</td><td>9.38</td><td>15.21</td><td>24.54</td><td>42.74</td><td>65.39</td></tr><tr><td>ε = 0.1</td><td>3.61</td><td>3.72</td><td>4.12</td><td>4.92</td><td>6.01</td></tr></tbody></table>

All these running time measurements account also for the transfer of the resulting arrays from the device (GPU) to the host, and correspond to the cost for the MC method to reach the expected accuracy for fixed $\epsilon$'s and $\alpha$'s. We clearly observe significantly larger running times for smaller values of $\epsilon$ and large values of $\alpha$ as expected. This is due to the increase in scattering events and diffusion steps as the mean free time decreases.

## Comparison with the Henyey-Greenstein scattering kernel

According to the definition of the so-called Henyey-Greenstein scattering kernel (4.4) together with the asymptotic (4.5), our MC method should be a candidate to simulate numerically RTE with such
a scattering kernel in the highly peaked forward regime. In this section, we compare the solution to the RTE with the Henyey-Greenstein scattering kernel, for an anisotropy factor *g* close to one, with the solution to (4.11) with a singular kernel derived from (4.5), that is by setting $a \equiv (1 - g)/(2\pi)$ and $\alpha = 1$ in (4.9). Note that the value of the constant *a* changes with *g*, and as a consequence $\bar{\Lambda}_{\epsilon}$, *h*, and $\sigma_{\epsilon}$ vary accordingly. To illustrate this approximation, we still consider the setting depicted in Fig. 4.10 and the various observables introduced in the previous sections, but now at time $T = 300$. We observe in Fig. 4.11 the very good agreement between the two solutions. The reflected signal

<figure><img src="image_18.png" /><figcaption>Fig. 4.11: Comparison of the observables obtained using the Henyey-Greenstein scattering kernel and our singular kernel with $\alpha = 1$, $T = 300$, $\epsilon = 0.01$ for the transmitted observables (left panels), and $\epsilon = 0.1$ for the reflected ones (right panels).</figcaption></figure>

is well captured by our method despite fairly large values of $\epsilon$ and $h$. Also, let us mention that the computational cost is decreasing as the anisotropic parameter $g$ is getting close to 1, as the overall jump intensity decreases in this case in the highly peaked regime $g \to 1$.

The RTE with a Henyey-Greenstein scattering kernel is simulated through a standard MC method. Compared to our method, its computational costs to simulate the transmitted and reflected observables are the following:

<table><thead><tr><th>running time (s)</th><th>g = 0.97</th><th>g = 0.98</th><th>g = 0.99</th></tr></thead><tbody><tr><td>HG kernel</td><td>8.3</td><td>7.6</td><td>6.0</td></tr><tr><td>singular kernel, ε = 0.01</td><td>13.9</td><td>6.7</td><td>2.0</td></tr><tr><td>singular kernel, ε = 0.1</td><td>2.5</td><td>1.5</td><td>0.7</td></tr></tbody></table>

According to this table, under our MC method, and for the reflected observables, lower computational times are observed for the three considered values of *g* compared to a standard MC method for the Henyey-Greenstein scattering kernel. Our MC method provides therefore an efficient tool to simulate numerically a RTE with such a scattering kernel. For transmitted observables, *g* needs to be quite close to one to provide a significant advantage to our method.

Now, we illustrate the influence of a varying $\alpha$ characterizing the strength of the singularity. We consider two situations, one inspired from optical tomography, and the second one from wave propagation through atmospheric turbulence.

## A two-stage model with a sphere

We keep the setting depicted in Fig. 4.10, and add a defect with a different value of $\alpha$. This defect is modeled by ball of radius 3 centered at the origin, with $\alpha$ sets to $\alpha_1$ within the ball, and to 1 in
the exterior of the ball. See Fig. 4.12 for an illustration. This situation models a biological tissue in which statistical properties are changing and define a region of interest for imaging.

<figure><img src="image_19.png" /><figcaption>Fig. 4.12: Illustration of the setting with $\lambda = \mathbf{1}_{\{x_3 \in (-5, 40)\}}$ and $\alpha(\mathbf{x}) = \alpha_1 \mathbf{1}_{\{x \in B\}} + 1 \cdot \mathbf{1}_{\{x \notin B\}}$ where $B$ is a ball centered at 0 with radius 3. The source is here located at $(0, 0, -5)$.</figcaption></figure>

We illustrate in Fig. 4.13 the impact of the defect on the observables introduced above. The impact is stronger on transmitted observables and quite significant, giving then the possibility to identify the defect inside the scattering medium. Reflected quantities tend to be less sensitive to the presence of the defect since a fraction of the signal is backscattered before reaching it.

<figure><img src="image_20.png" /><figcaption>Fig. 4.13: Illustration of the transmitted (left-hand-side) and reflected (right-hand-side) observables with $T = 300$, and $\varepsilon = 0.01$. For the two top pictures we set $x_2 = 0$ (solid lines), and $x_2 = 1.5$ for the top-left and $x_2 = 0.5$ for top-right picture (dotted lines).</figcaption></figure>

## Non-Kolmogorov turbulences

We consider one more time the setting depicted in Fig. 4.10, with the difference that $\lambda = \mathbf{1}_{\{x_3 \in (0, 20)\}}$ and that $\alpha$ takes three different large values depending on the altitude parametrized by $x_3$, see Fig. 4.14,

$$ \alpha(x_3) = 5/3 \cdot \mathbf{1}_{\{x_3 \le 2\}} + 4/3 \cdot \mathbf{1}_{\{2 < x_3 \le 8\}} + 1.9 \cdot \mathbf{1}_{\{8 < x_3\}} $$

The value 5/3 corresponds to a standard Kolmogorov turbulence model, while other values are associated with non-Kolmogorov turbulence models [31, 195, 212]. In these models, it is considered that
for altitudes higher than 8km, the atmospheric turbulence yields larger statistical patterns (which can be modeled through singular scattering kernels) than around the ground (0-2km). Hence, we set $\alpha = 1.9$ for altitudes greater than 8km. The function $a$ is no longer constant in these models, and for our illustrations we chose

$$a(r) = 0.002 \cdot \exp(-r^2/(2 \times 0.8^2)).$$

In Fig. 4.15, one can notice that non-Kolmogorov turbulence yields quite different signals compared to Kolmogorov turbulence, in particular for reflected quantities. As we can see, the higher the $\alpha$, the more diffuse is the signal which then enhances reflected signals. This explains the increased reflections in the non-Kolmogorov case.

<figure><img src="image_21.png" /><figcaption>Fig. 4.14: Illustration of a three stages $\alpha$-profile for a non-Kolmogorov phase function.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 4.15: Illustration of the transmitted (left-hand-side) and reflected (right-hand-side) observables with $T = 300$, and $\varepsilon = 0.01$. For the top two pictures we illustrate $x_2 = 0$ (solid lines) and $x_2 = 2$ for the top-left and $x_2 = 5$ for top-right picture (dot lines).</figcaption></figure>

## 4.5 Perspectives

The original motivation for studying radiative transfer models with a spatially dependent scattering kernel (4.9) is driven by the resolution of some inverse problems. These include source or reflector
localization in turbulent atmosphere, as well as the estimation of the map $\alpha$, which has applications in medical imaging or turbulence atmosphere. Inversion techniques for transport equations often rely on the singularities of the solution $u$ of the RTE (see [20, 21]). However, the smoothness of $u$ in the context of a singular scattering kernel suggests that the inversion could be more challenging compared to the situation of an integrable scattering kernel.

Although MC methods provide a relatively good accuracy at low cost, are easy to implement, and highly parallelizable, achieving higher accuracy often come at a significant increase in the computational time. Some variance reduction techniques might then be necessary to lower this cost. They generally amount to bias the particle trajectories toward regions of higher importance, and often result in important gains in terms of efficiency [142, 192].

Nowadays machine learning techniques are increasingly used to enhance the performance of numerical schemes and tackle inverse problems. Another perspective is therefore to evaluate the potential of machine learning techniques in simulating numerically RTEs, enabling variance reduction techniques, and solving radiative transfer based inverse problems.
Chapter 5

Boundary effects and weak localization in
radiative transfer

This chapter describes the results obtained in Adel Messaoudi's Ph.D. thesis [158] under the super-
vision of Régis Cottereau (Laboratory of Mechanics and Acoustics in Marseille) and myself. The
pictures presented in this chapter are all borrowed from this Ph.D. thesis. The two first sections of
this chapter correspond to [16, 18] in the publication list on page 6.

Radiative transfer models are usually derived for wave problems defined on the entire space $\mathbb{R}^3$
under the weak scattering regime [12, 45, 180]. For instance, the scalar wave equation,

$$
\partial_{tt}^2 p_{\varepsilon}(t, \mathbf{x}) - c^2 \left(\frac{\mathbf{x}}{\varepsilon}\right) \Delta_{\mathbf{x}} p_{\varepsilon}(t, \mathbf{x}) = 0 \quad (t, \mathbf{x}) \in (0, \infty) \times \Omega, \quad (5.1)
$$

with $\Omega = \mathbb{R}^3$,

$$
c^2\left(\frac{\mathbf{x}}{\varepsilon}\right) := c_0^2\left(1 + \sqrt{\varepsilon V}\left(\frac{\mathbf{x}}{\varepsilon}\right)\right), \qquad (5.2)
$$

and equipped with appropriate initial conditions at $t = 0$, can be recast as the following first-order
hyperbolic system:

$$
\epsilon \partial_t \mathbf{u}_\epsilon + \mathcal{A}_\epsilon \mathbf{u}_\epsilon = 0, \quad \text{where} \quad \mathcal{A}_\epsilon := - \begin{pmatrix} 0 & c^2(\mathbf{x}/\epsilon) \\ \epsilon^2 \Delta_\mathbf{x} & 0 \end{pmatrix},
$$

and

$$
\mathbf{u}_{\epsilon}(t, \mathbf{x}) := \begin{pmatrix} p_{\epsilon}(t, \mathbf{x}) \\ \epsilon c^{-2}(\mathbf{x}/\epsilon) \partial_t p_{\epsilon}(t, \mathbf{x}) \end{pmatrix} \quad (t, \mathbf{x}) \in [0, \infty) \times \Omega.
$$

The random field $V$ in (5.2) accounts for variations of the wave-speed profile, and is assumed to exhibit
some mixing properties. The energy carried by $\mathbf{u}_\epsilon$ can be described using the Wigner transform,

$$
W_{\epsilon}(t, \mathbf{x}, \mathbf{k}) := W[\mathbf{u}_{\epsilon}, \mathbf{u}_{\epsilon}](t, \mathbf{x}, \mathbf{k}) = \int e^{i\mathbf{k} \cdot \mathbf{y}} \mathbf{u}_{\epsilon}(\mathbf{x} - \epsilon\mathbf{y}/2) \otimes \mathbf{u}_{\epsilon}(\mathbf{x} + \epsilon\mathbf{y}/2) \frac{d\mathbf{y}}{(2\pi)^3}.
$$

The asymptotic behavior of this energy density is given by

$$
W_0(t, \mathbf{x}, \mathbf{k}) = \lim_{\epsilon \to 0} W_{\epsilon}(t, \mathbf{x}, \mathbf{k}) = a(t, \mathbf{x}, \mathbf{k})\mathbf{B}(\mathbf{k}) + a(t, \mathbf{x}, -\mathbf{k})\mathbf{B}^T(\mathbf{k}), \quad (5.3)
$$

where

$$
\mathbf{B}(\mathbf{k}) := \frac{1}{2} \begin{pmatrix} 1/|\mathbf{k}|^2 & i/(c_0|\mathbf{k}|) \\ -i/(c_0|\mathbf{k}|) & 1/c_0^2 \end{pmatrix}, \qquad (5.4)
$$

and $\mathbf{B}^T(\mathbf{k})$ stands for its transposition. Here, $\mathbf{B}$ corresponds, in the Fourier domain, to an eigenvector of a dispersion matrix associated to $\mathcal{A}_\epsilon$ with $c_0^2$ instead of $c^2(\mathbf{x}/\epsilon)$. The amplitude $a(t, \mathbf{x}, \mathbf{k})$ satisfies the following radiative transfer equation (RTE):

$$
\partial_t a(t, \mathbf{x}, \mathbf{k}) + c_0 \hat{\mathbf{k}} \cdot \nabla_{\mathbf{x}} a(t, \mathbf{x}, \mathbf{k}) = - \Sigma(\mathbf{k}) a(t, \mathbf{x}, \mathbf{k}) + \int \sigma(\mathbf{k}, \mathbf{q}) a(t, \mathbf{x}, \mathbf{q}) \delta(c_0(|\mathbf{q}| - |\mathbf{k}|)) d\mathbf{q}, \quad (5.5)
$$
with $\hat{\mathbf{k}} := \mathbf{k}/|\mathbf{k}|$, initial conditions of the form

$$a(t = 0, \mathbf{x}, \mathbf{k}) = \mathcal{A}(\mathbf{k})\delta(\mathbf{x} - \mathbf{x}_0), \qquad (5.6)$$

and scattering kernel and attenuation respectively given by

$$\sigma(\mathbf{k}, \mathbf{q}) := \frac{\pi c_0^2 |\mathbf{k}|^2}{2(2\pi)^3} \hat{R}(\mathbf{k} - \mathbf{q}), \quad \text{and} \quad \Sigma(\mathbf{k}) := \int \sigma(\mathbf{k}, \mathbf{q}) \delta(c_0(|\mathbf{q}| - |\mathbf{k}|)) d\mathbf{q}. \qquad (5.7)$$

Here, $\hat{R}$ stands for the power spectral density of the random fluctuations $V$. It is worth mentioning that the limiting equation is deterministic while the original wave problem is of stochastic nature. Due to the high-frequency regime under consideration, averaging phenomena occur resulting in a deterministic model for the energy propagation.

Despite being a significant challenge for engineering applications, radiative transfer models involving boundaries have received less attention [2, 22, 150, 161, 181]. The main configuration studied in the literature corresponds to a half-space ($\Omega = \mathbb{R}^2 \times \mathbb{R}^*$), where radiative transfer models remain valid. With either Dirichlet or Neumann boundary conditions on $\partial\Omega = \{x_n = 0\}$, the resulting boundary condition for the Wigner transform is a geometric optic type reflection condition:

$$W(t, \mathbf{x}_{\perp}, x_n = 0, \mathbf{k}_{\perp}, k_n) = W(t, \mathbf{x}_{\perp}, x_n = 0, \mathbf{k}_{\perp}, -k_n). \qquad (5.8)$$

In this chapter the spatial variable $\mathbf{x}$ is split into two components, $\mathbf{x} = (\mathbf{x}_{\perp}, x_n)$. Here, the $x_n$-variable represents the coordinate along the unit normal vector to the plane $\mathbb{R}^2 \times \{0\}$, while $\mathbf{x}_{\perp}$ denotes its transverse coordinate. Similarly, the wavevector is decomposed as $\mathbf{k} = (\mathbf{k}_{\perp}, k_n)$. Therefore, the relation (5.8) corresponds to a reflection condition w.r.t. the $x_n$-variable.

In what follows three propagation domains are discussed: the half-space, the slab, and the rectangle, all considered with Neumann boundary conditions for simplicity. For each case, the energy density propagation is described using RTEs, similar to (5.5), with reflection conditions (5.8) along each boundaries. Additionally, we present interference effects resulting from the presence of these boundaries. These effects can take place along the boundaries or within the propagation domain. Our approach is based on a method of images, which allows to extend the propagation domain over the full physical space $\mathbb{R}^3$ and conveniently represents these interference effects using the standard Wigner transform.

## 5.1 Boundary effects for a half-space

For the half-space $\Omega = \mathbb{R}^2 \times \mathbb{R}^*$, the wave equation (5.1) is equipped with Neumann boundary conditions,

$$\partial_{x_n} p_{\varepsilon}(t, \mathbf{x}_{\perp}, x_n = 0) = 0 \quad (t, \mathbf{x}_{\perp}) \in (0, \infty) \times \mathbb{R}^2,$$

and initial conditions of the form

$$p_{\varepsilon}(t = 0, \mathbf{x}) = \frac{1}{\varepsilon^{3/2}} A\left(\frac{\mathbf{x} - \mathbf{x}_0}{\varepsilon}\right) \quad \text{and} \quad \partial_t p_{\varepsilon}(t = 0, \mathbf{x}) = \frac{1}{\varepsilon^{3/2}} B\left(\frac{\mathbf{x} - \mathbf{x}_0}{\varepsilon}\right), \qquad (5.9)$$

where $\mathbf{x}_0$ represents the source location in $\Omega$. Here, the scaling in $\varepsilon$ allows a nontrivial limit for the Wigner transform. The functions $A$ and $B$ are assumed to be even w.r.t the $x_n$-variable to provide convenient symmetry properties, and compactly supported so that the supports of the initial conditions do not cross the interface $\{x_n = 0\}$. In this chapter, we only consider the situation of initial conditions supported at a distance of order 1 from the boundaries. Situations where the source in close to the border (e.g. $x_{0,n} \to \varepsilon x_{0,n}$ for instance) and Dirichlet boundary conditions are addressed in [158].

### 5.1.1 The method of images

The basic principle of the method of images is to replace the wave propagation problem posed on a half-space with two new ones posed on the full-space:

$$\partial_{tt}^2 p_{\varepsilon}^{\pm}(t, \mathbf{x}) - c_{\#}^2(\varepsilon)(\mathbf{x}) \Delta_{\mathbf{x}} p_{\varepsilon}^{\pm}(t, \mathbf{x}) = 0 \quad (t, \mathbf{x}) \in (0, \infty) \times \mathbb{R}^3,$$
<figure><img src="image_23.png" /><figcaption>Fig. 5.1: Illustration of the method of images for the half-space bounded above by the dashed line representing {$x_n = 0$}. $p^-$ represents an up-going wave propagating in a full-space from the original source location (in green). $p^+$ represents the down-going wave propagating from a fictitious source location resulting from the method of images (in red). The rightmost picture illustrates the original wave $p$ propagating in the half-space, from the original source location, as a superposition of $p^+$ and $p^-$.</figcaption></figure>

with symmetric initial conditions,

$$p_{\varepsilon}^{-}(t=0, \mathbf{x}) = \frac{1}{\varepsilon^{3/2}} A\left(\frac{\mathbf{x}_{\perp} - \mathbf{x}_{0,\perp}}{\varepsilon}, \frac{x_n - x_{0,n}}{\varepsilon}\right), \quad p_{\varepsilon}^{+}(t=0, \mathbf{x}) = \frac{1}{\varepsilon^{3/2}} A\left(\frac{\mathbf{x}_{\perp} - \mathbf{x}_{0,\perp}}{\varepsilon}, \frac{x_n + x_{0,n}}{\varepsilon}\right),$$

and

$$\partial_t p_{\varepsilon}^{-}(t=0, \mathbf{x}) = \frac{1}{\varepsilon^{3/2}} B\left(\frac{\mathbf{x}_{\perp} - \mathbf{x}_{0,\perp}}{\varepsilon}, \frac{x_n - x_{0,n}}{\varepsilon}\right), \quad \partial_t p_{\varepsilon}^{+}(t=0, \mathbf{x}) = \frac{1}{\varepsilon^{3/2}} B\left(\frac{\mathbf{x}_{\perp} - \mathbf{x}_{0,\perp}}{\varepsilon}, \frac{x_n + x_{0,n}}{\varepsilon}\right).$$

The extended wave-speed profile reads

$$c_{\# , \varepsilon}(\mathbf{x}) := c\left(\frac{\mathbf{x}_{\perp}}{\varepsilon}, \frac{x_n}{\varepsilon}\right) \mathbf{1}_{\mathbb{R}_{-}}(x_n) + c\left(\frac{\mathbf{x}_{\perp}}{\varepsilon}, -\frac{x_n}{\varepsilon}\right) \mathbf{1}_{\mathbb{R}_{+}^{*}}(x_n) \quad \mathbf{x} \in \mathbb{R}^3.$$

The solution $p_{\varepsilon}^{-}$ is referred to as the up-going wave, and $p_{\varepsilon}^{+}$ to as the down-going wave. Note that $c_{\# , \varepsilon}$ is an even function w.r.t. the $x_n$-variable so that

$$p_{\varepsilon}^{+}(t, \mathbf{x}_{\perp}, x_n) = p_{\varepsilon}^{-}(t, \mathbf{x}_{\perp}, -x_n), \qquad (5.10)$$

for all $(t, \mathbf{x}_{\perp}, x_n) \in [0, \infty) \times \mathbb{R}^2 \times \mathbb{R}$. Finally, the extended wave problem is obtained by considering $p_{\varepsilon}^{\#}$ defined as

$$p_{\varepsilon}^{\#}(t, \mathbf{x}) := p_{\varepsilon}^{-}(t, \mathbf{x}) + p_{\varepsilon}^{+}(t, \mathbf{x}) \quad (t, \mathbf{x}) \in \mathbb{R}_{+} \times \mathbb{R}^{3}.$$

The symmetry relation (5.10) leads to the boundary conditions $\partial_{x_n} p_{\varepsilon}^{\#}(t, \mathbf{x}_{\perp}, x_n = 0) = 0$, so that, restricted to the half-space, the extended wavefield $p_{\varepsilon}^{\#}$ corresponds to the original one $p_{\varepsilon}$:

$$p_{\varepsilon}^{\#}(t, \mathbf{x}) = p_{\varepsilon}(t, \mathbf{x}),$$

for any $(t, \mathbf{x}) \in [0, \infty) \times \Omega$.

The wave equations associated to the up- and down-going waves can be recast as first order hyperbolic systems:

$$\varepsilon \partial_t \mathbf{u}_{\varepsilon}^{j} + \mathcal{A}_{\varepsilon} \mathbf{u}_{\varepsilon}^{j} = 0, \quad \text{where} \quad \mathcal{A}_{\varepsilon} := - \begin{pmatrix} 0 & c_{\# , \varepsilon}^{2}(\mathbf{x}) \\ \varepsilon^{2} \Delta & 0 \end{pmatrix},$$

with

$$\mathbf{u}_{\varepsilon}^{j}(t, \mathbf{x}) := \begin{pmatrix} p_{\varepsilon}^{j}(t, \mathbf{x}) \\ \varepsilon c_{\# , \varepsilon}^{-2}(\mathbf{x}) \partial_{t} p_{\varepsilon}^{j}(t, \mathbf{x}) \end{pmatrix} \quad (t, \mathbf{x}) \in [0, \infty) \times \mathbb{R}^{3},$$

for $j \in \{-, +\}$. A similar system can be obtained for $p_{\varepsilon}^{\#}$ by setting

$$\mathbf{u}_{\varepsilon}^{\#} := \mathbf{u}_{\varepsilon}^{-} + \mathbf{u}_{\varepsilon}^{+}.$$
The energy density carried by $\mathbf{u}_{\varepsilon}^{\sharp}$ is described through its Wigner transform

$$W_{\varepsilon}^{\sharp} := W[\mathbf{u}_{\varepsilon}^{\sharp}, \mathbf{u}_{\varepsilon}^{\sharp}],$$

and the energy density associated to the original problem is assumed to be given by the restriction of $W_{\varepsilon}^{\sharp}$ over the original half-space.

Note that $c_{\sharp, \varepsilon}^2$ can be written as (5.2), where the fluctuations around the background wave speed $c_0$ are of the form

$$V_{\sharp}\left(\mathbf{x}, \frac{\mathbf{x}}{\varepsilon}\right) = V_{\sharp}(\mathbf{x}, \mathbf{y})|_{\mathbf{y}=\frac{\mathbf{x}}{\varepsilon}} := V(\mathbf{y}_{\perp}, y_n)\mathbf{1}_{\mathbb{R}_{-}}(x_n) + V(\mathbf{y}_{\perp}, -y_n)\mathbf{1}_{\mathbb{R}_{+}^{*}}(x_n)|_{\mathbf{y}=\frac{\mathbf{x}}{\varepsilon}}.$$

For any fixed slow component $\mathbf{x}$, the power spectral density of $V_{\sharp}(\mathbf{x}, \cdot)$, w.r.t. the fast component $\mathbf{y}$, reads

$$\mathbb{E}[\hat{\nu}_{\sharp}(\mathbf{x}, \mathbf{p})\hat{\nu}_{\sharp}(\mathbf{x}, \mathbf{q})] = (2\pi)^3 \delta(\mathbf{p} + \mathbf{q})\hat{R}(\mathbf{p}), \quad (5.11)$$

where $\hat{V}_{\sharp}(\mathbf{x}, \mathbf{p})$ stands for the Fourier transform of $V_{\sharp}(\mathbf{x}, \cdot)$ w.r.t. the $\mathbf{y}$-variable. Due to the symmetries, the resulting power spectral density does not depend on the slow component $\mathbf{x}$, and the method of images plays no role in the limiting RTE.

### 5.1.2 Radiative transfer model and boundary effects

The Wigner transform for $\mathbf{u}_{\varepsilon}^{\sharp}$ can be expanded as

$$W_{\varepsilon}^{\sharp}(t, \mathbf{x}, \mathbf{k}) = W[\mathbf{u}_{\varepsilon}^{\sharp}, \mathbf{u}_{\varepsilon}^{\sharp}](t, \mathbf{x}, \mathbf{k}) = W_{\varepsilon}^{s}(t, \mathbf{x}, \mathbf{k}) + W_{\varepsilon}^{s}(t, \mathbf{x}_{\perp}, -x_n, \mathbf{k}_{\perp}, -k_n) \\ + W_{\varepsilon}^{c}(t, \mathbf{x}, \mathbf{k}) + W_{\varepsilon}^{c}(t, \mathbf{x}_{\perp}, -x_n, \mathbf{k}_{\perp}, -k_n),$$

where, using the symmetry relation (5.10),

$$W_{\varepsilon}^{s} := W[\mathbf{u}_{\varepsilon}^{-}, \mathbf{u}_{\varepsilon}^{-}] \quad \text{and} \quad W_{\varepsilon}^{c} := W[\mathbf{u}_{\varepsilon}^{-}, \mathbf{u}_{\varepsilon}^{+}].$$

The term $W_{\varepsilon}^{s}$ is referred to as self-Wigner transform, and $W_{\varepsilon}^{c}$ to as cross-Wigner transform.

To describe the energy density carried by $\mathbf{u}_{\varepsilon}^{\sharp}$, it is necessary to analyze the asymptotic behavior of both the self-Wigner transform and cross-Wigner transform. Thanks to the method of images the asymptotic behavior of the self-Wigner transform remains as for the full space (5.3). However, when considering the cross-Wigner transform, the same asymptotic analysis to that of the self-Wigner transform can be applied. This strategy results as well in a standard (linear) RTE, but this time with null initial conditions,

$$\lim_{\varepsilon \to 0} W_{\varepsilon}^{c}(t=0) = 0,$$

as the supports of the initial conditions of $\mathbf{u}_{\varepsilon}^{-}$ and $\mathbf{u}_{\varepsilon}^{+}$ are disjoint. Therefore, the contribution of the cross-Wigner transform is here negligible:

$$\lim_{\varepsilon \to 0} W_{\varepsilon}^{c}(t) = 0,$$

for any time $t > 0$. As a result, the energy density carried by $\mathbf{u}_{\varepsilon}^{\sharp}$ is given by

$$W^{\text{tot}}(t, \mathbf{x}, \mathbf{k}) := \lim_{\varepsilon \to 0} W_{\varepsilon}^{\sharp}(t, \mathbf{x}, \mathbf{k}) = W_{0}(t, \mathbf{x}, \mathbf{k}) + W_{0}(t, \mathbf{x}_{\perp}, -x_{n}, \mathbf{k}_{\perp}, -k_{n}) \\ = (a(t, \mathbf{x}, \mathbf{k}) + a(t, \mathbf{x}_{\perp}, -x_{n}, \mathbf{k}_{\perp}, -k_{n}))\mathbf{B}(\mathbf{k}) \\ + (a(t, \mathbf{x}, -\mathbf{k}) + a(t, \mathbf{x}_{\perp}, -x_{n}, -\mathbf{k}_{\perp}, k_{n}))\mathbf{B}^{T}(\mathbf{k}),$$

where the amplitude $a$ satisfies the RTE (5.5) defined over the full-space $\mathbb{R}^3$, and $\mathbf{B}$ is given by (5.4). The resulting asymptotic energy density $W^{\text{tot}}$ satisfies the reflection condition (5.8), and describes the energy propagation over the entire propagation domain. To capture interference effects we need to be a bit more careful.

The idea behind the interference phenomena is the following. Considering two mean-zero random variables $X$ and $Y$ with variance $s^2$, one can draw an analogy between the variance $\mathbb{E}[(X + Y)^2]$ and
the total energy density of the extended wave field $p_{\epsilon}^{\sharp}$, which is the superposition of $p_{\epsilon}^{-}$ and $p_{\epsilon}^{+}$. By expanding the square and the expectation, we obtain:

$$ \mathbb{E}[(X + Y)^2] = E[X^2] + \mathbb{E}[Y^2] + 2\mathbb{E}[XY] = 2s^2(1 + \cos(\phi)), $$

with $\phi \in [0, 2\pi)$. Therefore, the manifestation of the amplification phenomenon depends on the correlation between the two random variables. More precisely, for two independent random variables, no phenomenon is observed; for two perfectly positively correlated variables, the mean square of $X+Y$ is doubled. To observe some correlations between $p_{\epsilon}^{-}$ and $p_{\epsilon}^{+}$, it should necessarily be at the interface $\{x_n = 0\}$, where the two waves are crossing. These correlations can be captured by the cross-Wigner transform that can be rewrite as

$$ W_{\epsilon}^{c}(t, \mathbf{x}, \mathbf{k}) = \frac{2}{\epsilon(2\pi)} \iint e^{-2ik_n y_n / \epsilon} e^{2ip_n x_n / \epsilon} W_{\epsilon}^{s}(t, \mathbf{x}_{\perp}, y_n, \mathbf{k}_{\perp}, p_n) dp_n dy_n. \quad (5.12) $$

The presence of the highly oscillatory terms $e^{-2ik_n/\epsilon}$ and $e^{2ip_n x_n/\epsilon}$ suggests to focus our attention at the vicinity of the interface $\{x_n = 0\}$ (within one wavelength), leading to the two changes of variables:

$$ \mathbf{x} \to \mathbf{x}_{\epsilon} := (\mathbf{x}_{\perp}, \epsilon \tilde{x}_n / 2) \quad \text{and} \quad \mathbf{k} \to \mathbf{k}_{\epsilon} := (\mathbf{k}_{\perp}, \epsilon \tilde{k}_n / 2). $$

This first change of variable allows to focus near the boundary $\{x_n = 0\}$, and the second one corresponds to directions approximately parallel to this plane.

As a result, the asymptotic (nontrivial) contribution of $W_{\epsilon}^{c}$ is given by

$$ \lim_{\epsilon \to 0} \epsilon W_{\epsilon}^{c}(t, \mathbf{x}_{\epsilon}, \mathbf{k}_{\epsilon}) = \frac{2}{(2\pi)} \int e^{-i\tilde{k}_{n}y_{n}} e^{ip_{n}\tilde{x}_{n}} W_{0}(t, \mathbf{x}_{\perp}, y_{n}, \mathbf{k}_{\perp}, p_{n}) dp_{n} dy_{n}, $$

with $W_0$ defined by (5.3). This yields the total energy at the vicinity of the boundary $\{x_n = 0\}$:

$$ E_{\text{boundary}}(t, \mathbf{x}_{\perp}, \tilde{x}_n) = 2 \int a(t, \mathbf{x}_{\perp}, x_n = 0, \mathbf{k})(1 + \cos(k_n \tilde{x}_n)) \mathbf{D}(\mathbf{k}) d\mathbf{k}, $$

with

$$ \mathbf{D}(\mathbf{k}) := \begin{pmatrix} 1/|\mathbf{k}|^2 & 0 \\ 0 & 1/c_0^2 \end{pmatrix}. \quad (5.13) $$

Finally, for $\tilde{x}_n = 0$, that is exactly at the boundary, we have the relation:

$$ E_{\text{boundary}}(t, \mathbf{x}_{\perp}, \tilde{x}_n = 0) = 2 \int W^{\text{tot}}(t, \mathbf{x}_{\perp}, x_n = 0) d\mathbf{k}, $$

corresponding to a doubling of the energy described at the macroscopic level.

To conclude this section, it's important to note that the underlying idea to exhibit interference effects involves the analyze of two waves propagating in opposite directions along the mediator plane of the line joining the two sources. In the case of the half-space, the two waves cross at the boundary. In the remaining of this chapter, this idea is behind all the interference effects we describe.

## 5.2 Weak localization phenomena for a slab

In this section, a slab is defined as a domain contained in between two parallel planes:

$$ \Omega = \mathbb{R}^2 \times (0, H). $$

The wave equation (5.1) is equipped with two Neumann boundary conditions,

$$ \partial_{x_n} p(t, \mathbf{x}_{\perp}, x_n = 0) = \partial_{x_n} p(t, \mathbf{x}_{\perp}, x_n = H) = 0, \quad (t, \mathbf{x}_{\perp}) \in (0, \infty) \times \mathbb{R}^2. $$

The initial conditions (5.9) are again assumed to be compactly supported within $\Omega$ to do not cross the boundaries, and even w.r.t the $x_n$-variable to provide convenient symmetry properties.

As described below, the situation in a slab slightly differs from that in a half-space. Although similar interference effects can be observed at the boundaries of the slab, additional interferences occur
<figure><img src="image_24.png" /><figcaption>Fig. 5.2: Illustration of the method of images for the slab. The red crosses correspond to the source locations. Each wave propagates in a periodic full-space. Interference effects are observed within one wavelength of the boundaries (thick black lines), as well as along the dashed lines within the slab.</figcaption></figure>

between waves emitted in opposite directions and crossing each other within the slab after bouncing back and forth between the boundaries. We refer to these additional interference effects as weak localization phenomena in the sense of [47, 76]. Although this term is often associated with enhance backscattering produced by heterogeneous propagation media, this terminology is here associated to reflection produced by the boundaries of the propagation medium.

The methodology used to analyze this situation relies again on the method of images, but it is now combined with an appropriate periodization of the source and wave-speed profile.

### 5.2.1 The method of images

The main idea is here to replace the wave problem posed on the slab, with given initial conditions, by a family of problems set over the full-space (see Fig. 5.2). To this end, we introduce the following wave problem over the full-space:

$$ \partial_{tt}^2 q_{\varepsilon}(t, \mathbf{x}) - c_{\#}^2_{\varepsilon}(\mathbf{x}) \Delta_{\mathbf{x}} q_{\varepsilon}(t, \mathbf{x}) = 0 \quad (t, \mathbf{x}) \in (0, \infty) \times \mathbb{R}^3, $$

with initial conditions (5.9), and wave-speed profile

$$ c_{\#}^{2}(\mathbf{x}) = \sum_{j \in \mathbb{Z}} c\left(\frac{\mathbf{x}_{\perp}}{\varepsilon}, \frac{x_n - 2jH}{\varepsilon}\right) 1_{[0,H]}(x_n - 2jH) + c\left(\frac{\mathbf{x}_{\perp}}{\varepsilon}, \frac{-x_n + 2jH}{\varepsilon}\right) 1_{[-H,0]}(x_n - 2jH) \quad \mathbf{x} \in \mathbb{R}^3. $$

The extended wavefield over the full-space is then given by

$$ p_{\varepsilon}^{\#}(t, \mathbf{x}) := \sum_{j \in \mathbb{Z}} q_{\varepsilon}(t, \mathbf{x}_{\perp}, x_n - 2jH) + q_{\varepsilon}(t, \mathbf{x}_{\perp}, -(x_n - 2jH)) \quad (t, \mathbf{x}) \in [0, \infty) \times \mathbb{R}^3, \quad (5.14) $$

and through symmetry relations corresponds to the original wavefield on the slab:

$$ p_{\varepsilon}^{\#}(t, \mathbf{x}) = p_{\varepsilon}(t, \mathbf{x}), $$

for all $(t, \mathbf{x}) \in [0, \infty) \times \Omega$.

A first order hyperbolic system can be associated to $p_{\varepsilon}^{\#}$ and $q_{\varepsilon}$ by considering

$$ \mathbf{u}_{\varepsilon}(t, \mathbf{x}) := \left( \begin{array}{c} q_{\varepsilon}(t, \mathbf{x}) \\ \varepsilon c_{\#}^{-2}_{\varepsilon}(\mathbf{x}) \partial_t q_{\varepsilon}(t, \mathbf{x}) \end{array} \right) \quad \text{and} \quad \mathbf{u}_{\varepsilon}^{\#}(t, \mathbf{x}) := \left( \begin{array}{c} p_{\varepsilon}^{\#}(t, \mathbf{x}) \\ \varepsilon c_{\#}^{-2}_{\varepsilon}(\mathbf{x}) \partial_t p_{\varepsilon}^{\#}(t, \mathbf{x}) \end{array} \right), $$

for $(t, \mathbf{x}) \in [0, \infty) \times \mathbb{R}^3$. According to (5.14) both wavefield are related as follows,

$$ \mathbf{u}_{\varepsilon}^{\#}(t, \mathbf{x}) = \sum_{j \in \mathbb{Z}} \mathbf{u}_{\varepsilon}(t, \mathbf{x}_{\perp}, x_n - 2jH) + \mathbf{u}_{\varepsilon}(t, \mathbf{x}, -x_n + 2jH) \quad (t, \mathbf{x}) \in [0, \infty) \times \mathbb{R}^3. $$
### 5.2.2 The self- and cross-Wigner transforms

The energy density carried by $p_{\varepsilon}^{\#}$, represented by the Wigner transform of $\mathbf{u}_{\varepsilon}^{\#}$, can be expanded as

$$W_{\varepsilon}^{\#} := W[\mathbf{u}_{\varepsilon}^{\#}, \mathbf{u}_{\varepsilon}^{\#}] = W_{\varepsilon}^{1, \#} + W_{\varepsilon}^{2, \#} + W_{\varepsilon}^{3, \#},$$

with

$$W_{\varepsilon}^{1, \#}(t, \mathbf{x}, \mathbf{k}) := \sum_{j \in \mathbb{Z}} W_{\varepsilon}(t, \mathbf{x}_{\perp}, x_n - 2jH, \mathbf{k}) + W_{\varepsilon}(t, \mathbf{x}_{\perp}, -x_n + 2jH, \mathbf{k}_{\perp}, -k_n),$$

$$W_{\varepsilon}^{2, \#}(t, \mathbf{x}, \mathbf{k}) := \sum_{j, \ell \in \mathbb{Z}} V_{\varepsilon}^{j\ell}(t, \mathbf{x}, \mathbf{k}) + V_{\varepsilon}^{j\ell}(t, \mathbf{x}_{\perp}, -x_n, \mathbf{k}_{\perp}, -k_n),$$

$$W_{\varepsilon}^{3, \#}(t, \mathbf{x}, \mathbf{k}) := \sum_{\substack{j, \ell \in \mathbb{Z} \\ \ell \neq j}} W_{\varepsilon}^{j\ell}(t, \mathbf{x}, \mathbf{k}) + W_{\varepsilon}^{j\ell}(t, \mathbf{x}_{\perp}, -x_n, \mathbf{k}_{\perp}, -k_n),$$

where

$$W_{\varepsilon} := W[\mathbf{u}_{\varepsilon}, \mathbf{u}_{\varepsilon}],$$

refers to as the self-Wigner transform, and

$$V_{\varepsilon}^{j\ell}(t, \mathbf{x}, \mathbf{k}) := \int_{\mathbb{R}^3} e^{i\mathbf{k} \cdot \mathbf{y}} \mathbf{u}_{\varepsilon} \left( t, \mathbf{x}_{\perp} - \frac{\varepsilon \mathbf{y}_{\perp}}{2}, x_n - \frac{\varepsilon y_n}{2} - 2jH \right) \otimes \mathbf{u}_{\varepsilon} \left( t, \mathbf{x}_{\perp} + \frac{\varepsilon \mathbf{y}_{\perp}}{2}, -x_n - \frac{\varepsilon y_n}{2} + 2\ell H \right) \frac{d\mathbf{y}}{(2\pi)^3}$$

$$W_{\varepsilon}^{j\ell}(t, \mathbf{x}, \mathbf{k}) := \int_{\mathbb{R}^3} e^{i\mathbf{k} \cdot \mathbf{y}} \mathbf{u}_{\varepsilon} \left( t, \mathbf{x}_{\perp} - \frac{\varepsilon \mathbf{y}_{\perp}}{2}, x_n - \frac{\varepsilon y_n}{2} - 2jH \right) \otimes \mathbf{u}_{\varepsilon} \left( t, \mathbf{x}_{\perp} + \frac{\varepsilon \mathbf{y}_{\perp}}{2}, x_n + \frac{\varepsilon y_n}{2} - 2\ell H \right) \frac{d\mathbf{y}}{(2\pi)^3},$$

which both refer to as cross-Wigner transforms.

The cross-Wigner transform $V_{\varepsilon}^{j\ell}$ corresponds to waves whose initial conditions are located around $(\mathbf{x}_{0,\perp}, 2jH + x_{0,n})$ and $(\mathbf{x}_{0,\perp}, 2\ell H - x_{0,n})$ respectively. The mediator plane of the line joining these two sources, and along which interference effects take place, is $\{x_n = (j + \ell)H\}$. When restricted to the slab, $V_{\varepsilon}^{j\ell}$ provides nontrivial effects at the boundaries of the slab: $\{x_n = 0\}$ for $\ell = -j$, and $\{x_n = H\}$ for $\ell = 1 - j$.

Regarding the cross-Wigner transform $W_{\varepsilon}^{j\ell}$, it corresponds to waves whose initial conditions are located around either $(\mathbf{x}_{0,\perp}, 2jH + x_{0,n})$ and $(\mathbf{x}_{0,\perp}, 2\ell H + x_{0,n})$ in a first case, or around $(\mathbf{x}_{0,\perp}, -2jH - x_{0,n})$ and $(\mathbf{x}_{0,\perp}, -2\ell H - x_{0,n})$ in a second case. The mediator plane of the line joining these two sources, and along which interference effects take place, is $\{x_n = \pm x_{0,n} \pm (j+\ell)H\}$. When restricted to the slab, $W_{\varepsilon}^{j\ell}$ provides nontrivial effects along $\{x_n = x_{0,n}\}$ for $\ell = -j$ in the first case, and $\{x_n = H - x_{0,n}\}$ for $\ell = -1 - j$ in the second case. In contrast to $V_{\varepsilon}^{j\ell}$, the terms $W_{\varepsilon}^{j\ell}$ describe interference effects affecting only the coherent part of the energy. The reason for this difference, w.r.t. $V_{\varepsilon}^{j\ell}$, arise from the fact that when the two waves associated to $W_{\varepsilon}^{j\ell}$ intersect, they have experienced different random propagation media, while for $V_{\varepsilon}^{j\ell}$ they travel through the same medium by symmetry.

### 5.2.3 Boundary effects and weak localization phenomena

As for the half-space, without specific precaution, both cross-Wigner transforms $V_{\varepsilon}^{j\ell}$ and $W_{\varepsilon}^{j\ell}$ go to 0 as $\varepsilon \to 0$. In fact, the limits of $V_{\varepsilon}^{j\ell}$ and $W_{\varepsilon}^{j\ell}$ can be both described as solutions to standard (linear) transport equations with null initial conditions. As a result, the total energy density propagating over
the whole slab is only given by the contribution of $W_{\epsilon}^{1,\#}$, which is associated to self-Wigner transforms:

$$
\begin{align*}
W^{\text{tot}}(t, \mathbf{x}, \mathbf{k}) &:= \lim_{\epsilon \to 0} W_{\epsilon}^{\#}(t, \mathbf{x}, \mathbf{k}) \\
&= \lim_{\epsilon \to 0} W_{\epsilon}^{1,\#}(t, \mathbf{x}, \mathbf{k}) \\
&= \sum_{j \in \mathbb{Z}} W_0(t, \mathbf{x}_{\perp}, x_n - 2jH, \mathbf{k}) + W_0(t, \mathbf{x}_{\perp}, -x_n + 2jH, \mathbf{k}_{\perp}, -k_n) \\
&= \sum_{j \in \mathbb{Z}} (a(t, \mathbf{x}_{\perp}, x_n - 2jH, \mathbf{k}) + a(t, (\mathbf{x}_{\perp}, -x_n + 2jH, \mathbf{k}_{\perp}, -k_n)))\mathbf{B}(\mathbf{k}) \\
&\quad + \sum_{j \in \mathbb{Z}} (a(t, \mathbf{x}_{\perp}, x_n - 2jH, -\mathbf{k}) + a(t, \mathbf{x}_{\perp}, -x_n + 2jH, -\mathbf{k}_{\perp}, k_n))\mathbf{B}^{\text{T}}(\mathbf{k}),
\end{align*}
$$

according to (5.3). From this formulation, $W^{\text{tot}}$ satisfies the reflexion condition (5.10) at both $x_n = 0$ and $x_n = H$. The sum in $j$ accounts for the periodic behavior of the energy density within the slab.

Regarding the interference effects within one wavelength of the boundary $\{x_n = 0\}$, and using a representation similar to (5.12), we have now

$$
\begin{align*}
E_{\{x_n=0\}}(t, \mathbf{x}_{\perp}, \tilde{x}_n) &:= \int \lim_{\epsilon \to 0} W_{\epsilon}^{\#}(t, \epsilon \tilde{x}_n/2, \mathbf{k}) d\mathbf{k} \\
&= \int \lim_{\epsilon \to 0} W_{\epsilon}^{1,\#}(t, \epsilon \tilde{x}_n/2, \mathbf{k}) d\mathbf{k} \\
&\quad + \frac{1}{2} \sum_{j \in \mathbb{Z}} \int \lim_{\epsilon \to 0} \epsilon V_{\epsilon}^{j,\ell=-j}(t, \mathbf{x}_{\perp}, \epsilon \tilde{x}_n/2, \mathbf{k}_{\perp}, \epsilon \tilde{k}_n/2) \\
&\quad + \lim_{\epsilon \to 0} \epsilon V_{\epsilon}^{j,\ell=-j}(t, \mathbf{x}_{\perp}, -\epsilon \tilde{x}_n/2, \mathbf{k}_{\perp}, -\epsilon \tilde{k}_n/2) d\mathbf{k}_{\perp} d\tilde{k}_n,
\end{align*}
$$

yielding after some algebra to

$$
E_{\{x_n=0\}}(t, \mathbf{x}_{\perp}, \tilde{x}_n) = 2 \sum_{j \in \mathbb{Z}} \int a(t, \mathbf{x}_{\perp}, 2jH, \mathbf{k})(1 + \cos(k_n \tilde{x}_n)) \mathbf{D}(\mathbf{k}) d\mathbf{k},
$$

where **D** is given by (5.13). At exactly the boundary, we obtain again a doubling of the total energy,

$$
E_{\{x_n=0\}}(t, \mathbf{x}_{\perp}, \tilde{x}_n = 0) = 2 \int W^{\text{tot}}(t, \mathbf{x}_{\perp}, x_n = 0) d\mathbf{k}.
$$

The same type of result holds at the boundary $\{x_n = H\}$,

$$
\begin{align*}
E_{\{x_n=H\}}(t, \mathbf{x}_{\perp}, \tilde{x}_n) &= \int \lim_{\epsilon \to 0} W_{\epsilon}^{\#}(t, H + \epsilon \tilde{x}_n/2, \mathbf{k}) d\mathbf{k} \\
&= 2 \sum_{j \in \mathbb{Z}} \int a(t, \mathbf{x}_{\perp}, (2j+1)H, \mathbf{k})(1 + \cos(k_n \tilde{x}_n)) \mathbf{D}(\mathbf{k}) d\mathbf{k}.
\end{align*}
$$

yielding also a doubling of the total energy

$$
E_{\{x_n=H\}}(t, \mathbf{x}_{\perp}, \tilde{x}_n = 0) = 2 \int W^{\text{tot}}(t, \mathbf{x}_{\perp}, x_n = H) d\mathbf{k}.
$$

Regarding, the interference effects along the plane $\{x_n = x_{0,n}\}$, we have this time

$$
\begin{align*}
E_{\{x_n=x_{0,n}\}}(t, \mathbf{x}_{\perp}, \tilde{x}_n) &:= \int \lim_{\epsilon \to 0} W_{\epsilon}^{\#}(t, \mathbf{x}_{\perp}, x_{0,n} + \epsilon \tilde{x}_n/2, \mathbf{k}) d\mathbf{k} \\
&= \int \lim_{\epsilon \to 0} W_{\epsilon}^{1,\#}(t, \mathbf{x}_{\perp}, x_{0,n} + \epsilon \tilde{x}_n/2, \mathbf{k}) d\mathbf{k} \\
&\quad + \frac{1}{2} \sum_{j \in \mathbb{Z}^*} \iint \lim_{\epsilon \to 0} W_{\epsilon}^{j,\ell=-j}(t, \mathbf{x}_{\perp}, x_{0,n} + \epsilon \tilde{x}_n/2, \mathbf{k}_{\perp}, \epsilon \tilde{k}_n/2) d\mathbf{k}_{\perp} d\tilde{k}_n.
\end{align*}
$$
After some algebra, the total energy reads

$$
\begin{align*}
E_{\{x_n=x_{0,n}\}}(t, \mathbf{x}_{\perp}, \tilde{x}_n) &= \int W^{\text{tot}}(t, \mathbf{x}_{\perp}, x_n=x_{0,n}) d\mathbf{k} \\
&\quad + 2 \sum_{j \in \mathbb{Z}^*} \int \cos(k_n \tilde{x}_n) e^{-\Sigma(\mathbf{k})t} \mathbb{A}(\mathbf{k}) \delta(\mathbf{x}_{\perp} - \mathbf{x}_{0,\perp} - c_0 t \hat{\mathbf{k}}_{\perp}) \\
&\qquad \times \delta(2jH - c_0 t \hat{\mathbf{k}}_n) \mathbf{D}(\mathbf{k}) d\mathbf{k} \\
&\quad + 2 \sum_{j \in \mathbb{Z}^*} \int \sin(k_n \tilde{x}_n) e^{-\Sigma(\mathbf{k})t} \mathbb{A}(\mathbf{k}) \delta(\mathbf{x}_{\perp} - \mathbf{x}_{0,\perp} - c_0 t \hat{\mathbf{k}}_{\perp}) \\
&\qquad \times \delta(2jH - c_0 t \hat{\mathbf{k}}_n) \tilde{\mathbf{D}}(\mathbf{k}) d\mathbf{k},
\end{align*}
$$

with

$$
\tilde{\mathbf{D}}(\mathbf{k}) := \frac{1}{c_0|\mathbf{k}|} \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix},
$$

D given by (5.13), A by (5.6), and Σ(k) by (5.7). In this formulation, the Dirac masses δ(2jH − c
0
t
k
n
) account for the periodic onset of the interferences, which occurs at each time 2jH/c
0
for j ≥ 1. Also, each contribution is damped at rate Σ(k) (no multiple scattering) and freely transport along the plane {x
n
= x
0,n
} . These additional contributions are associated to the coherent part carried by the amplitude a , satisfying the RTE (5.5), along the plan {x
n
= x
0,n
} .

The same type of result holds for the plane { x_n = H - x_{0,n} } , where the total energy reads

$$
\begin{align*}
E_{\{x_n=H-x_{0,n}\}}(t, \mathbf{x}_{\perp}, \tilde{x}_n) &= \int W^{\text{tot}}(t, \mathbf{x}_{\perp}, x_n=H-x_{0,n}) d\mathbf{k} \\
&\quad + 2 \sum_{j \in \mathbb{Z}} \int \cos(k_n \tilde{x}_n) e^{-\Sigma(\mathbf{k})t} \mathbb{A}(\mathbf{k}) \delta(\mathbf{x}_{\perp} - \mathbf{x}_{0,\perp} - c_0 t \hat{\mathbf{k}}_{\perp}) \\
&\qquad \times \delta((2j+1)H - c_0 t \hat{\mathbf{k}}_n) \mathbf{D}(\mathbf{k}) d\mathbf{k} \\
&\quad + 2 \sum_{j \in \mathbb{Z}} \int \sin(k_n \tilde{x}_n) e^{-\Sigma(\mathbf{k})t} \mathbb{A}(\mathbf{k}) \delta(\mathbf{x}_{\perp} - \mathbf{x}_{0,\perp} - c_0 t \hat{\mathbf{k}}_{\perp}) \\
&\qquad \times \delta((2j+1)H - c_0 t \hat{\mathbf{k}}_n) \tilde{\mathbf{D}}(\mathbf{k}) d\mathbf{k}.
\end{align*}
$$

5.3 The case of a rectangle

The interference effects we observe for a rectangular domain [0, H₁] × [0, H₂] share similarities with those observed in the half-space and the slab. Nevertheless, due to the bounded nature of the propaga- tion domain, there are additional interference locations. In this section, we follow the same strategy, but do not provide any mathematical details on the interference effects as their descriptions are quite technical. Instead, we provide formal descriptions supported by some illustrations. A detailed treatment of this situation can be found in [158].

**The method of images** The method of images follows here a similar strategy to that used for the slab. The solution to the wave equation on the rectangle is extended over the full space R³, but in contrast to the slab, both coordinates are bounded and necessitate a periodic extension. The solution to the extended wave equation can be expressed as a superposition of waves whose initial conditions are located in a symmetric and periodic manner in order to enforce Neumann boundary conditions. We refer to Fig. 5.3 for an illustration. As a result, the extended solution reduces to the original wave problem when restricted to the rectangle.

**Self- and cross-Wigner transform** The analysis of the energy propagation is conducted using the Wigner transform of the extended solution. This Wigner transform can be expanded over two types of Wigner transform: the self- and cross-Wigner transform. As with the slab, the self-Wigner
<figure><img src="image_25.png" /><figcaption>Fig. 5.3: Illustration of the method of images for the rectangle $[0, H_1] \times [0, H_2]$. Each wave propagates from a different initial location (in red).</figcaption></figure>

transform accounts for energy propagation all over the rectangle, and can be described in terms of the standard RTE (5.5) over the full space:

$$ \int W^{\text{tot}}(t, \mathbf{x}, \mathbf{k}) d\mathbf{k} = \sum_{j_1, j_2 \in \mathbb{Z}^2} \sum_{\alpha_1, \alpha_2 \in \{+, -\}} \int a(t, \alpha_1(x_1 - 2j_1H_1), \alpha_2(x_2 - 2j_2H_2), \mathbf{k}) \mathbf{D}(\mathbf{k}) d\mathbf{k}, $$

where a satisfies (5.5) in $\mathbb{R}^2$ instead of $\mathbb{R}^3$, D is given by (5.13). Here, $W^{\text{tot}}$ satisfies reflexion conditions similar to (5.10) along each boundaries. The infinite sums in $j_1$ and $j_2$ account for the periodic behavior of the energy propagation when restricted to the rectangle w.r.t. both variables.

The cross-Wigner transforms capture interference effects over the propagation domain. As observed for the half-space and the slab, these effects occur along the mediator line of the line connecting the two source locations.

**Interference effects along the boundaries** As in the context of a half-space or a slab, constructive interference effects can be observed within one wavelength of each boundaries, resulting in a doubling of the total energy when looking exactly at the borders. In this case the whole energy is affected, not only its coherent part. We refer to Fig. 5.4 for an illustration.

**Interference effects along lines parallel to the boundaries** Similar to the slab, interference effects occur within one wavelength along lines parallel to the boundaries, passing through the source location and three other symmetric points. In this context, the two interacting waves have not experienced the same propagation medium, so that only the coherent part of the energy is affected. We refer to Fig. 5.5 for an illustration.

Due to the bounded nature of the propagation domain, additional symmetries between two source locations become possible.

**Interference effects along lines passing through the corners** The periodic extension of the rectangle in both coordinates allows to draw mediator lines, w.r.t. the lines connecting the source locations, that are now oblique when crossing the rectangle. As illustrated in Fig. 5.6, pairs of sources are able to produce interference effects along oblique lines passing through symmetric points, here the corners of the rectangle. In this context, the two interacting waves have travel through the same propagation medium, resulting in interferences that affect the whole energy, not only its coherent part. In contrast with the interference lines parallel to the boundaries, here, a given line supports interference effects only once over time. Despite the periodic nature of these effects, as they pass through the same corner, the lines supporting the interference effects rotate as the source locations producing these interferences are getting apart from each other.
<figure><img src="image_26.png" /><figcaption>Fig. 5.4: Illustration of the interference effects localized within one wavelength of the four boundaries (red dashed lines).</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 5.5: Illustration of the interference effects within one wavelength of the vertical lines {$x_1 = x_{0,1}$} and {$x_1 = H_1 - x_{0,1}$}, as well as the horizontal lines {$x_2 = x_{0,2}$} and {$x_2 = H_2 - x_{0,2}$} (the red dashed lines).</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Fig. 5.6: Illustration of an interference effect within one wavelength of the line passing through the lower left corner (red dashed lines).</figcaption></figure>

Interference effects along lines passing through the source location and symmetric points

Twelve additional symmetric points can be identified, all located along the boundaries and inside the rectangle (including the original source location). We refer to Fig. 5.7 for an illustration. As for the corners, and for the same reasons, oblique lines of interference are passing periodically through these points. However, a given line supports interference effects only once over time. Nonetheless, in this context, the two interacting waves have not experienced the same propagation medium, resulting in interference affecting only the coherent part of the energy.

<figure><img src="image_29.png" /><figcaption>Fig. 5.7: Illustration of an interference effect within one wavelength of the line passing through one of the twelve points (here $(H_1, H_2 - x_{0,2})$). The points are $x_0$ (the source location), and from the bottom to the top, $(x_{0,1}, 0)$, $(H_1 - x_{0,1}, 0)$, $(0, x_{0,2})$, $(H_1 - x_{0,1}, x_{0,2})$, $(H_1, x_{0,2})$, $(0, H_2 - x_{0,2})$, $(x_{0,1}, H_2 - x_{0,2})$, $(H_1 - x_{0,1}, H_2 - x_{0,2})$, $(H_1, H_2 - x_{0,2})$, $(x_{0,1}, H_2)$, $(H_1 - x_{0,1}, H_2)$.</figcaption></figure>

5.4 Perspectives

The understanding of the weak localization phenomena in the context of RTEs is a major challenge in nondestructive testing of concrete structures. These phenomena, driven by interference effects, offer the opportunity to extract valuable information from the enhancement intensity about the internal
structure and material properties. Numerical test should be performed to compare our theory with the work of [47, 76] conducted in homogeneous media. In particular, in heterogeneous media, the weak localization effects are exponentially damped, much like the total energy traveling over the whole domain. Therefore, quantifying the duration over which weak localization phenomena can be effectively leveraged in measurements before becoming negligible would be an interesting pursuit.

For engineering applications, the boundaries of the propagation medium play an important role. An interesting problem is to provide a precise description, in the context of elastic waves, of the surface modes (Rayleigh and Love) and their interplay with bulk modes. This relation could be incorporated into radiative transfer models, for applications in geophysics and nondestructive testing, and investigate how surface modes are related to the weak localization effects presented in this chapter.
Chapter 6

# Wave scattering by rough surfaces

Two aspects are discussed in this chapter. The first situation corresponds to waves propagating along a waveguide with rough-boundaries, and the second to waves that impinge upon a rough surface.

Acoustic wave propagation in waveguides has been studied for a long time due to its wide range of applications. One of its most important application is submarine detection, such as mines or archaeological artifacts detection. Additionally, it finds applications in underwater communication and the study of the ocean's structure or biology. Underwater acoustic waveguides are used to model propagation media such as a continental shelves. These environments are characterized by indexes of refraction that are both spatially and time dependent. Despite this complexity, the sound speed in water, approximately 1500m/s, remains sufficiently large with respect to the motions of water masses to consider this propagation medium as time independent. However, even small inhomogeneities in water conditions, such as the heave or rough ocean bottoms, can induce significant effects over large propagation distances.

When waves impinge upon a rough surface, they undergo scattering events due to the irregularities in the surface profile. These scattering events result in changes of the wave's direction, amplitude, and phase. Consequently, the study of wave scattering phenomena by rough surfaces involves unraveling the intricate details of these interactions to predict and analyze the behavior of waves in practical scenarios. In applications such as radar systems, electromagnetic wave scattering theory from rough surfaces is crucial to account for the effects of the roughness of the land or the sea surface. A similar theoretical framework can be applied in optics. For acoustic waves, one can refer to applications in ocean acoustic tomography or thermometry, where the sea bed and sea surface roughness play a critical role, but also in non-destructive testing using ultrasonic waves.

This chapter describes both of these scenarios. First, we discuss the situation where a wave propagates along a random surface in context of a waveguide. Then, we discuss the classical situation where a wave impinges upon a rough surface, producing reflected and transmitted components. This chapter covers the results obtained in [7, 19] of the publication list on page 6.

## 6.1 A Pekeris waveguide model with a rough surface and bottom topography

Mathematical studies regarding acoustic wave propagation in randomly perturbed waveguides have been carried out in many papers [70, 82, 91] under the assumption of rigid waveguide boundaries and random perturbations of its index of refraction. This section focuses on a two-dimensional acoustic waveguide model with a Pekeris profile, a randomly perturbed free surface, and uneven bottom topography (see Fig. 6.1). Our approach is built on a conformal transformation providing a smooth change of coordinates and allowing the use of the modal decomposition of the unperturbed waveguide (Fig. 6.1 left picture). In this context, a propagating wavefield can be decomposed over three kinds of mode: the propagating modes, which propagate over long distances along the waveguide; the evanescent modes, whose amplitudes decrease exponentially with the propagation distance along the waveguide; and the radiating modes, representing modes that can penetrate beneath the ocean floor.
<figure><img src="image_30.png" /><figcaption>Fig. 6.1: Illustration of two semi-infinite waveguides with a Pekeris profile. The left picture depicts an unperturbed waveguide, while the right picture illustrates a perturbed waveguide with a free surface and an uneven bottom topography. In both cases, the bulk properties of the waveguides are homogeneous with a propagation speed $c_1$ in the oceanic section of the waveguide, and $c_0$ in the bottom part of the waveguide.</figcaption></figure>

Using an asymptotic analysis, corresponding to a high-frequency regime, the impact of the randomly perturbed waveguide geometry on the mode coupling can be described in term of a diffusion model.

These mechanisms result in an effective attenuation of the mode amplitudes. It turns out that the surface and bottom fluctuations mainly affect the propagating-mode amplitudes in a similar manner. However, for the majority of the propagating modes, the attenuation mechanism is mainly due to the coupling between the propagating mode themselves. In contrast, for the highest propagating modes, the attenuation mechanism is mainly due to the coupling with the radiating modes, and it is significantly stronger than for the lower propagating modes.

### 6.1.1 Waveguide model

A two-dimensional linear acoustic wave model is considered:

$$ \Delta - \frac{1}{c_{\varepsilon}^2(z, x)} \partial_{tt}^2 p = \nabla \cdot F \quad (t, z, x) \in \mathbb{R} \times \mathbb{R} \times (0, \infty), \quad (6.1) $$

with $\Delta = \partial_{xx}^2 + \partial_{zz}^2$. Here, the coordinate $z$ represents the propagation axis along the waveguide, and the coordinate $x$ represents the transverse section of the waveguide (see Fig. 6.1). Let $d > 0$ be the average ocean depth. The velocity field is given, for $z \in \mathbb{R}$, by

$$ c_{\varepsilon}(z, x) = \begin{cases} c_1 & \text{if } V_s^{\varepsilon}(z) < x < d + V_b^{\varepsilon}(z), \\ c_0 & \text{if } x > d + V_b^{\varepsilon}(z), \end{cases} \quad (6.2) $$

with $c_0 > c_1$, and where $V_s^{\varepsilon}$ and $V_b^{\varepsilon}$ model respectively the free surface and the bottom topography (see Fig. 6.2 for an illustration). In our context, the free surface models the heave produced by the wind speed, and the uneven bottom topography models a sandy seabed with variations induced by water currents [126]. The forcing term $F$ is defined as:

$$ \mathbf{F}(t, z, x) = f(t) \Psi(x) \delta(z - L_S) \mathbf{e}_z, $$

where $\mathbf{e}_z$ stands for the unit vector pointing in the $z$-direction. This term models a source located in the plane $z = L_S$, emitting a signal $f(t)$ in the $z$-direction with transverse profile $\Psi$ supported in $(0, d)$. Furthermore, due to the continuity of the pressure field at the free surface $x = V_s^{\varepsilon}(z)$, which is tantamount to neglect the surface tension, the wave equation (6.1) is complemented by the following boundary conditions:

$$ p(t, z, V_s^{\varepsilon}(z)) = p_{\text{surface}} \quad (t, z) \in \mathbb{R} \times \mathbb{R}, $$

where $p_{\text{surface}}$ represents the atmospheric pressure. However, one can assume without loss of generality that:

$$ p_{\text{surface}} = 0, $$
<figure><img src="image_31.png" /><figcaption>Fig. 6.2: Illustration of the waveguide model. In this figure d represents the mean ocean depth, and the random fluctuations are given by the graphs of $V_s^\epsilon$ at the free surface, and $d+V_b^\epsilon$ at the bottom of the waveguide model. We illustrate the source term F emitting a wave $f(t)$ in the z-direction. $L/\epsilon^2$ characterizes the size of the section in which the random fluctuations occur. The right arrow (with a 0) pointing to the left at $z = L/\epsilon^2$ indicates that no wave is coming from the right at the end of the random section.</figcaption></figure>

which leads to Dirichlet boundary conditions at the free surface, and corresponding to a pressure-release condition. The central wavelength $\lambda$ of the source term is assumed to be of order one, and we consider a propagation distance given by $L/\epsilon^2$ to place ourselves in a high-frequency regime.

Waveguide models with a Pekeris profile (6.2) have been extensively studied for half a century [174], and has been widely used to model an ocean with a constant propagation speed profile (like in Figure 6.1). Such conditions can be found during the winter in Earth's mid latitudes and in water shallower than about 30 meters [138]. This profile is convenient for mathematical treatments, but it underestimates the real complexity of the medium. Nevertheless, the analysis carried out in this paper can be extended to more general propagation speed profiles. This model can also be used for electromagnetic waveguides, such as dielectric slabs or optical fibers with randomly perturbed boundaries [145, 149].

In (6.2), the perturbations of the waveguide geometry $V_s^\epsilon$ and $V_b^\epsilon$ are defined by

$$V_s^\epsilon(z) = \epsilon V_s(z) f_s(\epsilon^2 z) \quad \text{and} \quad V_b^\epsilon(z) = \epsilon V_b(z) f_b(\epsilon^2 z).$$

Here, $V_s$ and $V_b$ are assumed to be two independent mean-zero stationary bounded stochastic processes. The two functions $f_s$ and $f_b$, which specify the locations of the waveguide fluctuations, are assumed to be smooth with their supports included within the interval $(0, L)$. Note that the scaling $\epsilon \ll 1$, in front of $V_s$ and $V_b$, corresponds to small fluctuations of the waveguide geometry, and place ourselves in the weak scattering regime. In fact, wind speeds at approximately 5m/s result in a standard deviation of the surface roughness at about 0.1m. A similar order of magnitude applies to the roughness of sandy bottoms (see [126, Chapter 2] and [137]). Therefore, we have to wait for long propagation distances (of order $1/\epsilon^2$) to observe significant cumulative stochastic effects on the pressure wave. Let us note that the two random processes $V_s$ and $V_b$ are assumed to be independent. This assumption is justified because the surface and bottom standard deviations are small compared to the ocean depth, and then the bottom topography does not induce any ripples on the free surface. Power spectral densities for $V_s$ and $V_b$ can be found in the literature [126, Section 2.4 and Section 2.9]

$$I_s(u) := \int R_s(z)e^{iuz} dz = \frac{C_s}{|u|^{\alpha_s}} e^{-C'_s/u^{\mu_s}} \quad \text{and} \quad I_b(u) := \int R_b(z)e^{iuz} dz = \frac{C_b}{4\pi(u^2+1)^{\alpha_b/2}}, \quad (6.3)$$

with $C_s, C'_s, C_b > 0$, and

$$R_s(z) := \mathbb{E}[V_s(z+z_0)V_s(z_0)] \quad \text{and} \quad R_b(z) := \mathbb{E}[V_b(z+z_0)V_b(z_0)].$$

For instance, the cases $\alpha_s = 4, \alpha_s = 6$, and $\mu_s = 2$, $I_s$ is referred to as the Pierson-Neumann spectra, and for $\alpha_s = 5$ and $\mu_s = 4$ to the Pierson-Moscovitz spectra. Consequently, the trajectories of $V_s^\epsilon$
<figure><img src="image_32.png" /><figcaption>Fig. 6.3: Schematic representation of the change of coordinates from the unperturbed domain $D_0$ onto the perturbed domain $D_{V_s}$.</figcaption></figure>

and $V_b^\epsilon$ are at least $C^2$ with bounded derivatives, and therefore the surface fluctuations (the graph of $V_s^\epsilon$) produce no breaking wave.

To develop the asymptotic analysis, the random fluctuations $V_s$ and $V_b$ are assumed to be $\phi$-mixing processes: considering the $\sigma$-algebras

$$ \mathcal{F}_{0,z} := \sigma(V_s(u), V_b(u), \quad 0 \le u \le z) \quad \text{and} \quad \mathcal{F}_{z,\infty} := \sigma(V_s(u), V_b(u), \quad z \le u), $$

we have

$$ \sup_{\substack{z \ge 0 \\ A \in \mathcal{F}_{z+u, \infty} \\ B \in \mathcal{F}_{0,z}}} |\mathbb{P}(A|B) - \mathbb{P}(A)| \le \phi(u) \quad \text{with} \quad \phi \in L^1(\mathbb{R}) \cap L^2(\mathbb{R}). $$

This mixing property describes, through the function $\phi$, the loss of the statistical dependency between the past and the future of a random trajectories.

### 6.1.2 The conformal mapping

To investigate the asymptotic behavior of the acoustic pressure wave within the perturbed waveguide domain, we use a change of coordinates through a conformal transformation [179] from

$$ D_0 = \{u + iv \in \mathbb{C} : u > 0\} \quad \text{onto} \quad D_{V_s} = \{z + ix \in \mathbb{C} : x > V_s^\epsilon(z)\}. $$

This change of coordinates allows us to study the pressure wave $p(t, z, x)$ in a waveguide domain with a flat surface (see Fig. 6.3). In other words, one can focus our attention on $p_0(t, u, v)$ defined by

$$ p_0(t, u, v) = p(t, z(u, v), x(u, v)), $$

which represents the pressure wave $p(t, z, x)$ solution to (6.1) in a waveguide domain with a flat surface. Here, $z(u, v) + ix(u, v) = \Phi(u + iv)$ stands for the change of coordinates from $D_0$ onto $D_{V_s}$. In the new set of coordinates, the acoustic pressure wave $p_0$ satisfies Dirichlet boundary conditions at the flat surface of the waveguide,

$$ p_0(t, u, v = 0) = 0 \quad (t, u) \in \mathbb{R} \times \mathbb{R}. $$

The original pressure wave $p(t, z, x)$ can be obtained from $p_0(t, u, v)$ by inverting the change of coordinates,

$$ p(t, z, x) = p_0(t, \text{Re}(\Phi^{-1}(z + ix)), \text{Im}(\Phi^{-1}(z + ix))). $$

The main effect of this change of coordinates is to transfer the random fluctuations of the free surface to random fluctuations into the interior and at the bottom of the waveguide (see Fig. 6.4).
<figure><img src="image_33.png" /><figcaption>Fig. 6.4: Illustration of the conformal transformation applied to a waveguide with a free surface and uneven bottom topography. The two top figures depict waveguides with a free surface. In these figures the vertical dashed lines represent unperturbed medium parameters. In the top-left picture, the bottom is unperturbed, that is $V_b^\epsilon \equiv 0$, while in the top-right picture, the bottom of the waveguide is perturbed. The bottom picture illustrates the waveguide after the conformal transformation. In this figure, the perturbed vertical dashed lines mean perturbed medium parameters. Moreover, if the bottom was already perturbed before the conformal transformation, it is now also perturbed by the waveguide transformation.</figcaption></figure>

In the time-Fourier domain $p_0$ can be decomposed over a basis of modes,

$$ \hat{p}_0(\omega, u, v) = \underbrace{\sum_{j=1}^{N(\omega)} \hat{p}_j(\omega, u)\phi_j(\omega, v)}_{\text{propagating modes}} + \underbrace{\int_0^{k_{\omega}^2} \hat{p}_{\gamma}(\omega, u)\phi_{\gamma}(\omega, v)d\gamma}_{\text{radiating modes}} + \underbrace{\int_{-\infty}^0 \hat{p}_{\gamma}(\omega, u)\phi_{\gamma}(\omega, v)d\gamma}_{\text{evanescent modes}} \quad (6.4) $$

corresponding to the spectral decomposition of the (unperturbed) Pekeris operator defined by

$$ R_{\omega} = \frac{d^2}{dv^2} + k_{\omega}^2 n^2(v) \quad \text{with} \quad n(v) = \begin{cases} n_1 = c_0/c_1 > 1 & \text{if } v \in (0, d), \\ 1 & \text{if } v \ge d, \end{cases} $$

where $k_\omega = \omega/c_0$. The $(\phi_j)_{j \in \{1, \dots, N(\omega)\}}$ and $(\phi_\gamma)_{\gamma \in (-\infty, k^2(\omega))}$ are the (eigen)elements of the spectral measure associated to $R_\omega$, corresponding respectively to its discrete spectrum and continuous spectrum. The coefficients $\hat{p}_j$ and $\hat{p}_\gamma$ represent the amplitudes associated to each mode (propagating, radiative and evanescent modes).

The waveguide fluctuations result in a coupling between all these modes. In (6.4), the only modes that can propagate over long distances along the waveguide are the propagating modes. All the other modes play only a secondary role in the analysis.

### 6.1.3 Main result

Neglecting the backscattering effects, an asymptotic analysis can be applied based on a approximation-diffusion theorem. The effective mode-coupling mechanism can be described in term of a diffusion process, so that the expected pressure field of $\hat{p}_0(\omega, u, v)$, at the end of the perturbed section $u = L/\epsilon^2$, can be approximated for $\epsilon \ll 1$ by

$$ \mathbb{E}[\hat{p}_0(\omega, u = \frac{L}{\epsilon}, v)] \approx \sum_{j=1}^{N(\omega)} \frac{\hat{a}_{j,0}(\omega)}{\sqrt{\beta_j(\omega)}} e^{-T_j(\omega, L)+i\mathcal{P}_j(\omega, L)} e^{i\beta_j(\omega)L/\epsilon^2} \phi_j(\omega, v). $$

Here, the $\beta_j(\omega)$'s correspond to modal wavenumbers as introduced in Sect. 3.1, and the $\hat{a}_{j,0}(\omega)$'s represent the initial amplitudes of the propagating modes generated by the source term and entering
the random section. It turns out that the coherent transmitted wavefield at $u = L/\epsilon^2$ is composed of a superposition of propagating modes, where each mode is affected by an effective frequency-dependent attenuation and phase modulation through respectively $T_j(\omega, L)$ and $P_j(\omega, L)$. In particular, the attenuation coefficients are of the form

$$T_j(\omega, L) = C_j^s(\omega) \int_0^L f_s^2(z) dz + C_j^b(\omega) \int_0^L f_b^2(z) dz,$$

with $C_j^s(\omega), C_j^b(\omega) > 0$, accounting for the modal cumulative net scattering effects of the perturbed waveguide geometry.

Due to the exponential decay of the evanescent mode amplitudes w.r.t the propagation distance, these modes do not play a direct role in the transmitted wavefield. Nonetheless, they induce an indirect influence on the propagating modes through a mode-dependent phase modulation encoded in $P_j(\omega, L)$. The radiating modes do not play any direct significant role neither. Through the asymptotic analysis and the continuous nature of the associated spectrum, it can be shown that these modes are not affected by the perturbations of the propagation medium. Their contributions are therefore of the form

$$\int_0^{k_\omega^2} \frac{a_{\gamma,0}(\omega)}{\gamma^{1/4}} e^{i\sqrt{\gamma}L/\epsilon^2} \phi_\gamma(\omega, v) d\gamma,$$

involving a rapid phase term which cancels out the integral in the limit $\epsilon \to 0$. Nevertheless, they also play an indirect role in the attenuation coefficient $T_j(\omega, L)$ as described in the following result.

**Theorem 6.1.1** In the limit of a large number of propagating modes $k_\omega \gg 1$, knowing that $N(\omega) = [k_\omega d\theta/(c_1\pi)]$, we have for $\alpha \in \{s, b\}$:

* for $j \ll N(\omega)^{1/2}$,

$$C_j^\alpha(\omega) \underset{k_\omega \gg 1}{\sim} K_{1,\alpha} k_\omega^{3/2} \frac{j^2}{N(\omega)^2} \int_0^\infty \sqrt{u} I_\alpha(u) du;$$

* for $j \sim N(\omega)^{1/2}$,

$$C_j^\alpha(\omega) \underset{k_\omega \gg 1}{\sim} K_{2,\alpha} k_\omega^{3/2} \frac{j^2}{N(\omega)^2} \int_0^\infty \sqrt{u} I_\alpha \left( u - \frac{j^2 \pi \theta}{2N(\omega)d} \right) du;$$

* for $N(\omega)^{1/2} \ll j \lesssim \nu N(\omega)$, with $\nu \in (0, 1)$,

$$C_j^\alpha(\omega) \underset{k_\omega \gg 1}{\sim} K_{3,\alpha} k_\omega^2 \frac{j^3}{N(\omega)^3} \int_0^\infty I_\alpha(u) du;$$

* for $j \sim N(\omega)$,

$$K_{4,\alpha} k_\omega^{5/2} \le C_j^\alpha(\omega) \le K_{5,\alpha} k_\omega^{5/2};$$

where all the $K_{j,\alpha}$ are positive constants, the $I_\alpha$'s are given by (6.3), and $\theta = \sqrt{1 - 1/n_1^2}$.

This result provides the order of magnitude of the mode-dependent and frequency-dependent attenuation $T_j(\omega, L)$, which is responsible for the effective attenuation of the forward propagating-mode amplitudes. Let us make two remarks. First, the surface and bottom fluctuations affect the mode amplitudes mainly in a similar manner. Secondly, higher-order modes experience more significant attenuation as they interact more intensely with the random boundaries, leading to increased scattering. In fact, the higher the mode is the more it bounces on the random boundaries, and then the more it is scattered. While low-order modes have their attenuations mainly produced by the coupling between the propagating modes, high-order modes ($j \sim N(\omega)$) experience stronger attenuations due their coupling with the radiating modes. The reason is that only the highest order modes can efficiently couple with the radiating modes, resulting in a stronger attenuation as the radiating modes allow strong losses in the bottom of the waveguide [91]. This kind of result has also been reported in [7] for waveguides with a bounded cross-section, which do not support radiative modes and, therefore, does not exhibit these additional effects.
<figure><img src="image_34.png" /><figcaption>Fig. 6.5: Illustration of the basic physical setup: A source illuminates a rough surface producing reflected (in red) and transmitted waves (in green). Both waves exhibit a specular and diffusive components. The diffusive component, the speckle, corresponds to scattering of the incident wave by the rough interface.</figcaption></figure>

## 6.2 Reflection and transmission problems of high-frequency waves through a randomly perturbed interface: generalized Snell's laws

Despite the huge importance of this problem and the wide range of applications, there exists a limited number of rigorous results characterizing the wavefield, whether it is transmitted through or reflected from a rough interface [6, 13, 89, 15]. Many works in the physical literature consider a random interface problem in a perturbative situation, where the effects of small interface variations yield relatively small corrections in the transmitted and reflected wavefields. In this section, the interface fluctuations are modeled through a random field, enabling us to describe the transmitted and reflected fields, even though they are modified to leading order. It is worth mentioning that the problem of reflection from or transmission through a rough interface requires specialized techniques that differ from those associated with bulk propagation, that is the approximation-diffusion approach [70]. The basic configuration we consider here is illustrated in Fig. 6.5: a beam like wavefield (under the paraxial scaling) illuminates a rough surface with constant medium parameters, that are different above and below the interface.

The incident wave is typically decomposed into a reflected wave and a transmitted wave. These wavefields are further decomposed into two main components: the specular component and the speckle, also called incoherent or diffuse component, which results from scattering. Two central questions are then:

* i How does the fluctuations of the interface modify the specular components?

* ii How can we characterize the speckle components of the wavefields, and what parameters influence their relative magnitudes and spatial supports?

It turns out that the central scaling ratio that distinguishes various canonical scattering situations is the correlation length (or characteristic scale) of interface variations relative to the beam width. This is particularly relevant in the critical scaling scenario considered here, where the amplitude of the interface fluctuations is of the order of the wavelength. In particular, if the interface fluctuations occur on the scale of the probing cone width, the reflected and transmitted specular cones exhibit random arrival time properties. However, their main wave energies remain confined within the specular cones as the interface fluctuations does not induce strong coupling in between modes with different lateral wave vectors (traveling in relative oblique direction). On the other hand, when the interface fluctuations occur on a smaller scale relative to the width of the probing cone, such a coupling takes place. This generates homogenized specular reflected and transmitted cones with frequency-dependAttenuations related the interface elevation statistics. The missing energy from these effective specular cones have been converted into wide speckle cones that carry a total energy of order one. The radius
<figure><img src="image_35.png" /><figcaption>Fig. 6.6: Illustration of the physical setup. The plan $z = 0$ contains the source location, while the rough interface between $D_0$ and $D_1$ occurs around $z = z_{int}$. The reflected wave is observed at $z = 0$, and the transmitted wave is observed at $z = z_{tr}$.</figcaption></figure>

of these large cones are characterized by the ratio of the wavelength with the correlation length of the interface variations. In our situation, the roughness of the interface is not strong enough to generate the enhanced backscattering effects [121]. The scattering operator we obtain is similar to the one obtained in [183] under the Born (single scattering) approximation, and involves the statistics of the interface relative elevation.

There is a large literature addressing the important interface scattering problem. Most of the literature deals with physically motivated expansions, such as perturbative approaches or Kirchhoff-type approximations, which impose stringent conditions on the applicable scaling regimes [25, 57, 189]. Additionally, integral equation formulations are commonly employed, especially in characterizing near-field scattering [148, 187]. Sophisticated numerical methods have been developed [56, 184] to understand wave scattering from complicated geometries. Here, we provide statistical representations of such geometries and related scattering effects, as well as a statistical description of the speckles in term of a Gaussian random fields.

### 6.2.1 The physical model

**The wave equation** In this section, a three-dimensional linear wave propagation modeled is considered through the scalar wave equation:

$$ \Delta u - \frac{1}{c^2(\mathbf{x}, z)} \partial_{tt}^2 u = \nabla \cdot F(t, \mathbf{x}, z) \quad (t, \mathbf{x}, z) \in \mathbb{R} \times \mathbb{R}^2 \times \mathbb{R}, \qquad (6.5) $$

equipped with null initial conditions

$$ u(t = 0, \mathbf{x}, z) = \partial_t u(t = 0, \mathbf{x}, z) = 0 \quad (\mathbf{x}, z) \in \mathbb{R}^2 \times \mathbb{R}, $$

and continuity conditions at the interface that are specified below. The coordinate $z$ represents the main propagation axis while $\mathbf{x}$ represents the transverse directions, $F$ corresponds to the source term, and the Laplacian operator $\Delta = \Delta_{\perp} + \partial_{zz}^2$ acts on all spatial variables $\mathbf{x}$ and $z$. The propagation medium consists of two homogeneous subdomains separated by a randomly perturbed interface around $z = z_{int}$:

$$ D^0 := \{(\mathbf{x}, z) \in \mathbb{R}^2 \times \mathbb{R} \text{ s.t. } z < z_{int} + \sigma V(\mathbf{x}/l_c)\} \qquad (6.6) $$

and

$$ D^1 := \{(\mathbf{x}, z) \in \mathbb{R}^2 \times \mathbb{R} \text{ s.t. } z > z_{int} + \sigma V(\mathbf{x}/l_c)\}. \qquad (6.7) $$

We refer to Fig. 6.6 for an illustration of the physical setup. The term $V$ corresponds to a mean-zero random field with second order derivatives (for technical reasons), modeling the interface variations and characterizing its roughness. Away from the interface the velocity field is given by:

$$ c(\mathbf{x}, z) := \begin{cases} c_0 & \text{if } (\mathbf{x}, z) \in D^0, \\ c_1 & \text{if } (\mathbf{x}, z) \in D^1. \end{cases} $$
The forcing term

$$F(t, \mathbf{x}, z) := \Psi\left(\frac{t - \mathbf{k}_0 \cdot \mathbf{x}}{\lambda}, \frac{\mathbf{x}}{r_0}\right) \delta(z) \mathbf{e}_z,$$

where $\mathbf{e}_z$ denotes the unit vector pointing in the z-direction, models a source located in the plane $z = 0$ and emitting a quasi plan wave $\Psi$ in the spatial direction $(\mathbf{k}_0, (c_0^{-2} - |\mathbf{k}_0|^2)^{1/2})$ towards the random interface (requiring $|\mathbf{k}_0| < c_0^{-1}$). In a dimensionless situation, $\lambda$ represents the wavelength and $r_0$ the spatial beam width of the source term.

The above wave equation (6.5) only requires the two following continuity relations across the randomly perturbed interface:

$$u(z = z_{\text{int}}(\mathbf{x})^{+}) = u(z = z_{\text{int}}(\mathbf{x})^{-}) \quad \text{and} \quad \partial_z u(z = z_{\text{int}}(\mathbf{x})^{+}) = \partial_z u(z = z_{\text{int}}(\mathbf{x})^{-}),$$

where

$$z_{\text{int}}(\mathbf{x}) := z_{\text{int}} + \sigma V(\mathbf{x}/l_c).$$

Finally, no wave is assumed to come from above the source location nor below the interface. The only waves present in the system are generated by the source term.

**The parameter scaling** Here, the distance from the source to the interface is denoted by $z_{\text{int}}$, and the distance from the source to the plane where the transmitted wave is recorded by $z_{\text{tr}}$. Both of these distances are assumed to be of order $L \sim 1$ (the reference scale). The dimensionless parameter

$$\varepsilon := \frac{\lambda}{L} \ll 1,$$

is assumed to be small to provide a high-frequency regime. Moreover, a paraxial (or parabolic) scaling is considered by enforcing

$$\frac{r_0^2}{\lambda} \sim L,$$

which corresponds to a Rayleigh length of order the typical propagation distance. As already mentioned in this manuscript, in the homogeneous case, the Rayleigh length $L_R = \pi r_0^2 / \lambda$ corresponds to the distance from the beam waist to the place where the area of its cross-section is doubled. The interface fluctuations are here *critically scaled* in the sense that the magnitude of the interface fluctuations and the central wavelength are of the same order,

$$\sigma \sim \lambda.$$

Finally, $l_c$ denotes the correlation length of the interface fluctuations, which ranges from the central wavelength to the beam width:

$$\lambda \lesssim l_c \lesssim r_0.$$

Note that if $r_0 \ll l_c$ the situation is trivial and corresponds essentially to a planar interface. The ratio $\lambda/l_c$ characterizes how rough the surface is for the incident wave, and the cone radius for the reflected and transmitted speckle components. For simplicity, we set

$$L = L_R \sim 1, \quad \lambda = \varepsilon, \quad r_0 = \sqrt{\varepsilon}, \quad \sigma = \varepsilon, \quad \text{and} \quad l_c = \varepsilon^{\gamma}, \quad (6.8)$$

with

$$\gamma \in [1/2, 1].$$

In the forthcoming presentation, a distinction is made between $\gamma = 1/2$ (that is $l_c \sim r_0$) and $1/2 < \gamma \le 1$ (corresponding to $\lambda \lesssim l_c \ll r_0$). These two scenarios result in fundamentally different behaviors for the reflected and transmitted signals. While the former induces only random specular components (without speckle), the latter, due to fast interface oscillations compared to the beam radius ($l_c \ll r_0$), induces homogenization effects for the specularly reflected and the transmitted wave components. In this case, random speckles over cones larger than the ones covered by the specular and transmitted components are observed.
**Random fluctuations** In (6.6) and (6.7), the random fluctuations of the interface delimiting the two subspaces are defined by a mean-zero stationary random field $V$. These random fluctuations are assumed to satisfy *mixing* properties describing the loss of statistical dependency for $V$ over the interface. In our 2D context, the *mixing* property is defined as follows. We introduce

$$ \alpha(r) := \sup_{\substack{S,S' \subset \mathbb{R}^2 \\ d(S,S') > r}} \sup_{\substack{A \in \sigma(V(\mathbf{x}), \mathbf{x} \in S) \\ B \in \sigma(V(\mathbf{x}), \mathbf{x} \in S')}} |\mathbb{P}(A \cap B) - \mathbb{P}(A)\mathbb{P}(B)|, $$

where

$$ d(S, S') = \inf_{\substack{s \in S \\ s' \in S'}} |s - s'| $$

is the distance between two nonempty subsets $S$ and $S'$, and $\sigma(V(\mathbf{x}), \mathbf{x} \in S)$ is the $\sigma$-field generated by the family $V_{|S} := (V(\mathbf{x}))_{\mathbf{x} \in S}$. Roughly speaking, the value $\alpha(r)$ quantifies the degree of statistical dependency for the random field $V$ over pair of regions at distance at least $r$. The $\alpha$-mixing property consists in assuming

$$ \alpha(r) \to 0 \quad \text{as} \quad r \to \infty, $$

giving a vanishing of the statistical dependency between $V_{|S}$ and $V_{|S'}$ as the distance between $S$ and $S'$ tends to infinity.

### 6.2.2 Refection and transmission for an unperturbed interface

In the context of an unperturbed (flat) interface, the asymptotic specular reflected wave is given by

$$ \begin{aligned} U^{\text{ref}}(s, \mathbf{y}) &:= \lim_{\varepsilon \to 0} u^{\varepsilon}(t_{\text{obs,ref}}^{\varepsilon}(\mathbf{y}) + \varepsilon s, \mathbf{x}_{\text{obs,ref}} + \sqrt{\varepsilon} \mathbf{y}, z=0) \quad (6.9) \\ &= \frac{\mathcal{R}}{2(2\pi)^3} \int e^{-i\omega(s-\mathbf{y}\cdot\mathbf{q})} \hat{U}_0(\omega, \mathbf{q}, 2z_{\text{int}}) \hat{\Psi}(\omega, \mathbf{q}) \omega^2 d\omega d\mathbf{q}. \end{aligned} $$

Here, $u^\varepsilon$ corresponds to the scaled solution of (6.5) according to (6.8). Therefore, $U^{\text{ref}}$ corresponds to the asymptotic wave front observed at position $\mathbf{x}_{\text{obs,ref}} + \sqrt{\varepsilon} \mathbf{y}$ (depending on the source frame $\sqrt{\varepsilon} \mathbf{y}$) and time $t_{\text{obs,ref}}^{\varepsilon}(\mathbf{y}) + \varepsilon s$ (depending of the pulse duration $\varepsilon s$), where

$$ t_{\text{obs,ref}}^{\varepsilon}(\mathbf{y}) := t_{\text{obs,ref}} + \sqrt{\varepsilon} \mathbf{k}_0 \cdot \mathbf{y}, \quad \text{with} \quad t_{\text{obs,ref}} := \frac{2z_{\text{int}}}{c_0^2 \lambda_0}, $$

and

$$ \mathbf{x}_{\text{obs,ref}} := \frac{2\mathbf{k}_0 z_{\text{int}}}{\lambda_0}. $$

Here, $\lambda_0$ is defined by

$$ \lambda_0 = \frac{\sqrt{1 - c_0^2 |\mathbf{k}_0|^2}}{c_0}. $$

The travel time is influenced here by the offset $\mathbf{y}$ as the beam width $\sqrt{\varepsilon}$ is large compared to the pulse duration $\varepsilon$, and that the wave front travels obliquely relative to the vertical direction. Note that $\mathbf{x}_{\text{obs,ref}}$ is twice the lateral position

$$ \mathbf{x}_{\text{int}} := \frac{\mathbf{k}_0 z_{\text{int}}}{\lambda_0}, \quad (6.10) $$

where the incident pulse hits the interface. Therefore, the standard reflexion relation reads

$$ \theta_{\text{inc}} = \theta_{\text{ref}}^0 := \arctan \left( \frac{|\mathbf{k}_0|}{\lambda_0} \right), \quad (6.11) $$

where the incident and reflected angles are equal. We refer to Fig. 6.7 for an illustration of the geometrical properties of the reflection.

Regarding the reflected wave profile, the reflection coefficient $\mathcal{R}$ is given by

$$ \mathcal{R} := \frac{\lambda_0 - \lambda_1}{\lambda_0 + \lambda_1}, $$
<figure><img src="image_36.png" /><figcaption>Fig. 6.7: Illustration of the reflection and transmission in the context of a flat interface at $z = z_{int}$. The source is located at $\mathbf{x} = 0$ in the plan $z = 0$, and the emitted wave hits the interface at $\mathbf{x} = \mathbf{x}_{int}$, in the plan $z = z_{int}$, with an incident angle $\theta_{inc}$. The reflected wave, with angle $\theta_{ref}^0$, is observed in the plan $z = 0$ at $\mathbf{x} = \mathbf{x}_{obs,ref}$, while the transmitted wave, with angle $\theta_{tr}^0$, is observed in the plan $z = z_{tr}$ at $\mathbf{x} = \mathbf{x}_{obs,tr}$.</figcaption></figure>

where $\lambda_1$ is defined as $\lambda_0$ by just replacing $c_0$ by $c_1$, and

$$ \hat{\mathcal{U}}_0(\omega, \mathbf{q}, z) := e^{-i\omega z c_0 \mathbf{q}^T A_0 \mathbf{q}}, \qquad (6.12) $$

where $A_0$ is defined by

$$ A_0 := \frac{1}{2(1 - c_0^2 |\mathbf{k}_0|^2)^{3/2}} \begin{pmatrix} 1 - c_0^2 \mathbf{k}_{0,2}^2 & c_0^2 \mathbf{k}_{0,1} \mathbf{k}_{0,2} \\ c_0^2 \mathbf{k}_{0,1} \mathbf{k}_{0,2} & 1 - c_0^2 \mathbf{k}_{0,1}^2 \end{pmatrix} = \frac{1}{2c_0^3 \lambda_0^3} (\mathbf{I}_2 - c_0^2 \mathbf{k}_0^\perp (\mathbf{k}_0^\perp)^T). \qquad (6.13) $$

The term $\hat{\mathcal{U}}_0$ leads to the following homogeneous semi-group

$$ \check{\mathcal{U}}_0(\omega, \mathbf{y}, z) := \frac{\omega^2}{(2\pi)^2} \int e^{i\omega\mathbf{y}\cdot\mathbf{q}} e^{-i\omega z c_0 \mathbf{q}^T A_0 \mathbf{q}} d\mathbf{q}, \qquad (6.14) $$

satisfying the Schrödinger type equation

$$ i\partial_z \check{\mathcal{U}}_0(\omega, \mathbf{y}, z) + \frac{1}{k_{0,\omega}} \nabla_{\mathbf{y}} \cdot (A_0 \nabla_{\mathbf{y}} \check{\mathcal{U}}_0)(\omega, \mathbf{y}, z) = 0 \quad z > 0, \qquad (6.15) $$

where

$$ \check{\mathcal{U}}_0(\omega, \mathbf{y}, z = 0) = \delta(\mathbf{y}), \quad \text{and} \quad k_{0,\omega} = \omega/c_0, $$

which is characteristic of the paraxial approximation. Due to the initial lateral direction $\mathbf{k}_0$, we do not have the standard Schrödinger equation with a Laplacian term. Nevertheless, the standard Laplacian can be recovered by sending this direction to zero. Under this paraxial wave model, the limit (6.9) corresponds to the travel of the emitted pulse from the source location to the interface, and back from the interface to the source location (yielding the $z = 2z_{int}$).

The transmitted specular wave is given by

$$ \begin{aligned} U^{tr}(s, \mathbf{y}) &:= \lim_{\epsilon \to 0} u^{\epsilon} (t_{obs,tr}^{\epsilon}(\mathbf{y}) + \epsilon s, \mathbf{x}_{obs,tr} + \sqrt{\epsilon} \mathbf{y}, z = z_{tr}) \qquad (6.16) \\ &= \frac{\mathcal{T}}{2(2\pi)^3} \sqrt{\frac{\lambda_0}{\lambda_1}} \int e^{-i\omega(s-\mathbf{y}\cdot\mathbf{q})} \hat{\mathcal{U}}_1(\omega, \mathbf{q}, z_{tr} - z_{int}) \hat{\mathcal{U}}_0(\omega, \mathbf{q}, z_{int}) \hat{\Psi}(\omega, \mathbf{q}) \omega^2 d\omega d\mathbf{q}, \end{aligned} $$

with this time

$$ \mathbf{x}_{obs,tr} := \mathbf{x}_{int} + \mathbf{x}_{tr} := \mathbf{x}_{int} + \frac{\mathbf{k}_0(z_{tr} - z_{int})}{\lambda_1}, $$

where $\mathbf{x}_{int}$ is given by (6.10), and

$$ t_{obs,tr}^{\epsilon}(\mathbf{y}) := t_{obs,tr} + \sqrt{\epsilon}\mathbf{y}, \quad \text{with} \quad t_{obs,tr} := \frac{z_{int}}{c_0^2\lambda_0} + \frac{z_{tr} - z_{int}}{c_1^2\lambda_1}. $$
Here, $U^{tr}$ corresponds to the wave front observed in the frame of the source term at position $\mathbf{x}_{obs,tr}$, and time $t_{obs,tr}^{\varepsilon}(\mathbf{y})$. The time $t_{obs,tr}$, for the pulse to reach the plan $z = z_{tr}$, is the sum of the time to reach the interface and the time to go from the interface to the plan $z = z_{tr}$. The lateral position $\mathbf{x}_{obs,tr}$, where the pulse hits the plan $z = z_{tr}$, corresponds to a sum of vectors. The first vector represents the lateral position where the pulse hits the interface, and the second one the additional lateral displacement once the pulse goes through the interface and reaches the plan $z = z_{tr}$. We refer to Fig. 6.7 for an illustration for the geometrical properties of transmission. From this formulation, one has the usual Snell's law:

$$ \frac{\sin(\theta_{inc})}{c_0} = \frac{\sin(\theta_{tr}^0)}{c_1}, $$

with $\theta_{inc}$ defined by (6.11), and

$$ \theta_{tr}^0 := \arctan \left( \frac{|\mathbf{k}_0|}{\lambda_1} \right). $$

Regarding the transmitted pulse profile, the transmission coefficient is defined by

$$ \mathcal{T} := \frac{2\sqrt{\lambda_0 \lambda_1}}{\lambda_0 + \lambda_1}, $$

so that the conservation relation $\mathcal{R}^2 + \mathcal{T}^2 = 1$ holds, and $\hat{\mathcal{U}}_0$ is defined by (6.12). The term $\hat{\mathcal{U}}_1$ is defined in the same way as $\hat{\mathcal{U}}_0$:

$$ \hat{\mathcal{U}}_1(\omega, \mathbf{q}, z) := e^{-i\omega z c_1 \mathbf{q}^T A_1 \mathbf{q}}, \qquad (6.17) $$

where $A_1$ is defined as $A_0$ (6.13) by just changing $c_0$ by $c_1$. The term $\hat{\mathcal{U}}_1$, defined as (6.14) with $c_1$ and $A_1$ instead of $c_0$ and $A_0$, satisfies a corresponding paraxial wave equation (6.15). Therefore, the formulation (6.16) corresponds to the travel of the emitted pulse from the source location to the interface, and from the interface to the plane $z = z_{tr}$ under the paraxial approximation. The factor $\mathcal{T}$ accounts for the transmission at the interface.

### 6.2.3 Random specular components for $l_c \sim r_0$ ($\gamma = 1/2$)

From now on, the interface is randomly perturbed. In case of a correlation length of order the beam radius, the asymptotic reflected pulse front, as defined in (6.9), is given by the following result.

**Proposition 6.2.1** The family $(U^{\varepsilon,ref})_{\varepsilon}$ converges in law in $L^2(\mathbb{R} \times \mathbb{R}^2)$, as $\varepsilon \to 0$, to

$$ U^{ref}(s, \mathbf{y}) := \frac{\mathcal{R}}{2(2\pi)^5} \iiint \int e^{-i\omega(s-2\lambda_0\mathcal{V}(\mathbf{y}')-\mathbf{q}\cdot\mathbf{y})} e^{i\omega(\mathbf{q}'-\mathbf{q})\cdot\mathbf{y}'} \\ \times \hat{\mathcal{U}}_0(\omega, \mathbf{q}, z_{int}) \hat{\mathcal{U}}_0(\omega, \mathbf{q}', z_{int}) \hat{\Psi}(\omega, \mathbf{q}')\omega^4 d\omega d\mathbf{y}' d\mathbf{q}' d\mathbf{q}, $$

where $\mathcal{V}$ is a random field with the same law as $\mathcal{V}$.

The reflected wave front consists in a superposition of contribution observed at point $\mathbf{y}$. Each of these contributions results in waves emitted with initial lateral direction $\mathbf{q}'$, propagating according to $\hat{\mathcal{U}}_0(\mathbf{q}')$, and scattered into direction $\mathbf{q}$ at point $\mathbf{y}'$ before propagating back according to $\hat{\mathcal{U}}_0(\mathbf{q})$ (see Fig. 6.8 for an illustration). Note that a random time correction $2\lambda_0\mathcal{V}(\mathbf{y}')$ is observed for contributions scattered at point $\mathbf{y}'$. In $U^{ref}$, for a given frequency $\omega$, the scattering operator corresponds to

$$ K_{ref}(\omega, \mathbf{q}', \mathbf{q}) := \int e^{i\omega(\mathbf{q}'-\mathbf{q})\cdot\mathbf{y}'} e^{2i\omega\lambda_0\mathcal{V}(\mathbf{y}')} d\mathbf{y}'. \qquad (6.18) $$

A similar result holds for the transmitted wave front.

**Proposition 6.2.2** The family $(U^{\varepsilon,tr})_{\varepsilon}$ converges in law in $L^2(\mathbb{R} \times \mathbb{R}^2)$, as $\varepsilon \to 0$, to

$$ U^{tr}(s, \mathbf{y}) := \frac{\mathcal{T}}{2(2\pi)^5} \sqrt{\frac{\lambda_0}{\lambda_1}} \iiint \int e^{-i\omega(s-(\lambda_0-\lambda_1)\mathcal{V}(\mathbf{y}')-\mathbf{q}\cdot\mathbf{y})} e^{i\omega(\mathbf{q}'-\mathbf{q})\cdot\mathbf{y}'} \\ \times \hat{\mathcal{U}}_1(\omega, \mathbf{q}, z_{tr} - z_{int}) \hat{\mathcal{U}}_0(\omega, \mathbf{q}', z_{int}) \hat{\Psi}(\omega, \mathbf{q}')\omega^4 d\omega d\mathbf{y}' d\mathbf{q}' d\mathbf{q}, $$

where $\mathcal{V}$ is a random field with the same distribution as $\mathcal{V}$.
<figure><img src="image_37.png" /><figcaption>Fig. 6.8: Illustration of the contributions to the reflected and transmitted waves observed at points $\mathbf{x}_{\text{obs,ref}} + \sqrt{\varepsilon}\mathbf{y}$ and $\mathbf{x}_{\text{obs,tr}} + \sqrt{\varepsilon}\mathbf{y}$ respectively. The dashed light grey lines represent the specular cones.</figcaption></figure>

The transmitted wave front $U^{tr}$ consists as well in a superposition of diffracted contributions observed at point $\mathbf{y}$. Each of these contribution results in waves emitted with initial lateral direction $\mathbf{q}'$, propagating according to $\hat{U}_0$, scattered into direction $\mathbf{q}$ at point $\mathbf{y}'$ when going through the interface, and then propagates according to $\hat{U}_1$ until it reaches $z = z_{tr}$ (see Figure 6.8 for an illustration). For the transmitted wave, the scattering operator reads

$$ K_{tr}(\omega, \mathbf{q}', \mathbf{q}) := \int e^{i\omega(\mathbf{q}' - \mathbf{q}) \cdot \mathbf{y}'} e^{i\omega(\lambda_0 - \lambda_1)\nu(\mathbf{y}')} d\mathbf{y}'. \quad (6.19) $$

### 6.2.4 Effective specular components for $l_c \ll r_0$ ($\gamma > 1/2$)

The case $\gamma > 1/2$ corresponds to a correlation length smaller than the beam width. The incident wave is therefore strongly scattered and homogenization phenomena take place, as well as incoherent wave fluctuations. In this context the asymptotic behavior of the reflected wave front (6.9) is given by the following result.

**Proposition 6.2.3** The family $(U^{\varepsilon,ref})_{\varepsilon}$ converges in probability, as $\varepsilon \to 0$, in $S'(\mathbb{R} \times \mathbb{R}^2)$, the set of tempered distributions, to the deterministic pulse profile

$$ U^{\text{ref}}(s, \mathbf{y}) = \frac{\mathcal{R}}{2(2\pi)^3} \iiint \int e^{-i\omega(s-\mathbf{q}\cdot\mathbf{y})} \hat{U}_0(\omega, \mathbf{q}, 2z_{\text{int}}) \phi_V(2\omega\lambda_0) \hat{\Psi}(\omega, \mathbf{q})\omega^2 d\omega d\mathbf{q}. $$

Here,

$$ \phi_V(u) = \mathbb{E}[e^{iuV(0)}] \quad (6.20) $$

is the characteristic function of $V(0)$, and $\hat{U}_0$ is defined by (6.12).

In this result, for simplicity, the wave front does not have any dependence at the scale of the correlation length as the asymptotic profile would not depend on it. Let us make two remarks. First, the impact of the interface fluctuations on the specular reflected wave can be described as the reflection problem with a flat interface where the homogenized scattering property are captured through a convolved initial pulse profile:

$$ \frac{1}{2} \Phi *_s \Psi(s, \mathbf{y}), $$

where $*_s$ is the convolution operator acting on the s-variable, and

$$ \Phi(s) := \frac{1}{2\pi} \int e^{-i\omega s} \phi_V(2\omega\lambda_0) d\omega $$

corresponds to a scaled version of the probability density function of $V(0)$. Second, compared to the case $\gamma = 1/2$, the scattering operator (6.18) is now homogenized,

$$ \mathbb{E}[K_{\text{ref}}(\omega, \mathbf{q}', \mathbf{q})] = \int e^{i\omega(\mathbf{q}' - \mathbf{q}) \cdot \mathbf{y}'} \mathbb{E}[e^{2i\omega\lambda_0 \nu(\mathbf{y}')}] d\mathbf{y}' = \delta(\omega(\mathbf{q}' - \mathbf{q}))\phi_V(2\omega\lambda_0). $$
This homogenized version of the scattering operator acts as for a flat surface, it does not modify the incident direction $\mathbf{q}'$, and the effective contributions of the interface fluctuations produce only a low-pass filter in frequency through the characteristic function $\phi_V$.

The same strategy applies for the transmitted wave front, and we have the following result.

**Proposition 6.2.4** The family $(U^{\varepsilon, tr})_\varepsilon$ converges in probability, as $\varepsilon \to 0$, in $S'(\mathbb{R} \times \mathbb{R}^2)$ to the deterministic pulse profile:

$$U^{tr}(s, \mathbf{y}) = \frac{T}{2(2\pi)^5} \sqrt{\frac{\lambda_0}{\lambda_1}} \iiint \int e^{-i\omega(s-\mathbf{q}\cdot\mathbf{y})} \hat{U}_1(\omega, \mathbf{q}, z_{tr} - z_{int}) \hat{U}_0(\omega, \mathbf{q}, z_{int}) \\ \times \phi_V((\lambda_0 - \lambda_1)\omega) \hat{\Psi}(\omega, \mathbf{q})\omega^2 d\omega d\mathbf{x}'d\mathbf{q},$$

where $\phi_V$ is defined by (6.20), the $\hat{U}_j$'s by (6.12) and (6.17).

The profile $U^{tr}$ is similar to the one obtained for a flat surface (6.16), but again with a convolved initial condition accounting for the homogenization effects. The scattering operator (6.19) becomes here

$$\mathbb{E}[K_{tr}(\omega, \mathbf{q}', \mathbf{q})] = \int e^{i\omega(\mathbf{q}'-\mathbf{q})\cdot\mathbf{y}'} \mathbb{E}[e^{i\omega(\lambda_0-\lambda_1)\mathcal{V}(\mathbf{y}')}] d\mathbf{y}' = \delta(\omega(\mathbf{q}'-\mathbf{q}))\phi_V((\lambda_0-\lambda_1)\omega).$$

### 6.2.5 Incoherent wave fluctuations and generalized Snell's laws for $l_c \ll r_0$ ($\gamma > 1/2$)

We now aim at describing the incoherent wave fluctuations, resulting from diffraction, which account for the missing energy carried by the homogenized specular components.

The reflected speckle profile refers to the following reflected wavefield

$$S^{\varepsilon,ref}(\bar{s}, \bar{\mathbf{y}}, \mathbf{y}, \tilde{s}, \tilde{\mathbf{y}}) := \varepsilon^{-2(\gamma-1/2)} u^{\varepsilon,ref}(t_{\text{obs,ref}}^{\varepsilon}(\bar{s}, \bar{\mathbf{y}}, \mathbf{y}, \tilde{\mathbf{y}}) + \varepsilon \tilde{s}, \mathbf{x}_{\text{obs,ref}}^{\varepsilon}(\bar{\mathbf{y}}, \mathbf{y}) + \varepsilon^{\gamma} \tilde{\mathbf{y}}, z=0),$$

observed at position

$$\mathbf{x}_{\text{obs,ref}}^{\varepsilon}(\bar{\mathbf{y}}, \mathbf{y}) := \mathbf{x}_{\text{obs,ref}} + \varepsilon^{1-\gamma} \bar{\mathbf{y}} + \sqrt{\varepsilon} \mathbf{y}, \quad (6.21)$$

and time

$$t_{\text{obs,ref}}^{\varepsilon}(\bar{s}, \bar{\mathbf{y}}, \mathbf{y}, \tilde{\mathbf{y}}) := t_{\text{obs,ref}} + \varepsilon^{1-\gamma} \mathbf{k}_0 \cdot \bar{\mathbf{y}} + \sqrt{\varepsilon} \mathbf{k}_0 \cdot \mathbf{y} + \varepsilon^{\gamma} \mathbf{k}_0 \cdot \tilde{\mathbf{y}} + \varepsilon^{2(1-\gamma)} \bar{s}. \quad (6.22)$$

Here, the observation point and time exhibit extra terms in $\varepsilon^{1-\gamma}$ and $\varepsilon^{2(1-\gamma)}$ corresponding to the roughness parameter $\lambda/l_c$. These terms account for neighborhood and time windows larger than the ones of the specular components (see Fig. 6.9 for an illustration). They can even be of order one if $\lambda \sim l_c$ (that is $\gamma = 1$). The order of magnitude of the speckle can be understood as spreading a beam width of order $r_0$ over a two dimensional spatial window of order $r_0^2/l_c \sim \lambda/l_c$. Therefore, it is indeed expected to observe a speckle with amplitude of order

$$\left(\frac{r_0}{r_0^2/l_c}\right)^2 \sim \left(\frac{l_c}{r_0}\right)^2 \sim \varepsilon^{2(\gamma-1/2)}.$$

The following result describes the asymptotic behavior of the two-point statistics of the speckle profile:

$$C^{\varepsilon,ref}(\bar{s}, \bar{\mathbf{y}}, \mathbf{y}, \tilde{s}_1, \tilde{\mathbf{y}}_1, \tilde{s}_2, \tilde{\mathbf{y}}_2) := S^{\varepsilon,ref}(\bar{s}, \bar{\mathbf{y}}, \mathbf{y}, \tilde{s}_1, \tilde{\mathbf{y}}_1)S^{\varepsilon,ref}(\bar{s}, \bar{\mathbf{y}}, \mathbf{y}, \tilde{s}_2, \tilde{\mathbf{y}}_2). \quad (6.23)$$

**Proposition 6.2.5** We have

$$\lim_{\varepsilon \to 0} C^{\varepsilon,ref}(\bar{s}, \bar{\mathbf{y}}, \mathbf{y}, \tilde{s}_1, \tilde{\mathbf{y}}_1, \tilde{s}_2, \tilde{\mathbf{y}}_2) = \mathbf{C}^{\text{ref}}(\bar{s}, \bar{\mathbf{y}}, \tilde{s}_1 - \tilde{s}_2, \tilde{\mathbf{y}}_1 - \tilde{\mathbf{y}}_2),$$

in probability in $S'(\mathbb{R} \times \mathbb{R}^2 \times \mathbb{R}^2 \times \mathbb{R} \times \mathbb{R}^2 \times \mathbb{R} \times \mathbb{R}^2)$, where

$$\mathbf{C}^{\text{ref}}(\bar{s}, \bar{\mathbf{y}}, \tilde{s}, \tilde{\mathbf{y}}) := \frac{\mathcal{R}^2}{4(2\pi)^3} \iint e^{-i\omega(\tilde{s}-\mathbf{p}\cdot\tilde{\mathbf{y}})} \mathcal{A}(2\lambda_0, \omega, \mathbf{p}) |\hat{\Psi}|_2^2(\omega) \times \delta(\bar{s} - s_{\mathbf{p}}^{\text{ref}}) \delta(\tilde{\mathbf{y}} - \mathbf{y}_{\mathbf{p}}^{\text{ref}}) \omega^4 d\omega d\mathbf{p}, \quad (6.24)$$
<figure><img src="image_38.png" /><figcaption>Fig. 6.9: The left picture illustrates the zone where the reflected speckle profile can be observed (in light grey) compared to the specular reflected component (in grey). The right picture illustrates the spatial window over which the speckle profile is analyzed.</figcaption></figure>

with

$$ \mathcal{A}(v, \omega, \mathbf{p}) := \int \mathbb{E}[e^{iv(V(\mathbf{y}') - V(0))}] e^{-i\omega\mathbf{p}\cdot\mathbf{y}'} d\mathbf{y}', \quad (6.25) $$

and

$$ |\hat{\Psi}|_2^2(\omega) := \frac{1}{(2\pi)^2} \int |\hat{\Psi}(\omega, \mathbf{q})|^2 d\mathbf{q}. \quad (6.26) $$

We also have,

$$ \mathbf{y}_{\mathbf{p}}^{\text{ref}} := 2z_{\text{int}}c_0 A_0 \mathbf{p} \quad \text{and} \quad s_{\mathbf{p}}^{\text{ref}} := \mathbf{p} \cdot \mathbf{y}_{\mathbf{p}}^{\text{ref}} / 2 = z_{\text{int}}c_0 \mathbf{p}^T A_0 \mathbf{p} \ge 0. $$

The asymptotic correlation function exhibits stationarity w.r.t. the variables in $\tilde{s}, \tilde{\mathbf{y}}$ (evolving respectively at the scale of the pulse width and correlation length), but does not depend on the variable $\mathbf{y}$ (corresponding to variations at the beam width scale). This means that the two-point statistics in $\tilde{s}$ and $\tilde{\mathbf{y}}$ are identical over regions of order the beam width. In a similar manner, the intensity carried by the speckle profile satisfies

$$ \lim_{\epsilon \to 0} |S^{\epsilon, \text{ref}}(\tilde{s}, \tilde{\mathbf{y}}, \mathbf{y}, \tilde{s}, \tilde{\mathbf{y}})|^2 = \mathbf{C}^{\text{ref}}(\tilde{s}, \tilde{\mathbf{y}}, 0, 0), $$

which does not depend on the variable $\mathbf{y}$ and the small scale fluctuations in $\tilde{s}$ and $\tilde{\mathbf{y}}$. The intensity is therefore uniform over these scales.

Regarding the self-averaging property of the correlation function, the speckle profile is here correlated around a given $(\tilde{s}, \tilde{\mathbf{y}}, \mathbf{y})$. Hence, looking at the second order moment of the correlation function in $S'$ corresponds to look at the fourth order moment of the speckle, but with two elements around some $(\tilde{s}_1, \tilde{\mathbf{y}}_1, \mathbf{y}_1)$ and two others around some $(\tilde{s}_2, \tilde{\mathbf{y}}_2, \mathbf{y}_2)$. These two observation points being far apart from each other (in time and space) the correlated speckles around each of these points become statistically independent:

$$ \mathbb{E}[C^{\epsilon, \text{ref}}(\tilde{s}_1, \tilde{\mathbf{y}}_1, \mathbf{y}_1)C^{\epsilon, \text{ref}}(\tilde{s}_2, \tilde{\mathbf{y}}_2, \mathbf{y}_2)] \approx \mathbb{E}[C^{\epsilon, \text{ref}}(\tilde{s}_1, \tilde{\mathbf{y}}_1, \mathbf{y}_1)]\mathbb{E}[C^{\epsilon, \text{ref}}(\tilde{s}_2, \tilde{\mathbf{y}}_2, \mathbf{y}_2)]. $$

The term $\mathcal{A}(\cdot, \omega, \mathbf{p})$ corresponds to the distribution of the scattered directions $\mathbf{p}$ produced by the random interface. These directions $\mathbf{p}$ are related to variations at the scale of the correlation length through the characteristic function of the relative elevation $V(\mathbf{y}') - V(0)$. This effective scattering operator has already been formally derived in [183] to describe the mean intensity of the diffusive reflected wavefield under the Born approximation. From the Dirac masses in (6.24), the resulting speckle can be observed, for a given direction $\mathbf{p}$, at position $\mathbf{y}_{\mathbf{p}}^{\text{ref}}$ and time $s_{\mathbf{p}}^{\text{ref}} := \mathbf{p} \cdot \mathbf{y}_{\mathbf{p}}^{\text{ref}} / 2$. Both of these quantities depend on the distance $z_{\text{int}}$ from the source term to the interface location, and therefore influence the cone size of the reflected speckle and the time window to observe its evolution.
<figure><img src="image_39.png" /><figcaption>Fig. 6.10: Illustration of the incident angle $\theta_{inc}^0$ (dash-dot line) and the reflected angle $\theta_{ref}(\mathbf{p}) \in (0, \pi/2)$ for $\gamma = 1$ and $\mathbf{p} = \beta_1\mathbf{k}_0 + \beta_2\mathbf{k}_0^\perp$, with $(\beta_1, \beta_2) = r(\cos(\varphi), \sin(\varphi))$.</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>Fig. 6.11: Illustration of the generalized Snell’s laws of refraction and transmission.</figcaption></figure>

**Generalized Snell’s law of refraction** From the observation point $\mathbf{y}_{\mathbf{p}}^{ref}$, the standard reflection relation between the incident and reflected angle (6.11) can be generalized as follows. For a given frequency $\omega$ and nonnull slowness vector $\mathbf{k}_0$, ($\theta_{inc} > 0$), we have

$$ \frac{\tan(\theta_{ref}(\mathbf{p}))}{\tan(\theta_{inc})} = \frac{|\mathbf{x}_{int} + \varepsilon^{1-\gamma}\mathbf{y}_{\mathbf{p}}^{ref}|}{|\mathbf{x}_{int}|} = \sqrt{\left(1 + \varepsilon_{\varepsilon} \frac{\mathbf{p} \cdot \mathbf{k}_0}{\cos^2(\theta_{inc})}\right)^2 + \varepsilon_{\varepsilon}^2 (\mathbf{p} \cdot \mathbf{k}_0^{\perp})^2}, \quad \varepsilon_{\varepsilon} := \varepsilon^{1-\gamma} \frac{c_0^2}{\sin^2(\theta_{inc})} $$

for $\mathbf{p}$ distributed according to the scattering operator $\mathcal{A}(2\lambda_0, \omega, \cdot)$. The case $\gamma \in (1/2, 1)$ yields the perturbation formula

$$ \theta_{ref}(\mathbf{p}) = \theta_{inc} + \varepsilon_{\varepsilon} \tan(\theta_{inc}) \mathbf{p} \cdot \mathbf{k}_0 + \mathcal{O}(\varepsilon^{2(1-\gamma)}|\mathbf{p}|^2). $$

We refer to Fig. 6.10 for an illustration of $\theta_{ref}(\mathbf{p})$ when $\gamma = 1$ (that is $\lambda \sim l_c$). From these illustrations, one can remark that the cone produced by the reflected speckle can be fairly large, depending on the norms of the scattered directions $\mathbf{p}$. For null slowness vector $\mathbf{k}_0$ ($\theta_{inc} = 0$), the refraction formula is simply given by $\tan(\theta_{ref}(\mathbf{p})) = \varepsilon^{1-\gamma}c_0|\mathbf{p}|$, that is

$$ \theta_{ref}(\mathbf{p}) = \arctan(\varepsilon^{1-\gamma}c_0|\mathbf{p}|). $$

In a similar manner, the transmitted speckle profile is given by

$$ S^{\varepsilon, tr}(\bar{s}, \bar{\mathbf{y}}, \mathbf{y}, \tilde{s}, \tilde{\mathbf{y}}) := \varepsilon^{-2(\gamma-1/2)} u^{\varepsilon, tr}(t_{obs,tr}^{\varepsilon}(\bar{s}, \bar{\mathbf{y}}, \mathbf{y}, \tilde{s}) + \varepsilon \tilde{s}, \mathbf{x}_{obs,tr}^{\varepsilon}(\bar{\mathbf{y}}, \mathbf{y}) + \varepsilon^{\gamma} \tilde{\mathbf{y}}, z = z_{tr}), $$

where

$$ \mathbf{x}_{obs,tr}^{\varepsilon}(\mathbf{y}, \bar{\mathbf{y}}) := \mathbf{x}_{obs,tr} + \varepsilon^{1-\gamma} \bar{\mathbf{y}} + \sqrt{\varepsilon} \mathbf{y}, $$
<figure><img src="image_41.png" /><figcaption>Fig. 6.12: Illustration to the zone where the transmitted speckle profile can be observed (in light grey) compared to the specular reflection component (in grey).</figcaption></figure>

and

$$t_{obs,tr}^{\epsilon}(\bar{s}, \bar{y}, \bar{y}, \tilde{y}) := t_{obs,tr} + \epsilon^{1-\gamma} \mathbf{k}_0 \cdot \bar{y} + \sqrt{\epsilon} \mathbf{k}_0 \cdot \bar{y} + \epsilon^{\gamma} \mathbf{k}_0 \cdot \tilde{y} + \epsilon^{2(1-\gamma)} \bar{s}.$$

The corresponding correlation function is defined as:

$$C^{\epsilon,tr}(\bar{s}, \bar{y}, \bar{y}, \tilde{s}_1, \tilde{y}_1, \tilde{s}_2, \tilde{y}_2) := S^{\epsilon,tr}(\bar{s}, \bar{y}, \bar{y}, \tilde{s}_1, \tilde{y}_1)S^{\epsilon,tr}(\bar{s}, \bar{y}, \bar{y}, \tilde{s}_2, \tilde{y}_2),$$

where its asymptotic behavior is given by the following result.

**Proposition 6.2.6 We have**

$$\lim_{\epsilon \to 0} C^{\epsilon,tr}(\bar{s}, \bar{y}, \bar{y}, \tilde{s}_1, \tilde{y}_1, \tilde{s}_2, \tilde{y}_2) = C^{tr}(\bar{s}, \bar{y}, \tilde{s}_1 - \tilde{s}_2, \tilde{y}_1 - \tilde{y}_2),$$

in probability in $S'(\mathbb{R} \times \mathbb{R}^2 \times \mathbb{R}^2 \times \mathbb{R} \times \mathbb{R}^2 \times \mathbb{R} \times \mathbb{R}^2)$, where

$$C^{tr}(\bar{s}, \bar{y}, \bar{y}, \tilde{s}, \tilde{y}) := \frac{\mathcal{T}^2 \lambda_0}{4\lambda_1(2\pi)^3} \iint e^{-i\omega(\tilde{s}-\mathbf{p}\cdot\tilde{\mathbf{y}})} \mathcal{A}(\lambda_0 - \lambda_1, \omega, \mathbf{p}) |\hat{\Psi}|_2^2(\omega) \delta(\bar{s} - s_{\mathbf{p}}^{tr}) \delta(\bar{y} - \mathbf{y}_{\mathbf{p}}^{tr}) \omega^4 d\omega d\mathbf{p}.$$

with $\mathcal{A}$ defined by (6.25), $|\hat{\Psi}|_2^2$ by (6.26), and

$$\mathbf{y}_{\mathbf{p}}^{tr} := 2(z_{tr} - z_{int})c_1 A_1 \mathbf{p} \geq 0, \quad \text{and} \quad s_{\mathbf{p}}^{tr} := \mathbf{p} \cdot \mathbf{y}_{\mathbf{p}}^{tr} / 2 = (z_{tr} - z_{int})c_1 \mathbf{p}^T A_1 \mathbf{p}.$$

Accordingly, the mean intensity carried by the speckle profile is given by

$$\lim_{\epsilon \to 0} |S^{\epsilon,tr}(\bar{s}, \bar{y}, \bar{y}, \tilde{s}, \tilde{y})|^2 = C^{tr}(\bar{s}, \bar{y}, 0, 0),$$

in probability, which does not depend on the variable $\mathbf{y}$ and the small scale fluctuations in $\tilde{s}$ and $\tilde{y}$. The intensity is therefore again uniform over these scales.

**Generalized Snell's law of transmission** From $\mathbf{y}_{\mathbf{p}}^{tr}$, we obtain the following relation between the transmission angle and the specular transmission angle $\theta_{tr}^0$:

$$\frac{\tan(\theta_{tr}(\mathbf{p}))}{\tan(\theta_{tr}^0)} = \frac{|\mathbf{x}_{obs,tr} - \mathbf{x}_{int} + \epsilon^{1-\gamma}\mathbf{y}_{\mathbf{p}}^{tr}|}{|\mathbf{x}_{obs,tr} - \mathbf{x}_{int}|} = \sqrt{\left(1 + \xi_{\epsilon} \frac{\mathbf{p} \cdot \mathbf{k}_0}{\cos^2(\theta_{tr})}\right)^2 + \xi_{\epsilon}^2 (\mathbf{p} \cdot \mathbf{k}_0^{\perp})^2}, \quad \xi_{\epsilon} := \epsilon^{1-\gamma} \frac{c_0^2}{\sin^2(\theta_{inc})},$$

for $\mathbf{p}$ distributed according to the scattering operator $\mathcal{A}(\lambda_0 - \lambda_1, \omega, \cdot)$. From this relation, the standard Snell's relation of transmission can be generalized as follow:

$$\frac{\sin(\theta_{tr}(\mathbf{p}))}{c_1} = \frac{\sin(\theta_{inc})}{c_0} \sqrt{\frac{\Xi}{1 + \sin^2(\theta_{tr}^0)(\Xi - 1)}}, \quad \Xi := \left(1 + \xi_{\epsilon} \frac{\mathbf{p} \cdot \mathbf{k}_0}{\cos^2(\theta_{tr})}\right)^2 + \xi_{\epsilon}^2 (\mathbf{p} \cdot \mathbf{k}_0^{\perp})^2,$$
<figure><img src="image_42.png" /><figcaption>Fig. 6.13: Illustration of the specular transmission angle $\theta_{tr}^0$ (dash-dot line) and the transmission angle $\theta_{ref}(\mathbf{p}) \in (0, \pi/2)$ for $\gamma = 1$, $c_0 = 1.5$, $c_1 = 1$, and $\mathbf{p} = \beta_1 \mathbf{k}_0 + \beta_2 \mathbf{k}_0^\perp$ with $(\beta_1, \beta_2) = r(\cos(\varphi), \sin(\varphi))$.</figcaption></figure>

leading for $\gamma \in (1/2, 1)$ to the perturbed relation

$$ \frac{\sin(\theta_{tr}(\mathbf{p}))}{c_1} = \frac{\sin(\theta_{inc})}{c_0} \left( 1 + \xi_{\varepsilon} \frac{\mathbf{p} \cdot \mathbf{k}_0}{\cos^2(\theta_{tr}^0)} + \mathcal{O}(\varepsilon^{2(1-\gamma)} |\mathbf{p}|^2) \right), $$

providing the approximation

$$ \theta_{tr}(\mathbf{p}) = \theta_{tr}^0 + \xi_{\varepsilon} \frac{\tan(\theta_{tr})}{\cos^2(\theta_{tr})} \mathbf{p} \cdot \mathbf{k}_0 + \mathcal{O}(\varepsilon^{2(1-\gamma)} |\mathbf{p}|^2). $$

We refer to Fig. 6.13 for an illustration of $\theta_{tr}(\mathbf{p})$ when $\gamma = 1$, that is in the context of a rough interface $\lambda \sim l_c$. From these illustrations, one can again remark that the cone produced by the transmitted speckle can also be fairly large, depending on the norms of the scattered directions $\mathbf{p}$. Finally, for a null slowness vector $\mathbf{k}_0 = 0$ ($\theta_{inc} = 0$), the transmission angle satisfies $\tan(\theta_{tr}(\mathbf{p})) = \varepsilon^{1-\gamma} c_1 |\mathbf{p}|$, so that

$$ \theta_{tr}(\mathbf{p}) = \arctan(\varepsilon^{1-\gamma} c_1 |\mathbf{p}|). $$

### 6.2.6 Gaussian statistics of the speckle patterns for $l_c \ll r_0$ ($\gamma > 1/2$)

We conclude this section with a direct statistical description of the speckles themselves in the Fourier domain at the scale of the pulse width and correlation length. For the reflected speckle, we consider

$$ \hat{S}^{\varepsilon, \text{ref}}(\bar{s}, \bar{y}, \mathbf{y}, \omega, \mathbf{p}) := \varepsilon^{-2(\gamma-1/2)} \iint e^{i\omega(\tilde{s}-\mathbf{p}\cdot\tilde{\mathbf{y}})} u^{\varepsilon, \text{ref}}(t_{\text{obs,ref}}^{\varepsilon}(\bar{s}, \bar{y}, \mathbf{y}, \tilde{\mathbf{y}}) + \varepsilon\tilde{s}, \mathbf{x}_{\text{obs,ref}}^{\varepsilon}(\bar{y}, \mathbf{y}) + \varepsilon^{\gamma}\tilde{\mathbf{y}}) d\tilde{s} d\tilde{\mathbf{y}}, $$

where $\mathbf{x}_{\text{obs,ref}}^{\varepsilon}$ and $t_{\text{obs,ref}}^{\varepsilon}$ are respectively given by (6.21) and (6.22). To regularize the Dirac masses in (6.24) around $\mathbf{y}_{\mathbf{p}}^{\text{ref}}$ and $\mathbf{s}_{\mathbf{p}}^{\text{ref}}$, we introduce

$$ \hat{S}_{\mathbf{y}}^{\varepsilon, \text{ref}}(\bar{s}, \bar{y}, \omega, \mathbf{p}) := \hat{S}^{\varepsilon, \text{ref}}(\bar{s}, \bar{y}, \mathbf{y}, \omega, \mathbf{p}) \frac{1}{\varepsilon^{3(\gamma-1/2)}} \varphi^{1/2} \left( 2 \frac{\bar{s} - s_{\mathbf{p}}^{\text{ref}}}{\varepsilon^{2(\gamma-1/2)}}, 2 \frac{\bar{y} - y_{\mathbf{p}}^{\text{ref}}}{\varepsilon^{2(\gamma-1/2)}} \right), \quad (6.27) $$

where we have added the square root of a symmetric mollifier, that acts as a windowing around $\bar{s} = s_{\mathbf{p}}^{\text{ref}}$ and $\bar{y} = y_{\mathbf{p}}^{\text{ref}}$. Recalling that $\bar{s}$ and $\bar{y}$ scale at respectively $\varepsilon^{2(1-\gamma)}$ and $\varepsilon^{1-\gamma}$, the resulting window scales at the pulse width $\varepsilon$ and correlation length $\varepsilon^{\gamma}$ respectively. The asymptotic behavior of the speckle is the following.

**Theorem 6.2.1** For $n \ge 1$ and any fixed $\mathbf{y}_1, \dots, \mathbf{y}_n, \in \mathbb{R}^2$, the family $(\hat{S}_{\mathbf{y}_1}^{\varepsilon, \text{ref}}, \dots, \hat{S}_{\mathbf{y}_n}^{\varepsilon, \text{ref}})_\varepsilon$ converges in law, as $\varepsilon \to 0$, in $S'_n$ (which is $n$ times the Cartesian product of $S'(\mathbb{R} \times \mathbb{R}^2 \times \mathbb{R} \times \mathbb{R}^2, \mathcal{C})$) to a limit $(\hat{S}_1^{\text{ref}}, \dots, \hat{S}_n^{\text{ref}})$ made of $n$ independent copies of a complex valued mean-zero Gaussian random field
$\hat{S}^{ref}$ (independent of the $\mathbf{y}_j$'s) with covariance functions of the form

$$ \mathbb{E}[\hat{S}^{ref}(\phi)\hat{S}^{ref}(\psi)] = \int \dots \int \hat{K}_{ref}(\bar{s}_1, \bar{s}_2, \bar{\mathbf{y}}_1, \bar{\mathbf{y}}_2, \omega_1, -\omega_2, \mathbf{p}_1, \mathbf{p}_2) \\ \times \overline{\phi(\bar{s}_1, \bar{\mathbf{y}}_1, \omega_1, \mathbf{p}_1)\psi(\bar{s}_2, \bar{\mathbf{y}}_2, \omega_2, \mathbf{p}_2)} d\bar{s}_1 d\bar{s}_2 d\bar{\mathbf{y}}_1 d\bar{\mathbf{y}}_2 d\omega_1 d\omega_2 d\mathbf{p}_1 d\mathbf{p}_2, \quad (6.28) $$

$$ \mathbb{E}[\hat{S}^{ref}(\phi) \overline{\hat{S}^{ref}(\psi)}] = \int \dots \int \hat{K}_{ref}(\bar{s}_1, \bar{s}_2, \bar{\mathbf{y}}_1, \bar{\mathbf{y}}_2, \omega_1, \omega_2, \mathbf{p}_1, \mathbf{p}_2) \\ \times \overline{\phi(\bar{s}_1, \bar{\mathbf{y}}_1, \omega_1, \mathbf{p}_1)\psi(\bar{s}_2, \bar{\mathbf{y}}_2, \omega_2, \mathbf{p}_2)} d\bar{s}_1 d\bar{s}_2 d\bar{\mathbf{y}}_1 d\bar{\mathbf{y}}_2 d\omega_1 d\omega_2 d\mathbf{p}_1 d\mathbf{p}_2, \quad (6.29) $$

for any $\phi, \psi \in S(\mathbb{R} \times \mathbb{R}^2 \times \mathbb{R} \times \mathbb{R}^2, \mathcal{C})$, where the kernel $\hat{K}_{ref}$ is given by

$$ \hat{K}_{ref}(\bar{s}_1, \bar{s}_2, \bar{\mathbf{y}}_1, \bar{\mathbf{y}}_2, \omega_1, \omega_2, \mathbf{p}_1, \mathbf{p}_2) = \frac{(2\pi)^3 \mathcal{R}^2}{4} A(2\lambda_0, \omega_1, \mathbf{p}_1) |\hat{\Psi}|_2^2(\omega_1) \hat{\phi}(\omega_1, \mathbf{p}_1) \delta(\omega_1 - \omega_2) \times \delta(\mathbf{p}_1 - \mathbf{p}_2) \delta(\bar{s}_1 - s_{\mathbf{p}_1}^{ref}) \delta(\bar{s}_2 - \bar{s}_1) \delta(\bar{\mathbf{y}}_1 - \mathbf{y}_{\mathbf{p}_1}^{ref}) \delta(\bar{\mathbf{y}}_2 - \bar{\mathbf{y}}_1). \quad (6.30) $$

In contrast to the correlation function (6.23), exhibiting a self-averaging property, the asymptotic behavior of the speckle (6.27) is described through a Gaussian random field w.r.t. the scattered direction $\mathbf{p}$. The term $\delta(\mathbf{p}_1 - \mathbf{p}_2)$ in (6.30) translates that two distinct directions produce independent components. This is of course expected as these two components will be observed at locations far apart from each other. The $\mathbf{p}$-variable being also related to fluctuations at the scale of the correlation length, the resulting Gaussian field is then described at this scale around the given time $s_{\mathbf{p}}^{ref}$ and location $\mathbf{y}_{\mathbf{p}}^{ref}$.

A similar result can be obtained for the transmitted speckle in the Fourier domain,

$$ \hat{S}^{\varepsilon, \text{tr}}(\bar{s}, \bar{\mathbf{y}}, \mathbf{y}, \omega, \mathbf{p}) := \varepsilon^{-2(\gamma-1/2)} \iint e^{i\omega(\tilde{\mathbf{s}}-\mathbf{p}\cdot\tilde{\mathbf{y}})} u^{\varepsilon, \text{tr}}(t_{\text{obs, tr}}^{\varepsilon}(\bar{s}, \bar{\mathbf{y}}, \mathbf{y}, \tilde{\mathbf{y}}) + \varepsilon\tilde{s}, \mathbf{x}_{\text{obs, tr}}^{\varepsilon}(\bar{\mathbf{y}}, \mathbf{y}) + \varepsilon^{\gamma}\tilde{\mathbf{y}}) d\tilde{s} d\tilde{\mathbf{y}}, $$

after introducing a window around $\mathbf{y}_{\mathbf{p}}^{tr}$ and $s_{\mathbf{p}}^{tr}$,

$$ \hat{S}_{\tilde{\mathbf{y}}}^{\varepsilon, \text{tr}}(\bar{s}, \bar{\mathbf{y}}, \omega, \mathbf{p}) := \hat{S}^{\varepsilon, \text{tr}}(\bar{s}, \bar{\mathbf{y}}, \mathbf{y}, \omega, \mathbf{p}) \frac{1}{\varepsilon^{3(\gamma-1/2)}} \varphi^{1/2} \left( 2 \frac{\bar{s} - s_{\mathbf{p}}^{tr}}{\varepsilon^{2(\gamma-1/2)}}, 2 \frac{\bar{\mathbf{y}} - \mathbf{y}_{\mathbf{p}}^{tr}}{\varepsilon^{2(\gamma-1/2)}} \right). $$

**Theorem 6.2.2** For $n \ge 1$ and any fixed $\mathbf{y}_1, \dots, \mathbf{y}_n \in \mathbb{R}^2$, the family $(\hat{S}_{\tilde{\mathbf{y}}_1}^{\varepsilon, \text{tr}}, \dots, \hat{S}_{\tilde{\mathbf{y}}_n}^{\varepsilon, \text{tr}})_\varepsilon$ converges in law, as $\varepsilon \to 0$, in $S'_n$ to $(\hat{S}_1^{tr}, \dots, \hat{S}_n^{tr})_\varepsilon$, made of $n$ independent copies of a complex mean-zero Gaussian random field $\hat{S}^{ref}$ with covariance functions similar to (6.28) and (6.29), but with kernel

$$ \hat{K}_{\text{tr}}(\bar{s}_1, \bar{s}_2, \bar{\mathbf{y}}_1, \bar{\mathbf{y}}_2, \omega_1, \omega_2, \mathbf{p}_1, \mathbf{p}_2) = \frac{(2\pi)^3 \mathcal{T}^2 \lambda_0}{4\lambda_1} A(\lambda_0 - \lambda_1, \omega_1, \mathbf{p}_1) |\hat{\Psi}|_2^2(\omega_1) \hat{\phi}(\omega_1, \mathbf{p}_1) \delta(\omega_1 - \omega_2) \\ \times \delta(\mathbf{p}_1 - \mathbf{p}_2) \delta(\bar{s}_1 - s_{\mathbf{p}_1}) \delta(\bar{s}_2 - s_{\mathbf{p}_1}) \delta(\bar{\mathbf{y}}_1 - \mathbf{y}_{\mathbf{p}_1}) \delta(\bar{\mathbf{y}}_2 - \mathbf{y}_{\mathbf{p}_1}). $$

## 6.3 Perspectives

In the second part of this chapter, the statistical properties of the reflected and transmitted wavefields are described in the case of a moderately rough interface (the central wavelength and the magnitude of the interface fluctuations are of the same order $\lambda \sim \sigma$). In this regime, effective specular and diffusive (speckle) components are observed, but no enhance backscattering effects. To observe such effects, it would be necessary for the rough surface to produce significant multiple scattering, which is achieved when $\lambda \ll \sigma$.

This novel framework precisely characterizes the statistics of the diffusive scattered waves from rough interfaces, both in transmission and reflection, in term of a Gaussian random field. This characterization can be used in remote sensing type imaging problems including: estimating the parameters of the rough (random) interface, and imaging objects hidden behind the interface by computing the empirical spectrum of the reflected speckle [88, 127]. The analysis can also be extended to Synthetic Aperture Radar (SAR) probing scenarios where the source cone is moving, typical mounted on an aircraft. In general imaging contexts one can exploit the memory effect, which means that the speckle pattern illuminating a hidden object is not completely changed, but rather shifted by a specific amount when the incoming angle of the source is shifted.
Chapter 7

Digressions in mathematical biology:
Stochastic and deterministic models for the
secondary metastatic emission process

The emergence of a metastatic phenotype is governed by a number of key mutations that was trans-
lated into mathematical models by [107,160], resulting in the derivation of a metastatic risk score. In
an early work, [136] established a link between the primary tumor size at surgery and the risk of re-
currence from a large cohort of breast cancer patients. In addition to such phenotypic characteristics,
specific genetic signatures of the primary tumor have been identified as associated with increased risk
of metastasis [204]. Nevertheless, these risk prediction models are static; they do not aim to represent
the time evolution of the disease. In contrast, dynamic models offer the ability to predict how a
system evolves over time, making it easier to incorporate data obtained at various observation points.
To account for the relationship between the metastatic process and the primary tumor, a dynamic de-
scription of primary tumor growth is integrated into metastatic models. The simplest dynamic model
for tumor growth is the exponential model, which describes growth under no restriction such as in
vitro settings. However, in many cases of interest, especially in vivo scenarios, sigmoid (S-shaped)
models with an initial exponential phase followed by a deceleration provide more suitable descriptions
of growth dynamics. The Gompertz model is a classical example that has been commonly used, as
well as power growth models that have been used for the description of clinical [111] and preclinical
tumor growth data [32].

A stochastic dynamic model for metastasis was proposed by [28]. Their approach characterizes the emission times of metastases as random events through a non-homogeneous Poisson process with an emission rate increasing according to the primary tumor size. The Poisson law allows for an interpretation of the emission process as being memoryless, and provide fairly explicit mathematical formulations. A variant of this model successfully described data on bone lesions in the case of a metastatic breast cancer [108]. Dynamic models for metastasis can also be used to infer parts of the process that cannot be directly observed experimentally or clinically. For instance, [122] introduced a PDE model to describe the size distribution of metastatic colonies. Using clinical information concerning visible metastatic colonies, they successfully characterized the micrometastatic state of a patient with hepatocellular carcinoma. Later, this size-structured model was integrated into mixed-effects models, providing prediction of the size evolution of metastases in animal models without [113] and with [33] surgical removal of the primary tumor. Both approaches incorporate secondary metastatic emission into the model, that is the ability for metastasis to spread further. Indirect evidence of secondary metastatic emission has been provided through cancer network models [50,164, 165] and the identification of self-seeding mechanisms [55,186].

In most of the models mentioned above, metastases do not have the ability to emit metastases
themselves. However, it is easy to think of situations where this property would make a difference in
the model. For instance, consider a scenario in which a single metastasis is emitted before the primary
tumor surgery<fnref n="1" />. If this metastasis cannot generate further metastases, its successful removal cures the patient. However, the second surgery might fail if the metastasis is able to seed additional metastasis as well. Of course, there are other mechanisms potentially leading to treatment failure (e.g. local recurrence, surgery impossible,...), but for simplicity they are not discussed here. This section presents the relationship between the Iwata PDE model, which accounts for secondary metastatic emission, and a stochastic model based on a cascade of Poisson processes. Under this stochastic framework, the cascade effect accounts for secondary metastatic emissions, and we illustrate how this stochastic approach can be used to assess metastatic risk post-surgery.

## 7.1 A simple probabilistic framework for metastatic risk

Forecasting the likelihood of metastatic disease at the time of the primary tumor's diagnosis is of major clinical importance since it is strongly linked to survival expectancy. One possibility to construct such a predictive model is by using large databases to correlate information regarding the presence of metastases and the primary tumor's characteristics at the time of diagnosis or surgery. As an illustration, [159] proposed a relationship between the primary tumor size at surgery and the probability of metastasis based on clinical data on breast cancer:

$$ \mathbb{P}(\text{no metastases}) = \exp(-c d^z), $$

where d is the largest diameter of the primary tumor at surgery, and c and z are parameters inferred from the data. Such a probability formula can be derived through the use of a nonhomogeneous Poisson process with an intensity function depending on a power growth model for the primary tumor.

**Power growth.** The growth diameter d of the largest primary tumor is the solution of the ordinary differential equation $d'(t) = a_{PG} d(t)^{\alpha}$. Here, $a_{PG}$ corresponds to the growth speed and $\alpha$ allows to describe different growth shapes: exponential growth ($\alpha = 1$), linear growth ($\alpha = 0$), and a range of sigmoid growth patterns in between ($0 < \alpha < 1$).

**Power law of emission.** The emission intensity $\lambda$ of the Poisson process depends on the current size of the primary tumor through a power law, $\lambda(t) = b d(t)^{\beta}$. In this context, b can be understand as the metastatic aggressiveness of the primary tumor. In [122], the parameter $\beta$ is linked to the type of primary tumor vascularization: $\beta = 3$ (dimension of space) for a uniform vascularization, and $\beta = 2$ (dimension of a surface)$^2$ for a surficial vascularization.

The stochastic nature of the metastatic emission process can also be described using various probability distributions, allowing to drop the memorylessness property. Nevertheless, the Poisson model offers the advantage of not necessitating extra statistical parameters, while ensuring a high level of analytical tractability and ease of numerical implementation.

## 7.2 Mathematical formalism and results

### 7.2.1 Size-structured model

To describe the micrometastatic state of cancer patients, [122] developed a size-structured model that describes the time evolution of a density function $\rho(t, x)$ representing the size distribution of metastatic colonies. Before describing this model, we need to introduce some notations.

The relationship between the size of the primary tumor $X_p$ and its emission intensity $\lambda_p$ is here given by a size-dependent emission law $\gamma_p$,

$$ \lambda_p(t) := \gamma_p(X_p(t)). $$

<footnote n="1">¹After a surgery the emission intensity is set to zero.</footnote>

<footnote n="2">²The interpretation of $\beta$ depends on the unit of the primary tumor measure. A surficial vascularization would correspond to $\beta = 2$ if the size is measured in diameter, but to $\beta = 2/3$ (the fractal dimension of a surface in space) if the size is measured in volume.</footnote>
The evolution of $X_p$ is given by a Gompertz law:

$$X'_{p} = g(X_{p}), \quad \text{with} \quad X_{p}(0) = x_{0},$$

with a set of clinical parameters summarized in Table 7.1. These parameters were estimated by [122] from clinical data on a hepatocellular carcinoma with multiple liver metastases. In a similar manner,

<table><caption>Table 7.1: Growth and emission laws derived by [122] from clinical data of a hepatocellular carcinoma with multiple liver metastases.</caption><thead><tr><th>Model</th><th>Parameter</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td rowspan="3">Growth (Gompertz law)<br>g(x) := a<sub>Gomp</sub> x log(x<sub>p</sub><sup>∞</sup>/x)</td><td>Initial size</td><td>x<sub>0</sub></td><td>1</td><td>cells</td></tr><tr><td>Growth rate</td><td>a<sub>Gomp</sub></td><td>0.00286</td><td>days<sup>-1</sup></td></tr><tr><td>Maximum size</td><td>x<sub>p</sub><sup>∞</sup></td><td>7.3 &middot; 10<sup>10</sup></td><td>cells</td></tr><tr><td rowspan="2">Emission (power law)<br>γ<sub>p</sub>(x) = b x<sup>β</sup></td><td>Rate constant</td><td>b</td><td>5.3 &middot; 10<sup>-8</sup></td><td>days<sup>-1</sup> cells<sup>-β</sup></td></tr><tr><td>Emission power</td><td>β</td><td>0.663</td><td>-</td></tr></tbody></table>

the size of a metastasis is given by $X_m$, which is the solution to

$$X'_{m}(t) = g(X_{m}(t)), \quad \text{with} \quad X_{m}(0) = x_{m}^{0},$$

and corresponding secondary emission intensity by

$$\lambda_{m}(t) := \gamma_{m}(X_{m}(t)).$$

The metastatic density function $\rho$ of the Iwata model solves the following PDE:

$$\left\{
\begin{array}{ll}
\partial_t \rho(x,t) + \partial_x[g(x)\rho(x,t)] = 0, & (x,t) \in (x_m^0, x_m^\infty) \times (0, +\infty), \\
g(x_m^0)\rho(x_m^0, t) = \gamma_p(X_p(t)) + \int_{x_m^0}^{x_m^\infty} \gamma_m(x)\rho(x,t)dx, & t \in (0, +\infty), \\
\rho(x,0) = 0, & x \in [x_m^0, x_m^\infty],
\end{array}
\right.
\quad (7.1)$$

where $x_m^0$ corresponds to the size of a newly emitted metastasis and $x_m^\infty$ corresponds to the theoretical maximum size of the metastases. The well-posedness for this model has been established under general conditions by [27].

Various *model observable* (MO) associated to this model can be defined: for a given function $f$, a MO is defined by

$$\mathrm{MO}_f(t) := \int_{x_m^0}^{x_m^\infty} f(x)\rho(x,t)dx.$$

Here are some examples.

* $f = 1$ corresponds to the number of metastases:

$$\mathrm{MO}_1(t) := \int_{x_m^0}^{x_m^\infty} 1 \cdot \rho(x,t)dx$$

* Similarly, the number of macrometastases is obtained with

$$f_{\text{macro}}(x) := \begin{cases} 1 & \text{if } x \ge c \\ 0 & \text{if } x < c \end{cases},$$

and the number of micrometastases with

$$f_{\text{micro}}(x) := \begin{cases} 0 & \text{if } x \ge c \\ 1 & \text{if } x < c \end{cases},$$

where $c$ stands for the detectability threshold (depending on the imaging modality).
• The total metastatic mass M is obtained with the identity function $f_{\text{Id}}(x) = x$ for all $x$:

$$ \mathrm{MO}_{f_{\text{Id}}}(t) := \int_{x_m^0}^{x_m^\infty} x \rho(x, t) dx. $$

The PDE model (7.1) is not convenient when it comes to evaluate MO numerically. Fortunately, all the MO can be characterized as being the solutions of a Volterra convolution equation [112].

**Theorem 7.2.1** For any $f \in L^\infty([x_m^0, x_m^\infty])$, $\mathrm{MO}_f$ is the unique solution of the following renewal equation,

$$ \mathrm{MO}_f(t) = \int_0^t \lambda_p(s) f(X_m(t-s)) ds + \int_0^t \lambda_m(s) \mathrm{MO}_f(t-s) ds. \quad (7.2) $$

This representation is the basis of an efficient numerical resolution algorithm [112] that can be used for parameter estimation [113].

## 7.2.2 Probabilistic framework for secondary metastatic emission

We now introduce the stochastic counter part of the Iwata model. This stochastic model is based on a cascade of independent Poisson processes (PPs), and define recursively the emission times with respect to a generational hierarchy.

• The emission times $(T^{(j)})_{j \ge 1}$ at which the first generation of metastases occur, those emanating from the primary tumor, are modeled as the jump times of a PP $(N^{(1)})_{t \ge 0}$ with intensity $\lambda_p$.

The emission times for the next generations of metastases are defined recursively.

• Let $k \ge 2$ and $n_1, \dots, n_{k-1} \ge 1$. The j-th emission time for the k-th generation of metastasis and filiation $n_1, \dots, n_{k-1}$ is defined by

$$ T^{(n_1, \dots, n_{k-1}, j)} := T^{(n_1, \dots, n_{k-1})} + \tilde{T}^{(n_1, \dots, n_{k-1}, j)}. $$

This corresponds to the time it takes for the offspring, with filiation $n_1, \dots, n_{k-1}$, to give birth to its j-th offspring. Here, $(\tilde{T}^{(n_1, \dots, n_{k-1}, j)})_{j \ge 1}$ are the jump times of a PP $(N_t^{(n_1, \dots, n_{k-1})})_{t \ge 0}$ with intensity $\lambda_m$.

For instance, $T^{(2,3,4)}$ corresponds to the emission time of the 4th offspring produced by the 3rd offspring of the 2nd offspring of the primary tumor. Also, we assume that

$$ \{(N_t^{(n_1, \dots, n_k)})_{t \ge 0}, \ k \ge 1, \ n_1, \dots, n_k \ge 1\} $$

forms a family of independent PPs. In other words, we assume that the primary tumors and all the metastases emit independently from each other.

From this construction, a stochastic model observable (SMO) for the k-th generation can be expressed by

$$ \mathrm{SMO}_f^{(k)}(t) := \sum_{n_1, \dots, n_k \ge 1} \mathbf{1}_{(T^{(n_1, \dots, n_k)} \le t)} f(X_m(t - T^{(n_1, \dots, n_k)})), $$

for $f \in L^\infty([x_m^0, x_m^\infty])$. Therefore, a global SMO with secondary emission can be defined by

$$ \mathrm{SMO}_f(t) := \sum_{k=1}^{\infty} \mathrm{SMO}_f^{(k)}(t). \quad (7.3) $$

A global SMO accounts for the contribution of $\mathrm{SMO}_f^{(k)}$ over all the generations. The following result establish the connection between the SMO and MO.

**Proposition 7.2.1** Let $f \in L^\infty([x_m^0, x_m^\infty])$. The SMO (7.3) is well defined in the sense that

$$ \mathbb{P}(0 \le \mathrm{SMO}_f(t) < +\infty \quad \forall t \ge 0) = 1, $$

and $\mathbb{E}[\mathrm{SMO}_f(t)]$ is finite for any time $t \ge 0$. Moreover, we have for any $t \ge 0$

$$ \mathrm{MO}_f(t) = \mathbb{E}[\mathrm{SMO}_f(t)]. $$
Let us remark that all the SMO can be seen as integrals w.r.t. a random measure:

$$ \text{SMO}_f(t) = \int f(x) \mathcal{M}_t(dx), $$

for any $t \ge 0$, with

$$ \mathcal{M}_t := \sum_{k \ge 1} \sum_{n_1, \dots, n_k \ge 1} \delta_{X_m(t-T^{(n_1, \dots, n_k)})}. $$

From this representation, we can establish the connection between the stochastic model and structure size model (Iwata model).

**Theorem 7.2.2 (Link to the structured population model) For all $t \ge 0$, the measure**

$$ \mu_t := \mathbb{E}[\mathcal{M}_t] $$

is $\sigma$-finite, absolutely continuous with respect to the Lebesgue measure, and its Radon-Nikodyn density is given by $\rho(\cdot, t)$,

$$ \frac{d\mu_t}{dx} = \rho(\cdot, t), $$

where $\rho$ is the solution to (7.1).

From Prop. 7.2.1, one can remark that all the SMO can be seen as a probabilistic representation of the corresponding MO. This remark is of interest to design rough confidence/prediction intervals from inferred parameters. In this respect, the following result concerns the variability of the SMO with respect to its mean MO.

**Proposition 7.2.2** Let $f \in L^\infty([x_m^0, x_m^\infty])$. The variance of the SMO,

$$ v_f(t) := \text{var}[\text{SMO}_f(t)] $$

is finite for any $t \ge 0$, and satisfies a renewal equation,

$$ v_f(t) = \int_0^t \lambda_p(s)(f(X_m(t-s)) + e_{m,f}(t-s))^2 + \int_0^t \lambda_m(s)v_f(t-s), \quad (7.4) $$

Here,

$$ e_{m,f}(t) := \mathbb{E}[\text{SMO}_{m,f}(t)], $$

where $\text{SMO}_{m,f}$ is defined as (7.3), but for a different cascade of PPs, which has only $\lambda_m$ for intensity (for both the first and subsequent generations).

The renewal equation (7.4) allows the use of an efficient numerical resolution algorithm [112] to compute the variance of $\text{SMO}_f$ and not rely on a Monte-Carlo approach.

## 7.3 Numerical illustrations

Let us illustrate this approach with the parameters of Table 7.1, and the observables that were used in [122] to parametrize their model, that is the number of metastases exceeding given thresholds $c$ (i.e. $f_{\text{macro}}$ with different thresholds). The corresponding numerical simulations are given in Fig. 7.1. In this illustration one can observe the prediction of the size-structured model (blue line). Some paths of the PP cascade model are also displayed to illustrate the variability inherent of a stochastic model. Note that these paths evolves mainly between the two red dashed lines accounting for the variance of the model.

In this context, the stochastic framework is a valuable tool to assess the impact of secondary metastatic emission following the surgery of the primary tumor, and evaluate the risk of a metastatic disease through a Monte-Carlo approach. To represent the impact of a surgery at time $t_{\text{surgery}}$, the emission intensity $\lambda_p(t)$ of the primary tumor is set to zero for $t \ge t_{\text{surgery}}$. Using the clinical parameters of Table 7.1, we assume that the primary tumor is surgically removed 500 days after its inception, where it has reached a tumor mass of 180 g, and evaluate the number of metastases
<figure><img src="image_43.png" /><figcaption>Fig. 7.1: Comparison of the size-structured model fit and stochastic variability of the PP cascade model. Expectation (blue solid line) corresponds to the size-structured model prediction used to fit the clinical data from [122] (computed via Eq. (7.2)). Random trajectories of the PP cascade model are displayed through thin grey lines, and the red dashed lines corresponds to $E[N_{\text{macro},t}] \pm 2 \cdot \sqrt{\text{var}[N_{\text{macro},t}]}$ (with var[$N_{\text{macro},t}$] computed via Eq. (7.4)). As in [122], we count time from inception of the first primary tumor cell, which was back-calculated from primary tumor data assuming a Gompertzian growth model (hence the first CT scan with metastatic disease is approximately 3 years post-inception).</figcaption></figure>

another 500 days later (Fig. 7.2). We consider both a model with secondary metastatic emission and a model without secondary emission (where metastases cannot emit themselves new metastasis). Knowing that every secondary emission is preceded by at least one primary emission, the probability of a metastatic disease is the same for both models. However, the model with secondary emission predicts a significantly higher number of metastases in average ($E[N_t] = 4.7$ with secondary emission vs. $E[N_t] = 1.2$ without). Not surprisingly, the model with secondary metastatic emission is more pessimistic than without secondary emission.

## 7.4 Perspectives

This approach combining a deterministic and a stochastic model could benefit from further developments for individualized risk predictions, and be included as a factor (among others) in comprehensive decision frameworks based on modern machine learning methods. Other (individual) factors play a key role in the context of cancer and metastatic diseases. For instance, circulating biomarkers, such as circulating tumor cells or circulating tumor DNA, can be a useful source of information. Once validated, a mathematical model can serve as a powerful tool for informed treatment decisions by integrating case-specific information into a consistent quantitative framework.
<figure><img src="image_44.png" /><figcaption>Fig. 7.2: Probability of a metastatic disease after surgery with (top panel) or without (bottom panel) secondary metastatic emission (each based on 10.000 simulations). The primary tumor is assumed to be surgically removed 500 days after its inception. The number of metastases is evaluated another 500 days later.</figcaption></figure>
Bibliography

[1] A. Abdulle, D. Cohen, G. Vilmart, and K. C. Zygalakis, *High weak order methods for stochastic differential equations based on modified equations*, SIAM J. Sci. Comput. **34** (2012), A1800–A1823.

[2] J.-L. Akian, *Space-time semiclassical measures for three-dimensional elastodynamics: boundary conditions for the hyperbolic set*, Asymptot. Anal. **78** (2012), 37–83.

[3] R. Alexandre, *A review of Boltzmann equation with singular kernels*, Kinet. Relat. Models **2** (2009), 551–646.

[4] ——, *Fractional order kinetic equations and hypoellipticity*, Anal. Appl. **10** (2012), 237–247.

[5] R. Alexandre, Y. Morimoto, S. Ukai, C.-J. Xu, and T. Yang, *Global existence and full regularity of the Boltzmann equation without angular cutoff*, Comm. Math. Phys. **304** (2011), 513–581.

[6] R. Alexandre and C. Villani, *On the Boltzmann equation for long-range interactions*, Comm. Pure Appl. Math. **55** (2002), 30–70.

[7] R. Alonso, L. Borcea, and J. Garnier, *Wave propagation in waveguides with random boundaries*, Commun. Math. Sci. **11** (2012), 233–267.

[8] V. I. Arnold, *Mathematical methods of classical mechanics*, Springer, 1989.

[9] S. Asmussen and J. Rosiński, *Approximations of small jumps of lévy processes with a view towards simulation*, J. Appl. Probab. **38** (2001), 482–493.

[10] F. Bailly, J. F. Clouet, and J. P. Fouque, *Parabolic and gaussian white noise approximation for wave propagation in random media*, SIAM J. Appl. Math. **56** (1996), 1445–1470.

[11] G. Bal, *On the self-averaging of wave energy in random media*, SIAM Multiscale model. simul. **2** (2004), 398–420.

[12] ——, *Kinetics of scalar wave fields in random media*, Wave Motion (2005), 132–157.

[13] G. Bal and Y. Gu, *Limiting models for equations with large random potential; a review*, Comm. Math. Sci. **13** (2013), 729–748.

[14] G. Bal, T. Komorowski, and L. Ryzhik, *Kinetic limits for waves in a random medium*, Kinet. Relat. Models **3** (2010), 529–644.

[15] ——, *Asymptotics of the phase of the solutions of the random Schrödinger equation*, Arch. Rat. Mech. Anal. **100** (2011), 613–664.

[16] G. Bal, G. C. Papanicolaou, and L. Ryzhik, *Diffusive energy scattering from weakly random surfaces*, J. Math. Phys. **40** (1999), 4813–4827.

[17] ——, *Radiative transport limit for the random Schrödinger equation*, Nonlinearity **15** (2002), 513–529.
[18] __________, *Self-averaging in time reversal for the parabolic wave equation*, Stoch. Dyn. **2** (2002), 507–532.

[19] G. Bal and O. Pinaud, *Dynamics of wave scintillation in random media*, Commun. Part. Diff. Eq. **35** (2010), 1176–1235.

[20] __________, *Imaging using transport models for wave-wave correlations*, Math. Mod. Meth. App. Sci. **21** (2011), 1071–1093.

[21] G. Bal and K. Ren, *Transport-based imaging in random media*, SIAM J. Appl. Math. **68** (2008), 1738–1762.

[22] G. Bal and L. Ryzhik, *Wave transport for a scalar model of the love waves*, Wave Motion **36** (2002), 49–66.

[23] __________, *Time splitting for wave equations in random media*, M2AN Math. Model. Numer. Anal. **38** (2004), 961–988.

[24] A. Bamberger, B. Engquist, L. Halpern, and P. Joly, *Parabolic wave equation approximations in heterogeneous media*, SIAM J. Appl. Math. **48** (1988), 99–128.

[25] J.P. Banon, Ø. S. Hetland, and I. Simonsen, *Physics of polarized light scattering from weakly rough dielectric surfaces: Yoneda and brewster scattering phenomena*, Phys. Rev. A **99** (2019), 023834.

[26] W. Bao, S. Jin, and P. A. Markowich, *On time-splitting spectral approximations for the Schrödinger equation in the semiclassical regime*, J. Comp. Phys. **175** (2002), 487–524.

[27] D. Barbolosi, B. Benabdallah, F. Hubert, and F. Verga, *Mathematical and numerical analysis for a model of growing metastatic tumors*, Math. Biosci **218** (2009), 1–14.

[28] R. Bartoszyński, L. Edler, L. Hanin, A. Kopp-Schneider, L. Pavlova, A. Tsodikov, A. Zorin, and A. Yakovlev, *Modeling cancer detection: tumor size as a source of information on unobservable stages of carcinogenesis*, Math. Biosci. **171** (2001), 113–142.

[29] F. G. Bass and I. M. Fuks, *Wave scattering from statistically rough surfaces*, International series in natural philosophy, Elsevier, 2013.

[30] A. Basson and D. Gérard-Varet, *Wall laws for fluid flows at a boundary with random roughness*, Commun. Pure Appl. Math. **61** (2008), 941–987.

[31] M. S. Belen'kii, S. J. Karis, C. L. Osmon, J.M. Brown II, and R. Q. Fugate, *Experimental evidence of the effects of non-kolmogorov turbulence and anisotropy of turbulence*, Proc. SPIE 3749, 18th Congress of the International Commission for Optics (1999).

[32] S. Benzekry, C. Lamont, A. Beheshti, A. Tracz, J.M.L Ebos, L. Hlatky, and P. Hahnfeldt, *Classical Mathematical Models for Description and Prediction of Experimental Tumor Growth*, PLOS Comput. Biol. **10** (2014), e1003800.

[33] S. Benzekry, A. Tracz, M. Mastri, R. Corbelli, D. Barbolosi, and J.M.L. Ebos, *Modeling spontaneous metastasis following surgery: an in vivo-in silico approach*, Cancer Res. **76** (2015).

[34] P. Blomgren, G. C. Papanicolaou, and H. Zhao, *Super-resolution in time-reversal acoustics*, J. Acoust. Soc. Am. **111** (2002), 230–248.

[35] L. Borcea and J. Garnier, *A ghost imaging modality in a random waveguide*, Inverse Probl. **34** (2018), 124007.

[36] L. Borcea, J. Garnier, and K. Sølna, *Wave propagation and imaging in moving random media*, SIAM Multiscale Model. Simul. **17** (2019), 31–67.

[37] L. Borcea, J. Garnier, and C. Tsogka, *A quantitative study of source imaging in random waveguides*, Commun. Math. Sci. **13** (2015), 749–776.