# Some aspects of high frequency wave propagation in randomly perturbed media

Christophe Gomez

<figure><img src="image_1.png" /></figure>

## To cite this version:

Christophe Gomez. Some aspects of high frequency wave propagation in randomly perturbed media. Mathematics [math]. Aix Marseille University, 2024. tel-04612870

HAL Id: tel-04612870

https://hal.science/tel-04612870

Submitted on 14 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Habilitation à diriger des recherches

Spécialité : Mathématiques Appliquées

Christophe Gomez

## Some aspects of high frequency wave propagation in randomly perturbed media

Soutenue publiquement le 12 Juin 2024

### Rapporteurs:

Liliana Borcea, Professeure,
University of Michigan

Tomasz Komorowski, Professeur,
Polish Academy of Sciences,
Institute of Mathematics

George Papanicolaou, Professeur,
Stanford University

### Composition du jury:

Liliana Borcea, Professeure,
University of Michigan

Sonia Fliss, Professeure,
ENSTA - Paris

Josselin Garnier, Professeur,
Ecole Polytechnique

Antoine Gloria, Professeur,
Sorbonne Université

Rémi Rhodes, Professeur,
Aix Marseille Université

Eric Savin, Professeur,
ONERA - CentraleSupelec
Abstract

This manuscript presents certain aspects of high-frequency wave propagation in randomly perturbed media. First, a particular attention is given to the role played by long-range correlations of the medium fluctuations. This property find applications in various areas of physics, engineering, and medical imaging for instance. The interactions between random fluctuations with long-range dependencies and oscillatory behaviors result in intriguing effects, leading to multiscale phenomena. Such multiscale properties contrast with scenarios involving random fluctuations with short-range correlations or mixing properties, where the stochastic effects are observed on a well-defined scale. These multiscale phenomena are investigated within the context of a stochastically forced nonlinear oscillator, the Schrödinger equation, the wave equation, and the radiative transfer equation. These properties are also investigated numerically, using a time-splitting scheme for the Schrödinger equation, and the design of a Monte-Carlo method in the context of radiative transfer.

Secondly, this manuscript describes the role played by boundaries or interfaces in the propagation
domain. Radiative transfer models are discussed for propagation media with boundaries: a half-space,
a slab, and a rectangle are considered. Each of these scenarios exhibits interference phenomena that
affect the energy propagation at specific locations, and for which the geometry of the propagation
media plays a critical role. The role played by rough boundaries and interfaces is discussed under
two scenarios. First, a waveguide model inspired by underwater acoustics is considered, in which
waves evolve along a randomly perturbed ocean surface and uneven bottom topography. Second, the
standard problem of transmission and reflection at an interface is revisited. This time a wave impinges
upon a random interface, and generalized Snell's laws for reflection and transmission are derived.

Finally, an aspect of mathematical biology is discussed, focusing on the stochastic and deterministic modeling of secondary metastatic emission. The relation between these two approaches allow the risk evaluation of a metastatic disease, even when unobservable at the time of cancer diagnosis.
# Résumé

Ce manuscrit présente certains aspects de propagation d'ondes hautes fréquences dans des milieux aléatoires. Dans un premier temps, une attention particulière est accordée au rôle joué par les corrélations à longue portée des fluctuations du milieu de propagation. Cette propriété trouve des applications dans divers domaines en physique, ingénierie et en imagerie médicale par exemple. Les interactions entre des fluctuations aléatoires ayant des dépendances à longue portée et des comportements oscillatoires produisent des effets surprenants, conduisant à des phénomènes multi-échelles. Ces propriétés multi-échelles contrastent avec les comportements effectifs observés pour des fluctuations aléatoires avec des propriétés de mélange ou de corrélations à courte portée. Dans ce dernier cas les effets stochastiques sont observés sur une échelle bien définie. Ces phénomènes multi-échelles sont étudiés dans le contexte d'un oscillateur non linéaire avec un terme source aléatoire, de l'équation de Schrödinger, de l'équation des ondes et de l'équation de transfert radiatif. Ces propriétés sont étudiées numériquement à l'aide un schéma *time-splitting* pour l'équation de Schrödinger et la mise en place d'une méthode de Monte-Carlo appropriée pour le transfert radiatif.

Dans un deuxième temps, ce manuscrit décrit le rôle joué par les frontières, ou interfaces, du domaine de la propagation. Les modèles de transfert radiatif sont abordés pour des milieux de propagation avec des frontières : un demi-espace, une couche et un rectangle. Chacun de ces cas présente des phénomènes d'interférence qui affectent la propagation de l'énergie à des endroits spécifiques, et pour lesquels la géométrie du milieu de propagation joue un rôle essentiel. Le rôle joué par des frontières, ou interfaces, rugueuses est présenté dans deux contextes différents. Tout d'abord, nous présentons un modèle de guide d'ondes inspiré de l'acoustique sous-marine, dans lequel les ondes évoluent le long de la surface de l'océan et d'un fond marin irrégulier. Ensuite, nous revisitons le problème classique de transmission et réflexion d'une onde incidente par une interface aléatoire. Des généralisations des lois de Snell-Descarte pour la réflection et la transmission sont proposées.

Enfin, un aspect de *mathbio*, concernant la modélisation stochastique et déterministe de l'émission secondaire de métastases, est présenté. La relation entre ces deux approches permet d'évaluer le risque métastatique, même lorsque celles-ci ne sont pas observables au moment du diagnostic du cancer.
Acknowledgments

First and foremost, I would like to express my gratitude to the referees of my habilitation thesis, Liliana Borcea, Tomasz Komorowski, and George Papanicolaou. It is a great honor that they took the time to evaluate my work, and I am especially thankful to Liliana for participating in the defense committee in person.

I am equally honored to have Sonia Fliss, Antoine Gloria, Rémi Rhodes, and Eric Savin on the defense committee. Special thanks go to Rémi for accepting the role of local mentor for the habilitation process and handling the associated administrative duties. I also thank Guillaume Bal for his administrative report recognizing the eligibility of my habilitation application.

I sincerely thanks Josselin Garnier, not only for being part of the defense committee but also for having been a very nice Ph.D. advisor. I am also deeply grateful to Lenya Ryzhik for hiring me as a postdoc at Stanford, where I spent two unforgettable years in California. I realize how lucky I was to begin my research career under such excellent conditions with wonderful advisors, both in terms of their human and scientific qualities.

I extend my gratitude to all my collaborators for the stimulating discussions and exciting projects we are working on. Thanks also to my colleagues at I2M, LMA, and the Qarma team for the enjoyable moments during coffee breaks, lunches, scientific-social events, informal or teaching-related discussions,... and even doing some maths sometimes. Special thanks to Caroline Bauzet, for providing me materials to go through all the administrative processes of the habilitation more easily, as well as Joelle Pieret, at the habilitation office, for her administrative flexibility, which allowed me to defend this habilitation on time.

Finally, a special thanks to my wife and kids for their support and our nice everyday life moments.
# Contents

<table><tr><td>Abstract</td><td>1</td></tr><tr><td>Résumé</td><td>2</td></tr><tr><td>Acknowledgments</td><td>3</td></tr><tr><td>List of publications</td><td>6</td></tr><tr><td>Introduction</td><td>7</td></tr><tr><td>1 Random perturbations with long-range correlations and nonlinear oscillators</td><td>16</td></tr><tr><td>1.1 Random fluctuations with long-range correlations</td><td>17</td></tr><tr><td>1.2 A non-central limit theorem</td><td>19</td></tr><tr><td>1.3 A randomly perturbed nonlinear oscillator</td><td>20</td></tr><tr><td>1.3.1 Diffusive limit for a quadratic Hamiltonian</td><td>20</td></tr><tr><td>1.3.2 Diffusive limit for the general Hamiltonian with one non-degenerate critical point</td><td>21</td></tr><tr><td>1.4 Perspectives</td><td>24</td></tr><tr><td>2 Random Schrödinger equation with long-range correlations</td><td>26</td></tr><tr><td>2.1 The random potential</td><td>27</td></tr><tr><td>2.2 Phase modulation</td><td>28</td></tr><tr><td>2.3 The Wigner transform</td><td>29</td></tr><tr><td>2.4 Loss of coherence</td><td>31</td></tr><tr><td>2.5 Radiative transfer regime</td><td>32</td></tr><tr><td>2.6 Asymptotic preserving time-splitting schemes</td><td>33</td></tr><tr><td>2.6.1 The phase modulation</td><td>34</td></tr><tr><td>2.6.2 Loss of coherence</td><td>35</td></tr><tr><td>2.6.3 Radiative transfer regime</td><td>35</td></tr><tr><td>2.7 Perspectives</td><td>36</td></tr><tr><td>3 Wave propagation in random media with long-range correlations</td><td>37</td></tr><tr><td>3.1 Random waveguides</td><td>39</td></tr><tr><td>3.2 Paraxial approximation and fractional Itô-Schrödinger equation</td><td>42</td></tr><tr><td>3.2.1 Paraxial scaling regime</td><td>43</td></tr><tr><td>3.2.2 Main result and properties</td><td>44</td></tr><tr><td>3.2.3 Stochastic integral and fractional Itô-Schrödinger equation</td><td>45</td></tr><tr><td>3.3 Effective fractional paraxial wave equation</td><td>47</td></tr><tr><td>3.3.1 The random fluctuations</td><td>47</td></tr><tr><td>3.3.2 Scaling regime</td><td>48</td></tr><tr><td>3.3.3 Main results</td><td>49</td></tr><tr><td>3.4 Perspectives</td><td>52</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>4</td>
      <td>Radiative transfer with nonintegrable singular scattering kernels</td>
      <td>54</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Regularizing effects</td>
      <td>56</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Diffusion limit</td>
      <td>57</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>A Monte-Carlo methods for radiative transfer with singular kernels in 2D propagation media</td>
      <td>58</td>
    </tr>
    <tr>
      <td>4.3.1</td>
      <td>The ACR method</td>
      <td>59</td>
    </tr>
    <tr>
      <td>4.3.2</td>
      <td>The AS method</td>
      <td>60</td>
    </tr>
    <tr>
      <td>4.3.3</td>
      <td>Numerical Illustrations</td>
      <td>61</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>A Monte-Carlo method for 3D radiative transfer equations with multifractional singular kernels</td>
      <td>64</td>
    </tr>
    <tr>
      <td>4.4.1</td>
      <td>Approximation and probabilistic representation</td>
      <td>65</td>
    </tr>
    <tr>
      <td>4.4.2</td>
      <td>The Monte-Carlo method</td>
      <td>67</td>
    </tr>
    <tr>
      <td>4.4.3</td>
      <td>Numerical illustrations</td>
      <td>68</td>
    </tr>
    <tr>
      <td>4.5</td>
      <td>Perspectives</td>
      <td>73</td>
    </tr>
    <tr>
      <td>5</td>
      <td>Boundary effects and weak localization in radiative transfer</td>
      <td>75</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Boundary effects for a half-space</td>
      <td>76</td>
    </tr>
    <tr>
      <td>5.1.1</td>
      <td>The method of images</td>
      <td>76</td>
    </tr>
    <tr>
      <td>5.1.2</td>
      <td>Radiative transfer model and boundary effects</td>
      <td>78</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Weak localization phenomena for a slab</td>
      <td>79</td>
    </tr>
    <tr>
      <td>5.2.1</td>
      <td>The method of images</td>
      <td>80</td>
    </tr>
    <tr>
      <td>5.2.2</td>
      <td>The self- and cross-Wigner transforms</td>
      <td>81</td>
    </tr>
    <tr>
      <td>5.2.3</td>
      <td>Boundary effects and weak localization phenomena</td>
      <td>81</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>The case of a rectangle</td>
      <td>83</td>
    </tr>
    <tr>
      <td>5.4</td>
      <td>Perspectives</td>
      <td>86</td>
    </tr>
    <tr>
      <td>6</td>
      <td>Wave scattering by rough surfaces</td>
      <td>88</td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>A Pekeris waveguide model with a rough surface and bottom topography</td>
      <td>88</td>
    </tr>
    <tr>
      <td>6.1.1</td>
      <td>Waveguide model</td>
      <td>89</td>
    </tr>
    <tr>
      <td>6.1.2</td>
      <td>The conformal mapping</td>
      <td>91</td>
    </tr>
    <tr>
      <td>6.1.3</td>
      <td>Main result</td>
      <td>92</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Reflection and transmission problems of high-frequency waves through a randomly perturbed interface: generalized Snell's laws</td>
      <td>94</td>
    </tr>
    <tr>
      <td>6.2.1</td>
      <td>The physical model</td>
      <td>95</td>
    </tr>
    <tr>
      <td>6.2.2</td>
      <td>Refraction and transmission for an unperturbed interface</td>
      <td>97</td>
    </tr>
    <tr>
      <td>6.2.3</td>
      <td>Random specular components for $l_c \sim r_0 (\gamma = 1/2)$</td>
      <td>99</td>
    </tr>
    <tr>
      <td>6.2.4</td>
      <td>Effective specular components for $l_c \ll r_0 (\gamma > 1/2)$</td>
      <td>100</td>
    </tr>
    <tr>
      <td>6.2.5</td>
      <td>Incoherent wave fluctuations and generalized Snell's laws for $l_c \ll r_0 (\gamma > 1/2)$</td>
      <td>101</td>
    </tr>
    <tr>
      <td>6.2.6</td>
      <td>Gaussian statistics of the speckle patterns for $l_c \ll r_0 (\gamma > 1/2)$</td>
      <td>105</td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>Perspectives</td>
      <td>106</td>
    </tr>
    <tr>
      <td>7</td>
      <td>Digressions in mathematical biology: Stochastic and deterministic models for the secondary metastatic emission process</td>
      <td>107</td>
    </tr>
    <tr>
      <td>7.1</td>
      <td>A simple probabilistic framework for metastatic risk</td>
      <td>108</td>
    </tr>
    <tr>
      <td>7.2</td>
      <td>Mathematical formalism and results</td>
      <td>108</td>
    </tr>
    <tr>
      <td>7.2.1</td>
      <td>Size-structured model</td>
      <td>108</td>
    </tr>
    <tr>
      <td>7.2.2</td>
      <td>Probabilistic framework for secondary metastatic emission</td>
      <td>110</td>
    </tr>
    <tr>
      <td>7.3</td>
      <td>Numerical illustrations</td>
      <td>111</td>
    </tr>
    <tr>
      <td>7.4</td>
      <td>Perspectives</td>
      <td>112</td>
    </tr>
  </tbody>
</table>
List of publications

[1] C. GOMEZ, *Time-reversal superresolution in random waveguides*, SIAM Multiscale Model. Simul., 7 (2009), pp. 1348–1386

[2] C. GOMEZ, *Wave propagation in shallow-water acoustic random waveguides*, Commun. Math. Sci., 9 (2011), pp. 81–125

[3] C. GOMEZ, *Radiative transport limit for the random Schroödinger equation with long-range correlations*, J. Math. Pures Appl., 98 (2012), pp. 295–327

[4] C. GOMEZ, *Wave decoherence for the random Schroödinger equation with long-range correlations*, Commun. Math. Phys., 320 (2013), pp. 37–71

[5] C. GOMEZ, *Loss of resolution for the time reversal of wave in underwater acoustic random channels*, Math. Mod. Meth. App. Sci., 23 (2013), pp. 2065–2110

[6] C. GOMEZ AND O. PINAUD, *Asymptotics of a time-splitting scheme for the random Schrödinger equation with long-range correlations*, Math. Model. Numer. Anal., 48 (2014), pp. 411–431

[7] C. GOMEZ, *Wave propagation in underwater acoustic waveguides with rough boundaries*, Commun. Math. Sci., 13 (2015), pp. 2005–2052

[8] C. GOMEZ, O. PINAUD, AND L. RYZHIK, *Hypoelliptic estimates in radiative transfer*, Commun. Part. Diff. Eq., 41 (2015), pp. 150–184

[9] C. GOMEZ, O. PINAUD, AND L. RYZHIK, *Radiative transfer with long-range interactions: regularity and asymptotics*, SIAM Multiscale Model. Simul., 15 (2017), pp. 1048–1072

[10] C. GOMEZ AND O. PINAUD, *Fractional white-noise limit and paraxial approximation for waves in random media*, Arch. Rat. Mech. Anal., 226 (2017), pp. 1061–1138

[11] C. GOMEZ AND N. HARTUNG, *Stochastic and deterministic models for the metastatic emission process: formalisms and crosslinks*, Cancer Systems Biology: Methods and Protocols, Springer 2018

[12] C. GOMEZ AND K. SØLNA, *Wave propagation in random waveguide with long-range correlations*, Commun. Math. Sci., 16 (2018), pp. 1557–1596

[13] C. GOMEZ AND O. PINAUD, *Monte Carlo methods for radiative transfer with singular kernels*, SIAM J. Sci. Comput., 40 (2018), pp. A1714–A1741

[14] C. GOMEZ, G.IYER, H. LE, AND A. NOVIKOV, *An oscillator driven by algebraically decorrelating noise*, Commun. Math. Phys., 402 (2023), pp. 231–284

[15] C. GOMEZ AND O. PINAUD, *A Monte Carlo method for 3D radiative transfer equations with multifractional singular kernels*, J. Comput. Phys., 489 (2023), 1122–1149

[16] A. MESSAUDI, R. COTTEREAU AND C. GOMEZ, *Boundary effects in radiative transfer of acoustic waves in a randomly fluctuating half-space*, SIAM Multiscale Model. Simul., 21 (2023), pp. 1299–1321

[17] C. GOMEZ, *An effective fractional paraxial wave equation for wave-fronts in randomly layered media with long-range correlations*, SIAM Multiscale Model. Simul., 21 (2023), pp. 1410–1456

[18] A. MESSAUDI, R. COTTEREAU AND C. GOMEZ, *Weak localization in radiative transfer of acoustic waves in a randomly-fluctuating slab*, submitted (hal-04174463).

[19] C. GOMEZ AND K. SØLNA, *Reflection and transmission problems for high-frequency waves at a randomly perturbed interface: generalized Snell's laws*, submitted (hal-04317227).
# Introduction

The present manuscript gathers all my scientific contributions dating back to my Ph.D. defense. The main topic revolves around the propagation of high-frequency waves in randomly perturbed media under various angles encompassing theoretical and numerical aspects. These researches cover scalar wave propagation, governs by the standard wave equation, as well as the random Schrödinger equation. Although the Schrödinger equation is known as the fundamental equation of quantum mechanics, describing the motion of wave functions at the atomic level, it can also be derived from the standard wave equation to elucidate wave-front propagation for high-frequency collimated waves at the macroscopic level. This latter approximation, known as paraxial (or parabolic) approximation, is widely used in practice as it strongly simplifies the description of physical phenomena at the macroscopic scale, along with the corresponding numerical simulations. In practical aspects the fine scales of variation occurring within the propagation media cannot be described exactly. These fluctuations can have various origins: *impurities* within the propagation medium, salinity or temperature variations, or even geometric perturbations of the medium itself for instance. Therefore, it is natural to treat some propagation media as inherently random. For a given situation, the fluctuations occurring in the propagation medium are deterministic, but they can be considered as a realization of a random process. This manuscript specifically focuses on *weak scattering* type regimes. Roughly speaking, for high-frequency waves, this regime corresponds to propagation media characterized by fast fluctuations occurring at the same rate as the wave frequency, while having low amplitudes. Despite the relatively small amplitudes of these fluctuations, the high-frequency waves and fast variations of the propagation media lead to strong interactions, resulting to significant cumulative stochastic effects on the propagating waves. The *microscopic* description provided by the wave equation at the wavelength scale (which is small for high-frequency waves), or by the Schrödinger equation, is in many cases too rich or too complex to be effectively exploited theoretically or numerically at large macroscopic scales. The derivation of macroscopic models from the microscopic dynamics aims to capture effective features or behaviors that prove useful for applications or theoretical studies on macroscopic scales. By mean of considerations on characteristic scales of the problem at hand, interesting macroscopic models that highlight the main characteristics can be identified. These models offer an effective statistical description of the effects on the propagating waves and can serve as the foundation for imaging techniques among other applications. Importantly, these derived macroscopic models do not depend on specific (unknown in practice) realizations of the medium fluctuations, but rather on their statistical properties through power spectral densities. Various models of power spectra can be found in the physical literature. Well known examples include the Kolmogorov power spectrum for standard atmospheric turbulence, the Pierson-Neuman or Pierson-Moscvitz spectra for the modeling of the swell in underwater acoustics, Gegenbauer scattering kernel in neutronics, and Henyey-Greenstein scattering kernel for imaging through biological tissues, among others.

It is worth mentioning that both the wave equation and the Schrödinger equation can be related to radiative transfer models that describe scattering properties and energy propagation through random media at the macroscopic scale. The concept of radiative transfer is one of the main aspect of this manuscript. The origin of this theory trace back to 1871, when the physicist John William Strutt (lord Rayleigh) conducted pioneering researches on the scattering of light by the atmosphere, shedding light on the explanation for the color of the sky. The theory of radiative transfer as a well-defined mathematical framework took shape in the early 20th century mainly due to the efforts of
Eugen von Lommel (1887), Orest Khvolson (1890), and Arthur Schuster (1905). A. Schuster plays a pivotal role in this theory when, in 1905, he formulated a radiative transfer problem to elucidate the origins of absorption and emission lines observed in stellar spectra [185]. In the 1930s, Subrahmanyan Chandrasekhar made remarkable contributions to the field of radiative transfer by developing a more comprehensive version of the radiative transfer equation that could be applied to a wide range of astrophysical problems [48]. Chandrasekhar's equation incorporates the intricate interplay of scattering, absorption, and emission of radiation within stellar atmospheres. Consequently, radiative transfer theory experienced a surge in prominence among the astrophysics community during the first half of the 20th century. Today, radiative transfer models find applications in many other areas such as neutronics, optics, geophysics, weather forecasting, and even the illumination of scenes in animated movies. Despite the stochastic nature of wave propagation problems, the associated radiative transfer models are deterministic, exhibiting therefore a property called *self-averaging*. When looking at the wave energy, this property results from highly-oscillating random phases that average out, leading to deterministic quantities. Rigorous derivations of the radiative transfer equation from the Schrödinger equation can be found in [18, 66, 92] and for the wave equation in [45, 180].

A significant part of the results in this manuscript is based on approximation-diffusion theorems for random ordinary differential equations or partial differential equations. Stochastic diffusion processes, defined as solutions to stochastic (partial) differential equations, can serve as a macroscopic model to describe effective wave scattering properties. The first results related to approximation-diffusion date back to 1966 with the pioneer works of Rafail Khasminskii [128, 129]. These asymptotic results were further extended to the general context of random fluctuations with mixing properties by George C. Papanicolaou and Werner Kohler in 1974 [171]. Additionally, martingale techniques were employed in this context by George C. Papanicolaou, Daniel W. Stroock and S. R. S Varadhan [173], as martingales and diffusion processes exhibit close connections [197]. This approach has found a wide range of applications in the field of wave scattering as discussed in [70]. The main applications of these models are to devised imaging and inverse problem techniques involving waves in complex media [35–38, 80, 88]. For instance, imaging functionals characterizing object localizations buried in strongly scattering propagation media can be designed by incorporating explicit quantities from macroscopic models accounting for wave scattering. Imaging methodologies have also been proposed by leveraging radiative transfer models [20, 21]. The use of macroscopic models to design these imaging functionals results in efficient inversion methodologies as they accurately represent scattering phenomena responsible for the degradation of the recorded signals at sensors.

There is an extensive literature describing physical systems perturbed through white-noises, Markov processes, or processes with mixing properties. These choices are related to the notion of memory-less perturbations, and under proper assumptions give rise to diffusion phenomena. Despite earlier findings, the notion of long-range dependence, as a distinct phenomena, started gaining recognition in the early 1960’s thanks to a series of paper by Benoit Mandelbrot and James Wallis [146, 147]. These works were spurred by empirical observations from Harold E. Hurst [118, 119] in the 1950’s regarding the water flow in the Nile river. Random fluctuations with long-range dependencies exhibit memory properties and can lead to anomalous diffusion phenomena as well as interesting macroscopic phenomena in physical systems such as for the Schrödinger equation [13, 15, 51, 106], the heat equation [131], turbulent transport [67, 132–135], or the wave equation [153, 154, 191]. This notion of long-range correlations plays also a central role in this manuscript, from Chap. 1 to 4, through the study of nonlinear oscillators, the Schrödinger equation, the wave equation, or radiative transfer models. It unveils interesting properties when studied at various propagation scales, and even within numerical schemes for the random Schrödinger equation. Analyzing the asymptotic behavior of random ordinary differential equations or partial differential equations involving long-range correlations presents inherent technical challenges. The standard martingale approach, typically applied for mixing random fluctuations, may not be suitable for situations involving long-range correlations. In such cases one often resorts to the rough path theory [74, 152] or a moment technique [15, 67] as an alternative approach. Monte-Carlo methods for radiative transfer models, related to randomly perturbed propagation media with long-range correlations, are also presented. The fundamental idea behind this method is to express the solution of the radiative transfer equation as an expectation of an appropriate Markov process, whose infinitesimal generator corresponds to the scattering properties (scattering kernel) of this equation. In the context of long-range correlations, the scattering properties exhibit a singular behavior giving rise to technical challenges. Although energy is scattered in every
directions, the singularity allocates an infinite weight on *infinitely small* changes of directions of the underlying Markov process [96], rendering the standard notion of mean-free path/time invalid.

Scattering of waves by rough surfaces is a fundamental phenomenon with profound implications across various fields of physics and engineering. Problems related to the role played by boundaries or interfaces are discussed in Chap. 5 and 6. The understanding of these phenomena is of significant importance across a wide spectrum of disciplines including optics, solid state physics, remote sensing, radar technology, environmental monitoring, communications, and non-destructive testing among others [29, 57, 120, 169]. The interplay between waves and rough surfaces results in complex phenomena of reflection, transmission, and diffraction, that necessitates the generalization of the standard Snell's laws of refraction and transmission. The influence of random boundaries has also been studied in various physical contexts. For instance, in fluid flows within a medium with random boundaries [30], water wave propagation with a free surface or a random depth [63, 81], and also in wave propagation in underwater acoustics with a perturbed sea surface [61, 137]. While not randomly perturbed, boundaries can also lead to technical difficulties in radiative transfer and interesting interference effects, such as intensity enhancement at boundaries and *weak localization effects* as described in [47, 76], for instance. The derivation of radiative transfer models usually takes place in the full physical space $\mathbb{R}^3$ [12, 45, 180], and some adjustments are required to deal with propagation media involving boundaries. In particular, the main tool to derive radiative transfer equations is the Wigner transform, which is not directly well-suited to address the presence of boundaries. This is a reason why such situations have been less studied in the literature [2, 22, 150, 161, 181]. For instance, in the case of a half-space, the standard radiative transfer equation is equipped with boundary conditions corresponding to a geometric optic type reflection of the energy. This condition is valid for both Dirichlet and Neumann boundary conditions. To account for interference effects at boundaries, the Wigner transform needs to be considered carefully.

One aspect of mathematical biology is also discussed (Chap. 7) while not being entirely disconnected from the ideas discussed for wave propagation. Mathematical models are valuable tools to gain insights into complex biological processes, such as population dynamics, disease dissemination, ecological interactions, and cellular behavior. These models enable quantitative analyses and predictions for biological systems, thus yielding a deeper understanding of these intricate phenomena. Metastasis are responsible for most cancer-related deaths and constitutes a pivotal point in course of the cancer disease [209]. However, knowing that metastases smaller than approximately $10^7$ cells remain undetectable by medical imaging and other diagnostic tools, the clinical absence of visible metastases may not accurately represent a patient's true metastatic status. Consequently, the estimation of metastatic risk in cancer patients without observable metastases holds significant clinical importance [170]. In this context, probabilistic techniques have the potential to derive risk scores from clinical data. For instance, probabilistic representation and Monte-Carlo methods can serve as valuable tools to handle *complex* PDE models, such as the Iwata model for metastatic emission [122]. After model inference from clinical data, the Monte-Carlo approach can be used in assessing forecast uncertainty and for diagnosis purposes.

The remaining of this section is dedicated to providing a chapter by chapter overview of the contents of this manuscript, along with the corresponding main results.

# Chapter 1

The notion of random fluctuations exhibiting long-range correlations plays a central role in this manuscript (from Chap. 1 to 4). A stationary process $v$ is said to have *long-range correlations or dependencies* if its two-point correlation function

$$R(t) := \mathbb{E}[V(t+s)V(s)],$$

for $(t, s) \in \mathbb{R} \times \mathbb{R}$, decays slowly enough at infinity,

$$R(t) \underset{|t| \to \infty}{\sim} \frac{R_0}{|t|^{\gamma}} \quad \gamma \in (0, 1],$$

to not be integrable:

$$\int_0^\infty |R(t)| dt = \infty.$$
This is in contrast with the notion of short-range correlations for which the correlation function is integrable at infinity. As outlined in this manuscript, the mere distinction between long-range and short-range correlations can result in drastically different behaviors in noisy physical systems. Some examples of stochastic processes or random fields are provided in this chapter, and one of them provides a construction that will be used in the following chapters.

The remaining of this chapter describes the long-time behavior of a stochastically forced nonlinear oscillator with 1 degree of freedom

$$ \ddot{x}(t) + f(x(t)) = \varepsilon V(t) \quad t \ge 0. $$

Here $f: \mathbb{R} \to \mathbb{R}$ is a given smooth function, and $V$ is a stochastic process with long-range correlations representing the noise. This situation provides a simple framework highlighting the interactions between random fluctuations with long-range correlations and oscillatory behaviors. The parameter $\varepsilon \ll 1$ represents the amplitude of these random fluctuations. It turns out that the noise itself has a nontrivial behavior over a time scale of order

$$ \frac{1}{\varepsilon^{1/\mathfrak{H}}} \quad \text{with} \quad \mathfrak{H} = 1 - \frac{\gamma}{2} \in (1/2, 1). $$

More precisely, the process defined by

$$ v^{\varepsilon}(t) = \varepsilon \int_{0}^{t/\varepsilon^{1/\mathfrak{H}}} V(s)ds \quad t \ge 0, $$

converges in law, in the space of continuous functions, to a fractional Brownian motion with Hurst index $\mathfrak{H}$ depending on the decay rate of the correlation function $\gamma$. This result is known as non-central limit theorem [152, 202, 203]. A one dimensional standard fractional Brownian motion $W_{\mathfrak{H}}$ with Hurst index $\mathfrak{H}$, on a probability space $(\Omega, \mathcal{T}, \mathbb{P})$, is a centered Gaussian process with covariance function

$$ \mathbb{E}[W_{\mathfrak{H}}(t)W_{\mathfrak{H}}(s)] = \frac{1}{2}(t^{2\mathfrak{H}} + s^{2\mathfrak{H}} - |t-s|^{2\mathfrak{H}}), $$

for any $t, s \ge 0$. The case $\mathfrak{H} = 1/2$ corresponds to a standard Brownian motion.

The non-central limit theorem is in contrast with the standard invariance principle, which relies on the standard central limit theorem for processes with short-range correlations or mixing properties [139, 152]. For the latter the appropriate time scale to observe nontrivial stochastic effects is $1/\varepsilon^2$. Long-range correlations provide nontrivial stochastic effects that become apparent at shorter time scale compared to situations involving short-range correlations or mixing properties:

$$ \frac{1}{\varepsilon^{1/\mathfrak{H}}} \ll \frac{1}{\varepsilon^2} \quad \text{for} \quad \mathfrak{H} \in (1/2, 1). $$

A more general non-central limit theorem is discussed in this chapter. One can naturally wonder how such statistical behavior affects the one of the nonlinear oscillator. It turns out that the energy of the system exhibits nontrivial diffusive effects only at time scales of order $1/\varepsilon^2$ despite fluctuations with long-range correlations, and not at shorter time scales. This scale aligns with the one at which diffusive effects become apparent when considering fluctuations with short-range correlations or mixing properties within the oscillator system. As a result, even though the random perturbations exhibit a degenerate behavior at the time scale $1/\varepsilon^2$, the oscillatory behavior of the system prevents the emergence of pathological behaviors. In the limit $\varepsilon \to 0$, the energy of the system can be described through a standard diffusion process similar (though not identical) to the one of the classical result [71] for a white-noise perturbation. More precisely, considering a smooth Hamiltonian

$$ H(x, y) = \frac{y^2}{2} + \int_{0}^{x} f(s)ds \quad x, y \in \mathbb{R}, $$

with exactly one non-degenerate critical point at $(0, 0)$, which is also the global minimum, and $V$ be a stationary Gaussian noise with long-range correlations, we observe the following convergence in law in $C(0, \infty)$,

$$ H(X^{\varepsilon}) \underset{\varepsilon \to 0}{\Rightarrow} \mathcal{X}. $$
Here, $X^\epsilon(t) = (x(t/\epsilon^2), \dot{x}(t/\epsilon^2))^T$ and $\mathcal{X}$ is a diffusion process driven by a standard Brownian motion (and not a fractional Brownian motion), with infinitesimal generator similar to the one of the standard Freidlin and Wentzell theory with white-noise perturbations of the oscillator [71–73].

This chapter covers the results obtained in [14] of the publication list on page 6, and was conducted during Hai Le’s Ph.D. program under the supervision of A. Novikov at Penn State University.

## Chapter 2

The first situation we present to describe the influence of long-range correlations on a wave propagation problem concerns the Schrödinger equation with a time-dependent random potential. More precisely, we consider the wave function $\phi$ which satisfies

$$i\partial_t\phi + \frac{1}{2}\Delta_x\phi - V(t, \mathbf{x})\phi = 0 \quad t \ge 0, \quad \mathbf{x} \in \mathbb{R}^d,$$

where $V$ exhibits long-range correlations w.r.t. the time variable. The long-time behavior of $\phi$ is studied under the weak scattering regime that consists of small random fluctuations, of order $\sqrt{\epsilon}$, with $\epsilon \ll 1$, as well as the time and space scaling

$$\phi_{\epsilon}(t, \mathbf{x}) = \phi\left(\frac{t}{\epsilon^s}, \frac{\mathbf{x}}{\epsilon^s}\right).$$

The parameter $s > 0$ represents how long and how far the wave function propagates. Under this scaling the Schrödinger equation now reads

$$i\epsilon^s \partial_t \phi_{\epsilon} + \frac{\epsilon^{2s}}{2} \Delta_x \phi_{\epsilon} - \sqrt{\epsilon} V \left( \frac{t}{\epsilon^s}, \frac{\mathbf{x}}{\epsilon^s} \right) \phi_{\epsilon} = 0. \quad (1)$$

To observe non-trivial cumulative stochastic effects on the wave function $\phi_\epsilon$ in the limit $\epsilon \to 0$, an appropriate choice of $s$ needs to be made. Considering short-range correlations or mixing properties for the potential $V$ leads to the unique choice $s = 1$ [15]. In case of long-range correlations the situation is much richer. Equipped with appropriate initial conditions, the latter scaled Schrödinger equation exhibits effective multiscale effects w.r.t. the propagation parameter $s$. In [15] the authors show that for an appropriate choice of $s = s_0$, corresponding to the non-central limit theorem scaling, the wave function exhibits a random phase modulation driven by a unique fractional Brownian motion for all the wavevectors supported by the wave function.

This chapter describes how the wave function is affected when it propagates over longer scales, that is when $s > s_0$. Over such scalings the wave function experiences a highly oscillating random phase modulation that is expected to break the wave function coherence. In this context the main tool to analyze this loss of coherence is based on a properly scaled Wigner transform of the wave function. This Wigner transform relates to the Fourier transform of the correlation function of the wave function at two nearby points. The resulting diffusive behaviors w.r.t. the momentum variable (the Fourier variable of the Wigner transform) quantifies the effective loss of coherence of the wave function. Three different regimes can be exhibited as the propagation parameter $s$ increases. As already mentioned, for $s > s_0$, the wave function produces fast phase modulations that affects significantly the structure of the wave function. One can show that for some $s = s_1 > s_0$ a loss of coherence of the wave field can be quantified through a stochastic fractional heat equation (SFHE) for the Wigner transform. At this stage the Wigner transform analyses the correlations of points that are far apart from each other. In other words, the loss of coherent impacts the large structures of the wave function in a first place. While the SFHE describes how the Wigner transform is impacted, its energy is conserved. Then, as the propagation scale increases $s > s_1$, the finer structures of the wave function are increasingly affected. This time, the random behavior of the SFHE averages out, homogenization phenomena take place, and the resulting equation governing the loss of coherence is a deterministic fractional heat equation. Once $s = 1$ is reached, the loss of coherence is described by a semi-classical limit and a radiative transfer equation. At this scale, the loss of coherence mechanisms are identical whether the random fluctuations exhibit long-range correlations or shot-range correlations.

At the numerical point of view, a time-splitting scheme for the scaled Schrödinger equation (1) can be considered. This numerical scheme provides a discretized solution in time (the spatial variable being handle through a Fourier transformation), which is able to capture, in the statistical sense,
all the mentioned asymptotic regimes. Surprisingly, in addition to all these regimes, the long-range correlations allow a time stepsize independent of $\epsilon$ for most of the configurations.

This chapter covers the results obtained in [3, 4, 6] of the publication list on page 6.

## Chapter 3

The analysis of the impact of random fluctuations with long-range correlations is extended in this chapter to the scalar wave equation for 2D and 3D propagation media. In this context the random fluctuations of the propagation media translate through random fluctuations of the wave-speed profile. This chapter presents some results that extend the ones obtained for 1D propagation media [84,85,154].

Three scenarios are considered in this chapter. The first situation concerns wave propagation in a planar waveguide, which constitutes an intermediate situation between 1D propagation media and propagation in a full space. The type of waveguide considered in this chapter supports a discrete number of modes, which simplifies in some sense the mathematical analysis. The regime under consideration corresponds to the one of a non-central limit theorem and a result similar to [153] is obtained. In this situations the effective stochastic effects lead to phase modulations for each of the propagating modes, but without any coupling between the modes. However, as for the Schrödinger equation, these modulations are driven by the same fractional Brownian motion. This latter characteristic is in contrast to what is observed in [70, Chap. 20] under medium fluctuations with short-range correlations, where the phase modulations are driven by correlated Brownian motions.

The second scenario considers again a regime similar to the non-central limit theorem but for the wave equation over the full ambient space $\mathbb{R}^3$. More precisely, two asymptotic regimes are simultaneously investigated under the same high-frequency limit: the paraxial approximation, where the wave is collimated and propagates along a privileged direction of propagation, and the fractional white-noise limit. It results that the frequency content of the asymptotic pulse front can be described through the fractional Itô-Schrödinger equation

$$d\Psi_{\omega}(z, \mathbf{x}) = \frac{i}{2k_{\omega}}\Delta_{\mathbf{x}}\Psi_{\omega}(z, \mathbf{x}) + ik_{\omega}\Psi_{\omega}(z, \mathbf{x})dW_{\mathfrak{D}}(z, \mathbf{x}) = 0,$$

where $\omega$ is a given frequency, $z$ corresponds to the variable of the main propagation axis, and $W_{\mathfrak{D}}$ to a fractional random field w.r.t. $z$. The variable $\mathbf{x}$ represents the transverse section w.r.t. the main propagation axis. The stochastic integral corresponds here to the fractional equivalent of the Itô-Stratonovich integral. Under short-range correlations a similar result has been obtained in [83] involving a standard Brownian random field and Itô-Stratonovich integral.

The last scenario corresponds this time to a regime similar to a central limit theorem, but involving a randomly layered propagation medium with long-range correlations. This scenario extends the result of [84, 85] for 1D propagation media. In this situation, we observe first a random travel time for the pulse characterized by a fractional Brownian motion that appears to have a standard deviation larger than the pulse width. This is in contrast with the standard O’Doherty-Anstey theory for random propagation media with short-range correlations for which both the standard deviation of the random travel time and the pulse width are of the same order [70, Chapter 8]. This unstable behavior of the travel time under long-range correlations may have a dramatic effect for applications in inverse problems based on travel time estimations, and a deeper understanding of the propagating waves is required. Second, in the present context, the pulse deformation can be characterized by a deterministic paraxial wave equation of the form

$$\partial_{t z}^2 \psi - \frac{c_0}{2} \Delta_{\mathbf{x}} \psi - a_0 D_t^{2+\gamma} \psi = 0 \quad \gamma \in (0, 1),$$

where the $z$-variable corresponds again to the main propagation axis, the $\mathbf{x}$-variable to the transverse section, $t$ to the time variable, $c_0$ to the background wave speed, and $a_0 > 0$ is a constant. Here, $D_t^{2+\gamma}$ stands for the Weyl fractional derivative with respect to time, and whose order depends explicitly on the power decay rate $\gamma \in (0, 1)$ of the correlation function of the medium fluctuations. In the Fourier domain, this equation can be recast as a Schrödinger equation of the form

$$i \frac{\omega}{c_0} \partial_z \tilde{\psi} + \frac{1}{2} \Delta_{\mathbf{x}} \tilde{\psi} + \tilde{a}_0 \omega |\omega|^{1+\gamma} \tilde{\psi} = 0,$$
where $\tilde{a}_0$ is a constant with positive imaginary part. This equation provides a frequency-dependent power-law attenuation of the form

$$|\omega|^\lambda \quad \text{with} \quad \lambda = 1 + \gamma \in (1, 2].$$

Such power law attenuations can be observed in various contexts such as geophysics, laser beam propagation through the atmosphere, or medical imaging for instance. Values for $\lambda$ ranging from 2 to 3 are typical of attenuations in biological tissues [117].

This chapter covers the results obtained in [10, 12, 17] of the publication list on page 6

## Chapter 4

Radiative transfer models have been used for more than a century to describe the propagation of the energy density through random media. When the momentum variable is restricted over the unit sphere $S^{d-1}$, the radiative transfer equation can be written as

$$\partial_t u + \hat{\mathbf{k}} \cdot \nabla_{\mathbf{x}} u = \int_{S^{d-1}} \Phi(|\hat{\mathbf{p}} - \hat{\mathbf{k}}|)(u(\hat{\mathbf{p}}) - u(\hat{\mathbf{k}}))\sigma(d\hat{\mathbf{p}}) \quad (t, \mathbf{x}, \hat{\mathbf{k}}) \in (0, \infty) \times \mathbb{R}^d \times S^{d-1},$$

where $\sigma(d\hat{\mathbf{p}})$ stands for the surface measure on $S^{d-1}$, and $\Phi$ for the scattering kernel. Under fluctuations with long-range correlations this model can be derived from the Schrödinger equation under a semi-classical limit [66]. For the full wave equation such a rigorous derivation has been obtained in [45] for short-range correlations, but it remains an open problem for long-range correlations. Radiative transfer models can be though considered in this situation by just considering scattering kernels with appropriate shapes. The scattering kernel being directly related to the power spectral density of the medium fluctuations, under long-range correlations the associated scattering kernel presents a nonintegrable singularity $\Phi(r) \propto r^{-(d-1+\alpha)}$ with $\alpha \in (0, 2)$. Despite this singularity the radiative transfer equation remains well-defined.

The first part of this chapter consists in presenting qualitative properties of the radiative transfer equation with a nonintegrable singular kernel. We present regularizing effects of such equation through the hypoellipticity property of the transport operator, which implies in particular that the solutions are infinitely differentiable in all variables. This property is obtained through hypoelliptic estimates for a kinetic equation of the form

$$\partial_t u + \hat{\mathbf{k}} \cdot \nabla_{\mathbf{x}} u = (-\Delta_{S^{d-1}})^\beta h \quad (t, \mathbf{x}, \hat{\mathbf{k}}) \in \mathbb{R} \times \mathbb{R}^d \times S^{d-1},$$

where $d \ge 2$, $\beta \ge 0$, $S^{d-1}$ denotes the unit sphere of $\mathbb{R}^d$ and $\Delta_{S^{d-1}}$ the Laplace-Beltrami operator on $S^{d-1}$. Assuming some fractional Sobolev regularity in the momentum variable $\hat{\mathbf{k}} \in S^{d-1}$, we obtain estimates for fractional derivatives of $u$ w.r.t the $(t, \mathbf{x})$-variables whose order depends on $\beta$ and the assumed regularity. The smoothness of $u$ is then obtained using a bootstrap argument. The diffusion limit is proved using probabilistic techniques as in the case of a regular scattering operator. The diffusion coefficient is therefore non-zero and finite.

The second part of this chapter is devoted to Monte-Carlo methods for radiative transfer equations with non-integrable singular scattering kernels. As opposed to the case where the scattering cross section is integrable and leads to a non-zero mean free time, the cross section is no longer integrable in the long-range situation and yields a vanishing mean free time. This gives rise to numerical difficulties as standard Monte-Carlo methods based on a naive regularization, by just introducing a cutoff around the singularity, exhibit large jump intensities and an increased computational cost to reach a proper accuracy. A particular care is then required when constructing the stochastic processes used in the Monte-Carlo methods. We propose a method inspired by the finance literature, and introduced by Asmussen-Rosiński [9] and Cohen-Rosiński [53]. This method is based on a small jumps - large jumps decomposition allowing us to treat the small jumps efficiently and reduce the computational burden.

This chapter covers the results obtained in [8, 9, 13, 15] of the publication list on page 6.

## Chapter 5

This chapter concerns as well some aspects of radiative transfer arising from acoustic wave propagation, but from a different perspective as the one discussed in the previous chapter. Rather than
investigating the effects produced by a singular scattering kernel, we discuss the role played by boundaries of the propagation medium. We present three scenarios: a half-space, a slab, and a rectangle. This chapter covers the results in [16,18] of the publication list on page 6, which are parts of Adel Messaoudi's Ph.D thesis conducted under the supervision of Régis Cotiereau (Laboratory of Mechanics and Acoustics in Marseille) and myself.

For each scenario two distinct aspects are discussed. While the standard derivations of radiative transfer models from the wave equation usually take place in the full physical space $\mathbb{R}^3$ [12, 45, 180], some adjustments are required to deal with medium involving boundaries. In particular, the main tool to derive radiative transfer equations is the Wigner transform. This tool is not directly well-suited to address the presence of boundaries. In this chapter, for each scenario, the wave propagation problem is extended to the full space thanks to the method of images together with a proper periodization of the medium parameters and the source. Such an extension allows the use of the standard Wigner transform and enable a standard asymptotic analysis [12,180].

The other aspect discussed in this chapter concerns interference effects resulting from reverberations at the boundaries. In a first time, we discuss the case of the half-space where boundary effects can be observed within one wavelength along the boundary. This results in a doubling of the intensity for Neumann boundary conditions, and a canceling of intensity for Dirichlet boundary conditions. In a second time, the situation of a slab is discussed. As for the half-space, interference effects are obtained within one wavelength of the two slab boundaries. However, additional interference effects are observed within one wavelength along two parallel plans w.r.t to the boundaries, and passing through the source location and one symmetric point. These extra effects are referred to as *weak localization* phenomena in the sense of [47, 76]. These effects result in intensity enhancement or canceling (depending on the boundary conditions) of the coherent energy only, not the one carried by the coda (the multiply scattered signal). The underlying idea of *weak localization* is the following. When coherent waves propagate through a medium delimited by boundaries, they follow different paths and bounce on these boundaries. However, some of these paths are the reciprocal versions of others, meaning that the waves can traverse a similar path in opposite directions. This symmetry between paths results in constructive interference effects within the considered structure. Finally, the case of a rectangle is discussed where radiative transfer models are related to non destructive testing applications [47, 76]. As for the slab, both types of interference effects can be observed within one wavelength of the boundaries, and along lines parallel to the rectangle boundaries passing through the source location and three symmetric points. These lines are fixed over time. In this context, additional interference effects can be observed. Constructive interferences affecting the whole energy (coherent and coda) can be observed along lines passing through the corners, and interference effects, affecting only the coherent part of the energy, appear along oblique lines passing through the source location and eleven symmetric points inside the rectangle. In both cases, these lines support interference effects only once over time.

# Chapter 6

This chapter discusses two distinct scenarios. First, we discuss the situation where the waves propagate along a waveguide with rough-boundaries, and second, the case where waves impinge upon a rough surface. In both scenarios, the waves interact with the interface in different manners, leading to distinct mathematical approaches.

In underwater acoustic the role of ocean swell or variations of the seabed topography can have a significant impact on acoustic signals. Wave propagation in waveguides with rough boundaries has been studied for a long time due to its wide range of applications, but mainly motivated by submarine detection and telecommunication [126, 138]. In this context, a pressure field can be decomposed over three kinds of mode: the propagating modes, which travel over long distances along the waveguide; the radiating modes, able to propagate deeply into the ocean bottom; and the evanescent modes, which decay exponentially w.r.t to the propagation distance along the waveguide. We describe in the first part of this chapter the effects produced by a randomly perturbed free surface and an uneven bottom topography on the propagating modes. Using an asymptotic analysis based on a separation of scales technique and an approximation-diffusion theorem, an asymptotic form of the distribution for the forward mode amplitudes can be obtained. This distribution allows to characterize effective attenuations on the propagating-mode amplitudes induced by the surface and bottom fluctuations.
It appears that both the surface and bottom fluctuations affect the propagating-mode amplitudes mainly in similar manners. However, this attenuation is stronger for the highest propagating modes, with losses into the ocean bottom, due to a strong coupling with the radiating modes.

The second scenario discussed in this chapter corresponds to the standard problem of reflexion and transmission of a high-frequency pulse at a rapidly oscillating rough interface with general mixing properties. Under the paraxial (parabolic) scaling, the specular and speckle (diffusive) components of the reflected and transmitted signals are precisely characterized. The specular components correspond to the reflected and transmitted components resulting from the standard laws of reflexion and transmission for a flat (unperturbed) interface, and producing what we refer to as specular cones. The speckles correspond to incoherent (random) wave fluctuations resulting from diffraction and mainly observed away from the specular cones. A *critically scaled* interface is considered, in the sense that the amplitudes of the interface fluctuations and the central wavelength are of the same order. In this context, if the correlation length of the interface fluctuations and the beam width are of the same order, random specular components are observed, but no speckle component. The reflected and transmitted fields are confined to the specular cones. The situation with a correlation length smaller than the beam width leads to homogenization effects providing deterministic specular components similar to the case of a flat interface, but with effective initial conditions accounting for the scattering effects. However, in this situation, there are also relatively broad cones (containing the specular cones) where the wavefields form speckle patterns. The width of these cones is characterized by the ratio $\lambda/l_c$, where $\lambda$ is the central wavelength and $l_c$ the correlation length of the interface fluctuations. The two-point correlation functions for these speckle patterns are presented and exhibit self-averaging properties. Furthermore, we present a central limit theorem type result for the speckle patterns, showing that they can be modeled as Gaussian random fields. These descriptions allow the derivation of generalized Snell's laws of refraction and transmission depending on an effective scattering operator at the interface.

This chapter covers the results obtained in [7, 19] of the publication list on page 6.

# Chapter 7

This chapter presents an aspect of mathematical biology related to metastatic spreading. This work has been conducted in collaboration with Niklas Hartung during his Ph.D. program, at the Institut of Mathematics of Marseille, under the supervision of Guillemette Chapuisat and Florence Hubert (see [11] in the publication list on page 6).

Although the detection of metastases radically influences the prognosis and treatment decisions for cancer patients, the presence of clinically undetectable micrometastases hampers the consistent classification of the disease as localized or metastatic. This chapter discusses mathematical modeling efforts that could help to estimate the metastatic risk in such a scenario. We focus on two approaches both accounting for secondary metastatic emissions, that is the ability for metastases to emit some metastases themselves. The first one relies on a deterministic framework to describe the micrometastatic state using a size-structured density function in a partial differential equation model. The second approach is based on a stochastic framework to describe metastatic emission events occurring at random times formalized through a cascade of Poisson processes. We highlight an inherent crosslink between the stochastic and deterministic frameworks and discuss its implication for assessing metastatic risk post-surgery.

The remaining of the manuscript consists in providing more detail presentations of the results outlined in this introduction.
Chapter 1

Random perturbations with long-range
correlations and nonlinear oscillators

It is well known that a properly scaled stochastic process with long-range correlations (or dependen-
cies) converges to a fractional Brownian motion (fBm), whose Hurst index is determined by the rate
at which the correlation function decays [152, 202, 203]. This result is known as non-central limit
theorem. When the Hurst parameter is not 1/2, the normalized limit of the noise has memory prop-
erties and is a non-Markovian process. An interesting example to illustrate how oscillatory behaviors
effectively erase memory effects of random fluctuations is a stochastically forced nonlinear oscillator
with 1 degree of freedom:

$$
\ddot{x}(t) + f(x(t)) = \epsilon V(t), \quad x_0 \in \mathbb{R}, \quad \dot{x}_0 = y_0 \in \mathbb{R}. \tag{1.1}
$$

Here, $f: \mathbb{R} \rightarrow \mathbb{R}$ is a given smooth function, and $V$ is a stochastic process with long-range correlations as defined precisely below. Our interest is here to study the asymptotic long-time behavior of $x$. This chapter is based on the results obtained in [14] of the publication list on page 6.

Although the rescaled noise converges to a fBm, a process that exhibits memory properties, the oscillatory behavior of the system has the ability to affect the memory properties of the random fluctuations, resulting in an effective diffusive behavior with no memory. The aim of this chapter is to provide a precise description of how the oscillatory behavior is able to erase all the memory of the random fluctuations. A similar question was investigated in [133] for a passive tracer advected by a periodic shear flow. In this case, it appears that there exists a parameter regime where the time-rescaled dynamics is Markovian, and the memory effect of the noise is erased. However, there is also a regime (namely the *very long-time* behavior when the Hurst index of the driving noise is larger than 1/2) where the memory effect persists. In contrast, for the oscillator (1.1), the memory effect *never* persists, and the effective long-time statistical behavior is always Markovian. The main reason is that the oscillatory nature of the deterministic dynamics counterbalances the slowly decaying correlations, and kills the memory effects. To study (1.1) we cannot use directly the limit theorem for additive functionals of fBm’s used by [133]. Instead, we recast (1.1) as a fast-slow system and use an approximation-diffusion framework.

To recast (1.1) as a stochastically perturbed Hamiltonian system, we first introduce the Hamilto-
nian

$$
H(x, y) := \frac{1}{2} y^2 + \int_0^x f(s) \, ds, \tag{1.2}
$$

so that, setting $X(t) := (x(t), y(t))^T = (x(t), \dot{x}(t))^T$, we have

$$
\dot{X}(t) = \nabla^{\perp} H(X(t)) + \varepsilon V(t) \mathbf{e}_2, \quad X(0) = \begin{pmatrix} x_0 \\ y_0 \end{pmatrix} \in \mathbb{R}^2, \qquad (1.3)
$$

with

$$
\nabla^{\perp} := \begin{pmatrix} \partial_y \\ -\partial_x \end{pmatrix}, \quad \text{and} \quad \mathbf{e}_2 := \begin{pmatrix} 0 \\ 1 \end{pmatrix}.
$$
Here, the Hamiltonian $H$ is assumed to have exactly one non-degenerate critical point. To study the long-time behavior of this system, we consider the time rescaled process

$$X^{\epsilon}(t) := X(t/\epsilon^2),$$

satisfying

$$\dot{X}^{\epsilon}(t) = \frac{1}{\epsilon^2} \nabla^{\perp} H(X^{\epsilon}(t)) + \frac{1}{\epsilon} V\left(\frac{t}{\epsilon^2}\right) \mathbf{e}_2, \quad X_0^{\epsilon} = X(0) = \begin{pmatrix} x_0 \\ y_0 \end{pmatrix} \in \mathbb{R}^2. \quad (1.4)$$

In the absence of noise (meaning that $V \equiv 0$), the process $X^\epsilon$ travels very fast along the level sets of the Hamiltonian. When $V$ is a white-noise, the asymptotic behavior of $X^\epsilon$ is described through the averaging principle of Freidlin and Wentzell [71–73] by a diffusion process across these level sets. To capture this limiting behavior, the fast motion is filtered out by projecting $X^\epsilon$ onto the *Reeb graph* of the Hamiltonian. This has the effect of identifying all closed trajectories of the Hamiltonian system, where the fast motions take place, into single points. In this context it is convenient to describe $X^\epsilon$ through action-angle coordinates. The angular coordinate of $X^\epsilon$ evolves very fast providing no meaningful limit as $\epsilon \to 0$. On the other hand, the action coordinate of $X^\epsilon$ exhibits a non-trivial limit as $\epsilon \to 0$, which results from the interaction between the noise and the averaged angular coordinate. To study this behavior, the Hamiltonian itself is used as a proxy for the action coordinate and the convergence is obtained for $(H(X^\epsilon))_\epsilon$.

This chapter is organized as follows. Some examples of stochastic processes with long-range correlations are provided in Sect. 1.1. One them is used in several chapters of this manuscript. Sect. 1.2 describes how long-range dependencies are related to memory effects through fBm's and a *non-central* limit theorem. Finally, Sect. 1.3 describes how the memory properties of the noise are erased by the oscillatory behavior of the Hamiltonian system. The case of the quadratic Hamiltonian ($f(x) = x$) is presented as it provides explicit formulations and a simple understanding of the phenomena. For general Hamiltonian (1.2), action-angle coordinates are introduced as well as their corresponding approximation-diffusion theorem. From this result, the asymptotic diffusive behavior of $(H(X^\epsilon))_\epsilon$ is derived.

## 1.1 Random fluctuations with long-range correlations

A stationary random process $V$ is said to have long-correlations if its two-point correlation function

$$R(t) := \mathbb{E}[V(t+t')V(t)] \quad (t, t') \in \mathbb{R} \times \mathbb{R}$$

decays slowly enough at infinity, for instance

$$R(t) \underset{|t| \to \infty}{\sim} \frac{R_0}{|t|^{\gamma}} \quad \text{with} \quad \gamma \in (0, 1], \qquad (1.5)$$

to not be integrable, in the sense that

$$\int_0^\infty |R(t)| dt = \infty.$$

There exists several basic examples of continuous-time processes with long-range correlations that can be found in the literature (see [84, 154] for instance). Two simple examples are the following: the fractional white-noise model corresponding to

$$V(t) = W_{\mathfrak{H}}(t) - W_{\mathfrak{H}}(t+l_c),$$

and the fractional Ornstein-Uhlenbeck model for which

$$V(t) = W_{\mathfrak{H}}(t) - \frac{1}{l_c} \int_{-\infty}^{z} e^{(y-t)/l_c} W_{\mathfrak{H}}(y) dy.$$

In these two examples, $W_{\mathfrak{H}}$ is a fBm with Hurst index $\mathfrak{H} \in (1/2, 1)$, and the asymptotic behavior at infinity of their correlation functions is given by (1.5) with

$$\gamma = 2 - 2\mathfrak{H}, \quad \text{and} \quad R_0 = \mathfrak{H}(2\mathfrak{H} - 1)l_c^2.$$
More sophisticated models can also be found in [154]. These two Gaussian models are convenient for an analysis based on moment techniques [15,66,99] or the rough-path theory [152,154]. When studying the asymptotic behaviors of certain randomly perturbed physical systems arising from the Schrödinger equation or the wave equation for instance, it is not always clear how these two approaches can be effectively applied. While the application of the moment technique may rely on purely technical aspects, the rough-path approach has to face the inherent infinite-dimensional nature of the functional spaces to which the solutions belong. Additionally the absence of a Hilbert-Schmidt structure in the PDE's further complicates the process of reducing the original problem to a finite-dimensional one. This latter aspect needs further investigations. Even for the simple oscillator presented in this chapter, a nonlinear Hamiltonian gives rise to an oscillatory behavior that needs to be expended over the infinite basis of Fourier modes. The resulting structure of the system involves a nonlinear relationship between the noise V with a degenerate term. This nonlinear relationship brings technical difficulties to apply both the moment technique or the rough-path approach. To address these technical difficulties, a specific definition of the random fluctuations is considered, exhibiting long-range correlations, but allowing the use of the perturbed-test-function method and the martingale approach [139]. These approaches can be applied together to handle non-necessarily Markovian noise, such as *mixing processes*. The noise we consider to apply this approach is obtained, as in [15], by super-imposing Ornstein-Uhlenbeck type processes written in the form

$$V(t) = \int_{-\infty}^{t} \int_{S} e^{-\mu|p|^{2\beta}(t-u)} \mathcal{B}(du, dp). \qquad (1.6)$$

Here, $S \subseteq \mathbb{R}$ is a bounded symmetric open interval containing 0. Also, $\mu$ and $\beta$ are positive constants, and $\mathcal{B}$ is a Gaussian random measure that is white in time and colored in space:

$$\mathcal{B}(du, dp) = \sqrt{2\mu r(p)} |p|^{\beta} \mathbf{1}_{S}(p) \xi(du, dp),$$

where $\xi$ is a 2D Gaussian white noise, and $r: S \setminus \{0\} \to [0, \infty)$ is defined by

$$r(p) = \frac{\lambda(p)}{|p|^{2\alpha}}.$$

Here, $\lambda: S \to \mathbb{R}_+$ is a smooth bounded even function such that $\lambda(0) \neq 0$. The correlation function of $V$ is therefore given by

$$R(t) = \int_S r(p) e^{-\mu|p|^{2\beta}|t|} dp, \qquad (1.7)$$

so that the condition

$$\int_S r(p) dp \in (0, \infty),$$

that is

$$\alpha < \frac{1}{2},$$

is required for this correlation function to be well defined. This construction gives rise to the limit

$$\lim_{t \to \infty} t^{\gamma} R(t) = R_0,$$

where $R_0$ and $\gamma$ are defined by

$$\gamma = \frac{1 - 2\alpha}{2\beta} \quad \text{and} \quad R_0 = \lambda(0) \int \frac{e^{-\mu|p|^{2\beta}}}{|p|^{2\alpha}} dp.$$

Long-range correlations for $V$ are therefore obtained by assuming

$$2(\alpha + \beta) \geq 1.$$

In Fig. 1.1, we illustrate the difference of statistical behavior between long-range correlations and short-range correlations, the latter having an integrable two-point correlation function. From these pictures, one can observe that slowly decaying correlations produce longer excursions of the random trajectories, due to the persistence of the correlations, than for rapidly decaying correlations. In this latter case, the trajectories cannot really produce correlation patterns and look almost like the ones of a white-noise.
<figure><img src="image_3.png" /><figcaption>Fig. 1.1: Illustration of three realizations of the random process defined by (1.6) with long-range correlations for the left-picture, and short-range correlations for the right-picture. Here, $\mu = 1$, $\alpha = 1/4$, $a(p) = 1_{(-10,10)}(p)$, $\beta = 1/2$ (that is $\gamma = 1/2$) for the left-picture illustrating the long-range correlations, and $\beta = 1/6$ (that is $\gamma = 3/2$) for the right-picture illustrating short-range correlations.</figcaption></figure>

## 1.2 A non-central limit theorem

It is known for quite some time that under long-range correlations the normalized integrated noise converges to a fBm [152,202,203]. A version of this result, known as non-central limit theorem, is the following.

**Proposition 1.2.1** Let $V$ be a stationary Gaussian process with covariance function $R$ given by

$$R(t) = \frac{L(t)}{t^{\gamma}} \qquad (1.8)$$

for some $\gamma \in (0,2)$, and a slowly varying function $L$ at infinity<fnref n="1" />. If $\gamma \in (1,2)$, we additionally suppose

$$\int_{0}^{\infty} R(t) dt = 0. \qquad (1.9)$$

Let

$$\sigma(\varepsilon) := \begin{cases} L(\varepsilon^{-2})^{1/2} \varepsilon^{\gamma} & \text{if } \gamma \neq 1, \\ L(\varepsilon^{-2})^{1/2} \varepsilon |\ln(\varepsilon)|^{1/2} & \text{if } \gamma = 1, \end{cases} \quad \text{and} \quad u^{\varepsilon}(t) := \frac{1}{\sigma(\varepsilon)} \int_{0}^{t} V\left(\frac{s}{\varepsilon^2}\right) ds.$$

Then, as $\varepsilon \to 0$, the family of processes $(u^{\varepsilon})_{\varepsilon}$ converges in law in $C(0, \infty)$ to $\sigma_{\mathfrak{H}}W_{\mathfrak{H}}$, where $W_{\mathfrak{H}}$ is a fractional Brownian motion with Hurst index $\mathfrak{H} = 1 - \gamma/2$, and

$$\sigma_{\mathfrak{H}}^2 := \begin{cases} \frac{1}{\mathfrak{H}|2\mathfrak{H} - 1|} & \text{if } \mathfrak{H} \neq \frac{1}{2}, \\ \frac{1}{\mathfrak{H}} & \text{if } \mathfrak{H} = \frac{1}{2}. \end{cases}$$

Note that despite short-range correlations, the additional assumption (1.9) leads to a negatively correlated fBm with Hurst index $H < 1/2$. Note also that (1.9) is not compatible with the definition (1.6) in view of (1.7). Without this additional assumption, a standard invariance principle would applied with the scaling $\sigma(\varepsilon) = L(\varepsilon^{-2})^{1/2}\varepsilon$, and a standard Brownian motion as limiting process.

A more refined version of this result can be found in [154]. Under additional assumptions the limiting process can be a general Hermite process, in particular with non-Gaussian statistics. This is in contrast to processes with short-range correlations (or mixing properties) and the standard invariance principle where the limits inherently exhibit Gaussian statistics. As described below, when introducing oscillatory behaviors, the memory effects of the random fluctuations are lost and standard diffusion phenomena (memoryless Markov processes) emerge despite the memory effects of the system perturbations.

Note that the above result, is not exactly what is stated in the introduction. However, both are connected through the change of variable $\varepsilon \to \varepsilon^{1/(2-\gamma)}$.

<footnote n="1">¹A function $L$ is said to be slowly varying at infinity if for every $s > 0$ we have $L(st)/t \to 1$ as $t \to \infty$.</footnote>
## 1.3 A randomly perturbed nonlinear oscillator

The remaining of this chapter discusses two results showing that even though the driving noise has memory properties, the time correlations are erased by the oscillatory dynamics of (1.4), and $(H(X^\epsilon))_\epsilon$ converges to a diffusion driven by standard Brownian motion (without memory). Two scenarios are discussed below. The first one concerns a quadratic Hamiltonian, for which explicit formulations can be carried out, highlighting explicitly the role of the oscillations in the loss of memory properties. In this particular case, the limiting diffusion is a scaled version of a 2D square Bessel process. The second scenario extends this result to more general smooth Hamiltonian with exactly one non-degenerate critical point.

### 1.3.1 Diffusive limit for a quadratic Hamiltonian

For a quadratic Hamiltonian the level sets are simply circles centered at 0 and whose radii depend on the associated energy. This simple characterization of the level sets provides a simple description of the role played by the oscillations of the oscillator. The first result is as follows.

**Proposition 1.3.1** Let $H$ be the quadratic Hamiltonian,

$$H(X) = \frac{|X|^2}{2}, \qquad (1.10)$$

for $X = (x, y) \in \mathbb{R}^2$, and suppose that the noise $V$ is a stationary Gaussian process whose covariance function is of the form (1.8). For $\gamma \in (0, 1]$, we further assume that $L$ has the slow increase property

$$\lim_{t \to \infty} \frac{L'(t)}{L(t)/t} = 0.$$

Then, the family $(H(X^\epsilon))_\epsilon$ converges in law to $H(W)$, as $\epsilon \to 0$, where $W$ is a 2D Brownian motion with $W_0 = X_0$ and $2 \times 2$ covariance matrix

$$D := \begin{pmatrix} \Gamma_c & 0 \\ 0 & \Gamma_c \end{pmatrix}, \quad \text{with} \quad \Gamma_c := \int_0^\infty R(t) \cos(t) \, dt. \qquad (1.11)$$

It is worth noticing that, even for $\gamma \in (0, 1]$, the coefficient $\Gamma_c$ in (1.11) is well-defined since the cos function is periodic and mean zero. This resulting cos function is due to the quadratic shape of the Hamiltonian (1.10) and the resulting circular level sets. As a result, despite the nonintegrability of the correlation function and the memory effects, the diffusive (memoryless) limit is valid thanks to the oscillatory behavior of the system.

The Hamiltonian $H$ being here rotationally invariant, Prop. 1.3.1 is based on a time-dependent rotation in space that prove convergence of $(H(X^\epsilon))_\epsilon$ by studying integrals of the form

$$\frac{1}{\epsilon} \int_{0}^{t} V\left(\frac{s}{\epsilon^{2}}\right) \cos\left(\frac{s}{\epsilon^{2}}\right) d s \quad \text { and } \quad \frac{1}{\epsilon} \int_{0}^{t} V\left(\frac{s}{\epsilon^{2}}\right) \sin\left(\frac{s}{\epsilon^{2}}\right) d s.$$

At the limit these processes provide a 2D standard Brownian motion with diagonal covariance matrix $D$. This result has also been obtained in [152] under slightly different assumptions, and generalized in [155].

To reiterate our main point, note that Prop. 1.3.1 corresponds to

$$H\left(X_0 + \frac{1}{\epsilon^2} \int_0^t \nabla^\perp H(X^\epsilon(s)) ds + \frac{1}{\epsilon} \int_0^t V\left(\frac{s}{\epsilon^2}\right) e_2 ds\right) \underset{\epsilon \to 0}{\Longrightarrow} H(W(t)), \qquad (1.12)$$

in law in $C(0, \infty)$. However, from Prop. 1.2.1, the noise term satisfies

$$\frac{1}{\sigma(\epsilon)} \int_{0}^{t} V\left(\frac{s}{\epsilon^{2}}\right) d s \underset{\epsilon \to 0}{\Longrightarrow} \sigma_{\mathfrak{S}} W_{\mathfrak{S}}(t),$$

for a fBm $W_{\mathfrak{S}}$ with Hurst index $\mathfrak{S}$. Interestingly, while $H(W)$ is a Markov process with no memory, $W_{\mathfrak{S}}$ is a non-Markovian process with memory. Note further that:
• When $\gamma \in (0, 1]$, we have $\varepsilon \ll \sigma(\varepsilon)$ and so the term $\frac{1}{\varepsilon} \int_0^t V(s/\varepsilon^2) ds$ appearing in (1.12) diverges as $\varepsilon \to 0$.

• On the other hand, when $\gamma \in (1, 2)$, together with (1.9), we have this time $\sigma(\varepsilon) \ll \varepsilon$, and so the term $\frac{1}{\varepsilon} \int_0^t V(s/\varepsilon^2) ds$ vanishes.

In both cases, the oscillatory behavior of the Hamiltonian system contributes non-trivially by preventing from a blow up of the fluctuations when $\gamma \in (0, 1]$, or enhancing the effects of the fluctuations for $\gamma \in (1, 2)$.

### 1.3.2 Diffusive limit for the general Hamiltonian with one non-degenerate critical point

The case where the Hamiltonian (1.2) has exactly one non-degenerate critical point, but is not necessarily quadratic, is now discussed. As the orbits of the Hamiltonian are more complex, a more involved analysis is required. Using the action-angle coordinates, the role of the oscillatory behavior can still be highlighted though the angle variable, but in a more complex framework. In this case, for technical reasons discussed in Sect. 1.1, we need to work with the stationary Gaussian noise given by (1.6). The asymptotic behavior of $(H(X^\varepsilon))_\varepsilon$ can be described as follows.

**Theorem 1.3.1** Let $H: \mathbb{R}^2 \to \mathbb{R}$ be a smooth Hamiltonian with exactly one non-degenerate critical point at $(0, 0)$, which satisfies (1.19) and (1.20) under the action-angle coordinates introduced below, and let $V$ be the stationary Gaussian noise given by (1.6). Then

$$H(X^\varepsilon) \underset{\varepsilon \to 0}{\longrightarrow} \mathcal{X}, \qquad (1.13)$$

in law in $C(0, \infty)$, where $\mathcal{X}$ is a diffusion process driven by a standard Brownian motion, and infinitesimal generator

$$L := \frac{1}{2\Lambda(\mathcal{X})} \partial_{\mathcal{X}}(\Sigma(\mathcal{X})\partial_{\mathcal{X}}). \qquad (1.14)$$

Here, the coefficients $\Lambda$ and $\Sigma$ are defined by

$$\Lambda(\mathcal{X}) := \oint_{\{H=\mathcal{X}\}} \frac{dl}{|\nabla H|},$$

and

$$\Sigma(\mathcal{X}) := 2 \int_0^\infty ds R(s) \oint_{\{H=\mathcal{X}\}} \partial_y H(\tilde{X}^{x,y}(s)) \partial_y H(x, y) \frac{dl(x, y)}{|\nabla H(x, y)|},$$

where $\tilde{X}^{x,y}$ is the solution to the unperturbed Hamiltonian system:

$$\partial_t \tilde{X}^{x,y}(t) = \nabla^\perp H(\tilde{X}^{x,y}(t)), \quad \text{with} \quad \tilde{X}^{x,y}(t=0) = (x, y).$$

Thm. 1.3.1 follows immediately from the more general Thm. 1.3.2 below concerning a fast-slow system.

Let us make some comments on this result before diving into the ideas behind the proof. First of all, the main difference between Thm. 1.3.1, and the standard Freidlin-Wentzell theory [71] is in the coefficient $\Sigma$. In the present context, $\Sigma$ involves an average of the term $\partial_y H(\tilde{X}^{x,y}(s))\partial_y H(x, y)$ along the level sets, with a shift of one factor to the point $\tilde{X}^{x,y}(s)$, while in [71] the coefficient $\Sigma$ depends only on $(\partial_y H(x, y))^2$. This shifted average term plays the role of the cos function in (1.11) for the quadratic Hamiltonian, so that the diffusion coefficient is still well defined despite long-range correlations of the noise term. The reason for this modification, compared to [71], lies in the fact that the oscillations of the Hamiltonian system and the random fluctuations in (1.4) are strongly coupled in our case, they are both evolving on the scale $1/\varepsilon^2$, and fully interact with each. This strong interaction produces an effective diffusion coefficient involving the correlation of pairs of points belonging to the same orbit: $\partial_y H(\tilde{X}^{x,y}(s))$ and $\partial_y H(x, y)$. If the oscillatory behavior of the Hamiltonian system is slower than the random fluctuations, then their coupling is weak and the same generator as in [71] would be obtained. This dichotomy of limiting behaviors between strongly and weakly coupled setups is well-known in other approximation-diffusion limits of random differential equations with periodic components, (for
instance Thm. 6.4 and 6.5 in [70] can be compared). However, in the present chapter, the oscillations of the Hamiltonian system need to behave at the same scale as the random fluctuations, in other words the oscillations need to be strong enough, for the coefficient Σ to be well-defined.

When $H$ is not quadratic, the convergence of $H(X^\epsilon)$ cannot be reduced to Gaussian integrals as in Prop. 1.3.1, and so the proof of Thm. 1.3.1 is more involved. One way to analyze integrable Hamiltonian systems is to use a set of *action-angle coordinates*. These coordinates separate the slow and fast motion, and preserve the Hamiltonian structure. The Liouville-Arnold theorem [8, Sect. 50 pp. 279] asserts that there exists a symplectic canonical transformation $\varphi: X = (x, y) \mapsto (I, \theta) \in \mathbb{R}\times\mathbb{T}$, where the action variable $I$ and the angle variable $\theta$ satisfy

$$K(I) = H(x, y) \quad \text{and} \quad \{I, \theta\} = 1,$$

where $K$ is a smooth enough one variable increasing function such that $K(0) = 0$, $\{\cdot, \cdot\}$ stands for the standard Poisson bracket and is defined by

$$\{g, h\} := \partial_x g \partial_y h - \partial_y g \partial_x h$$

for our Hamiltonian system with one degree of freedom. In the action-angle coordinates the Hamiltonian is a function of the action coordinate alone. Writing $\varphi = (\varphi_1, \varphi_2)$ the action-angle coordinates are defined as

$$I(t) := I(X(t)) = \varphi_1(X(t)), \quad \text{and} \quad \theta(t) := \theta(X(t)) = \varphi_2(X(t)).$$

In the absence of random fluctuations (that is when $V \equiv 0$), the Hamiltonian system $\dot{X}(t) = \nabla^\perp H(X(t))$ becomes

$$\dot{I}(t) = 0 \quad \text{and} \quad \dot{\theta}(t) = \omega(I(t)), \quad \text{with} \quad \omega(I) = K'(I).$$

In the presence of random fluctuations, the Hamiltonian system (1.3) now becomes

$$\dot{I}(t) = \epsilon V(t) a(I(t), \theta(t)) \quad \text{and} \quad \dot{\theta}(t) = \omega(I(t)) + \epsilon V(t) b(I(t), \theta(t)), \quad (1.15)$$

where

$$a = e_2 \cdot \nabla \varphi_1 \circ \varphi^{-1} \quad \text{and} \quad b = e_2 \cdot \nabla \varphi_2 \circ \varphi^{-1}.$$

As an example, when the Hamiltonian is quadratic the action-angle coordinates read

$$x = \sqrt{\frac{I}{\pi}} \cos \left( 2\pi\theta - \frac{\pi}{2} \right) \quad \text{and} \quad y = \sqrt{\frac{I}{\pi}} \sin \left( 2\pi\theta - \frac{\pi}{2} \right).$$

One can then see that the action variable is a multiple of the Hamiltonian,

$$K(I) = \frac{I}{2\pi} \quad \text{and} \quad \omega(I) = \frac{1}{2\pi},$$

while the angle variable corresponds to the angle of a trajectory describing a circle with period 1 in this case. One can show that the functions $a$ and $b$ are in this case of the form

$$a(I, \theta) = 2\sqrt{\pi I} \sin \left( 2\pi\theta - \frac{\pi}{2} \right) \quad \text{and} \quad b(I, \theta) = \frac{1}{2\sqrt{\pi I}} \cos \left( 2\pi\theta - \frac{\pi}{2} \right).$$

In the general case, to study the long-time behavior of (1.15), we consider

$$I^{\epsilon}(t) = I(t/\epsilon^2) \quad \text{and} \quad \theta^{\epsilon}(t) = \theta(t/\epsilon^2),$$

satisfying

$$\dot{I}^{\epsilon}(t) = \frac{1}{\epsilon} V\left(\frac{t}{\epsilon^2}\right) a(I^{\epsilon}(t), \theta^{\epsilon}(t)) \quad \text{and} \quad \dot{\theta}^{\epsilon}(t) = \frac{\omega(I^{\epsilon}(t))}{\epsilon^2} + \frac{1}{\epsilon} V\left(\frac{t}{\epsilon^2}\right) b(I^{\epsilon}(t), \theta^{\epsilon}(t)), \quad (1.16)$$

with $I^\epsilon(0) = I_0$ and $\theta^\epsilon(0) = \theta_0$. As illustrated in Prop. 1.3.1 with the cos function in $\Gamma_c$, a key feature to handle the long-range correlations in the asymptotic analysis is the mean-zero property in $\theta$ for $a$ and $b$. In fact, if the system possesses a Fourier mode of order zero, that is with no oscillatory
behavior, a noise with long-range correlations will charge this component and cause the system to blow up as $\epsilon \to 0$. Due to the symplectic property of the transformation $\varphi$,

$$a(I, \theta) = -\partial_{\theta} \varphi_1^{-1}(I, \theta)$$

is clearly mean-zero with respect to $\theta$ as being the derivative of a periodic function, but regarding $b$ it is not as clear. Following the approach of [78], the fast and slow component of the angular variable are separated, and $\theta^\epsilon$ is split into two parts,

$$\theta^{\epsilon} = \psi^{\epsilon} + \tau^{\epsilon}.$$

Here, we set

$$\dot{t}^{\epsilon}(t) = \frac{\omega(I^{\epsilon}(t))}{\epsilon^2} + \frac{1}{\epsilon}V\left(\frac{t}{\epsilon^2}\right)\langle b(I^{\epsilon}(t), \cdot) \rangle \quad (1.17)$$

with initial condition $\tau^\epsilon(0) = 0$, and

$$\langle g \rangle := \int_{0}^{1} g(\theta) \, d\theta.$$

Regarding the slow motion variables $(I^\epsilon, \psi^\epsilon)$, the system (1.16) becomes

$$\dot{I}^{\epsilon}(t) = \frac{1}{\epsilon}v\left(\frac{t}{\epsilon^2}\right)A(I^{\epsilon}(t), \psi^{\epsilon}(t), \tau^{\epsilon}(t)) \quad \text{and} \quad \dot{\psi}^{\epsilon}(t) = \frac{1}{\epsilon}V\left(\frac{t}{\epsilon^2}\right)B(I^{\epsilon}(t), \psi^{\epsilon}(t), \tau^{\epsilon}(t)) \quad (1.18)$$

with

$$A(I, \psi, \tau) := a(I, \psi + \tau) \quad \text{and} \quad B(I, \psi, \tau) := b(I, \psi + \tau) - \langle b(I, \cdot) \rangle,$$

which are both mean-zero with respect to $\tau$. The above equations are coupled with the initial conditions $I^\epsilon(0) = I_0$, and $\psi^\epsilon(0) = \theta_0$.

Thm. 1.3.2 below obtains the limiting behavior of (1.18) as $\epsilon \to 0$ under the following assumptions on $K$:

* The function $K$ is smooth, and

$$\inf_{I \ge 0} K'(I) = \inf_{I \ge 0} \omega(I) > \omega_0 > 0, \quad (1.19)$$

for some strictly positive number $\omega_0$;

* There exists $r > 0$, and positive constants $c_{1,r}, c_{2,r} > 0$, such that for any $I \in (0, r)$

$$c_{1,r}I \le K(I) \le c_{2,r}I \quad \text{and} \quad |\omega'(I)| \le \frac{c_{2,r}}{I}. \quad (1.20)$$

These conditions imply that $K(0) = 0$ and that $K$ is an increasing function in $I$. These requirements are not too restrictive as any Hamiltonian satisfies (1.19) and (1.20) near non-degenerate critical points.

**Theorem 1.3.2** *Assuming (1.19) and (1.20), the family $(I^\epsilon, \psi^\epsilon)_\epsilon$, defined by (1.17)-(1.18), converges in law in $C([0, \infty), \mathbb{R}^2)$, as $\epsilon \to 0$, to a diffusion process $(I(t), \psi(t))_{t \ge 0}$ where $(I_t)_{t \ge 0}$ is the unique weak solution to the SDE*

$$dI(t) = \int_{0}^{1} a(I(t), \tau) dW(t, I(t), \tau) d\tau \\ + \left[ \int_{0}^{\infty} R(s) \int_{0}^{1} \partial_{I} (a(I, \tau + \omega(I)s))_{|I=I(t)} a(I(t), \tau) \right. \\ \left. + a(I(t), \tau + \omega(I(t))s) \partial_{I} a(I(t), \tau) d\tau ds \right] dt,$$

with initial condition $I(0) = I_0$. Here, $W$ is a real valued Brownian field with covariance function

$$\mathbb{E}[W(t_1, y, \phi_1)W(t_2, y, \phi_2)] = t_1 \wedge t_2 \int_0^\infty R(s) \\ \times \int_0^1 \phi_1(\tau + \omega(y)s)\phi_2(\tau) + \phi_1(\tau)\phi_2(\tau + \omega(y)s) \, d\tau \, ds,$$
for any $\phi_1, \phi_2 \in L_0^2(\mathbb{T})$ where

$$L_0^2(\mathbb{T}) := \{\phi \in L^2(\mathbb{T}) : \langle \phi \rangle = 0\}.$$

Also, we have

$$d\psi(t) = \int_0^1 b(I(t), \tau) dW(t, I(t), \tau) d\tau + \left[ \int_0^\infty R(s) \int_0^1 \partial_I (b(I, \tau + \omega(I)s))_{|I=I(t)} a(I(t), \tau) + b(I(t), \tau + \omega(I(t))s) \partial_I a(I(t), \tau) d\tau ds \right] dt,$$

with $\psi(0) = \theta_0$.

Note that the action variable $I$ does not depend on the slow angular motion $\psi$. However, the distribution of $\psi$ is completely determined by the motion of the action variable $I$ and the Brownian field.

The proof of Thm. 1.3.1 follows from a direct application of the Itô formula applied to the limiting diffusion process $I$.

For the quadratic Hamiltonian, stochastic differential equations for $(I(t))_{t \ge 0}$ and $(\psi(t))_{t \ge 0}$ can be directly derived leading to

$$dI_t = 2\sqrt{mI_t} dB_t^1 + 2m dt \quad \text{and} \quad d\psi_t = \frac{m}{2\pi\sqrt{I_t}} dB_t^2,$$

with

$$m = \pi \int_0^\infty R(u) \cos(u) du,$$

and where $B^1$ and $B^2$ are two independent standard Brownian motions. One can easily remark that $(I_t/m)_{t \ge 0}$ is a 2D squared Bessel process, so that $(I_t/2\pi)_{t \ge 0}$ has the same law as $(|W_t|^2/2)_{t \ge 0}$, for $W$ being a 2D Brownian motion with covariance matrix (1.11). Then, we recover the result of Prop. 1.3.1, but we also access to the effective behavior of the slow component of the angle variable.

## 1.4 Perspectives

The present analysis, which leads to diffusion behaviors, relies on the perturbed-test-function technique together with a martingale approach, limiting the class of noises that can be considered. It can be expected to extend the class of noise by the use of the more general, but more involved, rough path framework [74,144]. Further investigations are required in this direction in order to overcome the technical difficulty of the dimensionality as well as dealing with nonlinear relations involving the noise. In the present context, these mentioned nonlinear relations can be seen in (1.18) together with (1.17) where the blowing term $V(t/\epsilon^2)/\epsilon$ is involved. The exact same difficulty is observed for the wave propagation problem discussed in Sect. 3.3 of Chap. 3. At present, we are only able to analyze the scenario where the Hamiltonian is smooth with exactly one non-degenerate critical point. The Reeb graph of such Hamiltonian is simply the half line, and the limiting process in Thm. 1.3.1 is a diffusion on the half line. For general Hamiltonians, Reeb graph may have multiple vertices, and the limiting process should be a diffusion on this Reeb graph. When the driving noise is a standard Brownian motion (as in the Freidlin-Wentzell theory), this diffusion is characterized through its generator on each edge of the Reeb graph, and certain gluing conditions on the vertices. While the proofs in the Freidlin-Wentzell theory rely heavily on the Markov property, and do not apply to our situation, the form of the infinitesimal generator (1.14) is similar to that in [71, Eq. (1.8)]. As a result, similar gluing conditions are expected in our context. Explicitly, for a sufficiently general class of Hamiltonians, we expect that the process $(H(X^\epsilon))_\epsilon$ converges in law to a diffusion process on the Reeb graph. The generator of this diffusion should be

$$L_j h_j(\mathcal{X}) = \frac{1}{2\Lambda_j(\mathcal{X})} \partial_{\mathcal{X}}(\Sigma_j(\mathcal{X}) \partial_{\mathcal{X}} h_j)$$
on each edge $I_j$. Considering $C_j(\mathcal{X})$ the connected components of the level set $\{H = \mathcal{X}\}$ corresponding to a point $\mathcal{X} \in I_j$, $\Sigma_j$ and $\Lambda_j$ should be

$$ \Sigma_j(\mathcal{X}) = 2 \int_{u=0}^{\infty} du R(u) \oint_{C_j(\mathcal{X})} \partial_y H(\tilde{X}_u^{x,y}) \partial_y H(x, y) \frac{dl(x, y)}{|\nabla H(x, y)|}, $$

and

$$ \Lambda_j(\mathcal{X}) = \oint_{C_j(\mathcal{X})} \frac{dl(x, y)}{|\nabla H(x, y)|}. $$

To describe the gluing conditions on the vertices, consider an interior vertex $O_k$ corresponding to the saddle point $(\mathbf{x}_k, \mathbf{y}_k)$ and the level set $\{H = \mathcal{X}_k\}$, at which the collection of edges $\{I_j\}$ meet. Now define

$$ \beta_{kj} = \lim_{\mathcal{X} \to \mathcal{X}_k} C_j(\mathcal{X}), $$

where the limit is taken along the edge $I_j$. The gluing condition at the vertex $O_k$ can then be written as

$$ \sum_{j: I_j \sim O_k} \pm \beta_{kj} h_j'(H(\mathbf{x}_k, \mathbf{y}_k)) = 0, $$

where the sign before $\beta_{kj}$ is positive if $\mathcal{X} > \mathcal{X}_k$ along the edge $I_j$, and negative if $\mathcal{X} < \mathcal{X}_k$ along the edge $I_j$.
Chapter 2

Random Schrödinger equation with
long-range correlations

The Schrödinger equation with a time-dependent random potential has been studied for a long time [17, 64–66, 115, 193], with a particular interest for applications in wave propagation in random media under the so-called paraxial (or parabolic) approximation [11, 18, 34, 52, 80, 87, 151, 196, 201]. The paraxial approximation allows a simplified treatment of wave propagation phenomena at both the theoretical and numerical point of views. The main idea is to provide an approximation of a solution u to the scalar wave equation

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

with

$$
\frac{1}{c^2(z, \mathbf{x})} = \frac{1}{c_0^2} (1 + V(z, \mathbf{x})),
$$

where $\Delta = \partial_{zz}^2 + \Delta_x$ is the Laplacian operator, and $V$ is a mean-zero random field describing fluctuations of the wave-speed profile around its background value $c_0$.

Considering the z-axis as the main propagation axis, that is along which the wave is mainly propagating, the paraxial approximation consists in approximating $u$ by

$$
u(t, z, \mathbf{x}) \approx \int \phi_{\omega}(z, \mathbf{x}) e^{i\omega(z/c_0 - t)} d\omega.
$$

From this representation, the wave is assumed to propagate in only one direction along the z-axis, neglecting the backscattering phenomena. This latter aspect is the key point that simplifies theoretical and numerical treatments of wave propagation in random media. The envelop $\phi_{\omega}$ therefore satisfies the following random Schrödinger equation

$$
i k_{\omega} \partial_z \phi_{\omega} + \frac{1}{2} \Delta_{\mathbf{x}} \phi_{\omega} + k_{\omega}^2 V(z, \mathbf{x}) \phi_{\omega} = 0, \quad (2.1)
$$

with wavenumber $k_\omega = \omega/c_0$. In this Schrödinger equation, the time variable is played by the
$z$-variable of the main propagation axis on which the random fluctuation $V$ depends. A rigorous
treatment of this approximation will be discussed in Sect. 3.2 of Chap. 3.

The random Schrödinger equation involving a time-dependent potential,

$$
i \partial_t \phi + \frac{1}{2} \Delta_x \phi - \sqrt{\varepsilon} V(t, \mathbf{x}) \phi = 0 \quad (t, \mathbf{x}) \in (0, \infty) \times \mathbb{R}^d, \quad (2.2)
$$

have been considered in [14, 17] to derive a radiative transfer model under short-range correlations, but also in [15] under both short-range and long-range correlations. In this latter work the asymptotic behavior of the wave function itself is explored. Under long-range correlations the authors prove a non-central limit theorem type result (described below) using a moment technique. This chapter
presents other scaling regimes exhibiting non-trivial effects under long-range correlations. In contrast with Chap. 1, where no spatial effect can take place, the multiscale properties of the potential w.r.t. the variable $\mathbf{x}$ enable intermediate non-trivial regimes from the non-central type result up to the radiative transfer regime under the so-called semi-classical limit. This chapter is based on the results obtained in [3, 4, 6] of the publication list on page 6.

To highlight the effective phenomena produced by long-range correlations over a range of scaling regimes, we introduce the time and space scaling:

$$t \to t/\epsilon^s \quad \text{and} \quad \mathbf{x} \to \mathbf{x}/\epsilon^s, \quad s \in (0, 1],$$

where $s$ represents the propagation scale parameter. The wave function associated to this scaling is

$$\phi_{\epsilon}(t, \mathbf{x}) := \phi\left(\frac{t}{\epsilon^s}, \frac{\mathbf{x}}{\epsilon^s}\right),$$

with corresponding Schrödinger equation

$$i\epsilon^s \partial_t \phi_{\epsilon} + \frac{\epsilon^{2s}}{2} \Delta_{\mathbf{x}} \phi_{\epsilon} - \sqrt{\epsilon} V \left( \frac{t}{\epsilon^s}, \frac{\mathbf{x}}{\epsilon^s} \right) \phi_{\epsilon} = 0, \quad \text{and} \quad \phi_{\epsilon}(0, \mathbf{x}) = \phi_0(\mathbf{x}/\epsilon^{s-s_c}). \quad (2.3)$$

The parameter $s_c$ accounts for low frequency oscillations of the wave function. Such low frequencies are the key points when describing how the loss of coherence affects the wave function over the different propagation scales. Note that the standard semi-classical limit holds for $s = 1$ and $s_c = 0$. An anisotropy parameter $\varrho \in [0, 1)$ can also be added into the fluctuation model,

$$\epsilon^{\frac{1-\varrho}{2}} V \left( \frac{t}{\epsilon^{s+\varrho}}, \frac{\mathbf{x}}{\epsilon^s} \right).$$

Recalling that for the paraxial approximation the time variable $t$ plays the role of the main propagation axis (see (2.1)), $\varrho > 0$ corresponds to faster fluctuations along this axis than in the transverse direction. In other words, this condition corresponds to a statistical anisotropy of the medium fluctuations.

The organization of this chapter is as follows. After a brief description of the random potential under consideration in Sect. 2.1, the different non-trivial regimes describing the loss of coherence are presented. Sect. 2.2 reviews the result obtained in [15] regarding the onset of the random phase modulation. In Sect. 2.3 a scaled Wigner transform is introduced. Its scaling enables the description of the loss of coherence, in Sect. 2.4 and 2.5, over different scaling regimes and how this loss of coherence propagates from the large scales supported by the wave function to the finer ones. Finally, Sect. 2.6 introduces a time-splitting scheme that are able to capture statistically all the scaling regimes we describe under long-range correlations.

## 2.1 The random potential

Throughout this chapter the random potential is defined following the lines of Chap. 1 (see also [15]) by

$$V(t, \mathbf{x}) := \int_{-\infty}^{t} \int_{S} e^{i\mathbf{p} \cdot \mathbf{x}} e^{-\mu|\mathbf{p}|^{2\beta}(t-u)} \mathcal{B}(du, d\mathbf{p}),$$

where $S \subseteq \mathbb{R}^d$ is a bounded symmetric open set containing 0. Here, $\mu$ and $\beta$ are constants, and $\mathcal{B}$ is a Gaussian random measure which is white in time and colored in space:

$$\mathcal{B}(du, d\mathbf{p}) := \sqrt{2\mu r(\mathbf{p})} |\mathbf{p}|^{\beta} \mathbf{1}_{S}(\mathbf{p}) \xi(du, d\mathbf{p}),$$

where $\xi$ is a $\mathbb{R}^{d+1}$ Gaussian white noise, and $r: S \setminus \{0\} \to [0, \infty)$ is defined again by

$$r(\mathbf{p}) := \frac{\lambda(\mathbf{p})}{|\mathbf{p}|^{d-1+2\alpha}}. \quad (2.4)$$

Here, $\lambda: S \to \mathbb{R}_+$ is a smooth bounded even function such that $\lambda(0) \neq 0$. The correlation function of $V$ is this time given by

$$R(t, \mathbf{x}) := \mathbb{E}[V(t + s, \mathbf{x} + \mathbf{y})V(s, \mathbf{y})] = \int_S e^{i\mathbf{p} \cdot \mathbf{x}} r(\mathbf{p}) e^{-\mu|\mathbf{p}|^{2\beta}|t|} d\mathbf{p}.$$
Again, the condition

$$ \alpha < \frac{1}{2} $$

is required for the process and correlation function to be well defined. In this chapter the condition $\beta \in (0, 1/2]$ is also required for technical reasons.

One can observe from (2.4) that the random field $V(t)$, for any fixed time $t$, exhibits long-range correlations in space. In fact, the function $r$ is singular at $\mathbf{p} = 0$ (therefore unbounded), so that its inverse Fourier transform is not integrable. However, for any fixed $\mathbf{x} \in \mathbb{R}^d$,

$$ \lim_{t \to \infty} t^{\gamma} R(t, \mathbf{x}) = R_0, $$

where

$$ \gamma = \frac{1 - 2\alpha}{2\beta} \quad \text{and} \quad R_0 = \lambda(0) \int \frac{e^{-\mu|\mathbf{p}|^{2\beta}}}{|\mathbf{p}|^{d-1+2\alpha}} d\mathbf{p}. $$

Long-range correlations (in time) are therefore considered here by assuming

$$ 2(\alpha + \beta) > 1, $$

but short-range correlations (still in time) can also be obtained by assuming $2(\alpha + \beta) < 1$.

Note that Eq. (2.2) is a time evolution problem, so we have to take care about the statistical properties of the random potential $V$ w.r.t. the temporal variable. It is worth mentioning that if $V$ is assumed to have rapidly decaying correlations in time, $(V(t_1), V(t_2))$ exhibits now rapidly decaying cross-spatial correlations, and the evolution problem (2.2) behaves like in the short-range case addressed in [15]. As a result, even if at each fixed time the spatial correlations are slowly decaying, the resulting time evolution problem behaves as if the random potential has rapidly decaying correlations in all variables.

As in Sect. 1.3.2 of Chap. 1, this specific form for $V$ is considered to carry out the asymptotic analysis.

## 2.2 Phase modulation

As in Chap. 1, a nontrivial scaling regime, corresponding to the non-central limit theorem, can be exhibited. In the context of this chapter, it is given by the following result (provided by [15,93]) that can be proved through a moment technique.

**Theorem 2.2.1** *Denoting*

$$ \mathfrak{H} := 1 - \gamma/2 \in (1/2, 1) \quad \text{and} \quad \mathfrak{H}_{\varrho} := \frac{\mathfrak{H}}{1 - \varrho\left(\frac{\alpha + \beta - 1/2}{\beta}\right)}, \quad (2.5) $$

for $\varrho \in [0, 1)$, the process $\hat{\xi}_{\mathfrak{H}, \varrho}(t, \mathbf{k})$ defined by

$$ \hat{\xi}_{\mathfrak{H}, \varrho}(t, \mathbf{k}) := \frac{1}{\varepsilon^{d_s}} \hat{\phi}_{\varrho}(t, \frac{\mathbf{k}}{\varepsilon^{s-s_c}}) e^{i|\mathbf{k}|^2 t / (2\varepsilon^{s-2s_c})}, \quad \text{with } s = 1/(2\mathfrak{H}_{\varrho}), \quad (2.6) $$

converges in law to

$$ \hat{\xi}(t, \mathbf{k}) = \hat{\phi}_0(\mathbf{k}) \exp \left( i \sqrt{D(\alpha, \beta, \mathbf{k})} W_{\mathfrak{H}}(t) \right), $$

for each $t \ge 0$ and $\mathbf{k} \in \mathbb{R}^d$. Here, $W_{\mathfrak{H}}$ is a standard fractional Brownian motion with Hurst index $\mathfrak{H}$,

$$ D(\alpha, \beta, \mathbf{k}) := \frac{a(0)}{(2\pi)^d \mathfrak{H}(2\mathfrak{H} - 1)} \int_0^\infty d\rho \frac{e^{-\mu\rho}}{\rho^{2\alpha-1}} \int_{\mathbb{S}^{d-1}} dS(u) e^{i|\mathbf{k}|\rho u \cdot \mathbf{e}_1}, \quad (2.7) $$

if $\beta = 1/2$, $\varrho = 0$, $s_c = 0$, and

$$ D(\alpha, \beta, \mathbf{k}) := D(\alpha, \beta) = \frac{a(0)\Omega_d}{(2\pi)^d \mathfrak{H}(2\mathfrak{H} - 1)} \int_0^\infty d\rho \frac{e^{-\mu\rho^{2\beta}}}{\rho^{2\alpha-1}}, \quad (2.8) $$

otherwise, where $\Omega_d$ stands for the surface area of the unit sphere in $\mathbb{R}^d$, and $\mathbf{e}_1 \in \mathbb{S}^{d-1}$.
In this result, the term $\hat{\xi}_{\hat{\mathfrak{g}}_\varrho, \varepsilon}$ is scaled according to the initial condition $\phi_\varepsilon(t=0, \mathbf{x}) = \phi_0(\mathbf{x}/\varepsilon^{s-s_c})$. The resulting highly-oscillating behavior of the solution $\phi_\varepsilon$ is filtered out through the term $e^{i|\mathbf{k}|^2t/(2\varepsilon^{s-2s_c})}$. As a result, $\hat{\xi}_{\hat{\mathfrak{g}}_\varrho, \varepsilon}$ exhibits mainly an effective random behavior through the term $\exp(i\sqrt{D(\alpha, \beta, \mathbf{k})}W_{\hat{\mathfrak{g}}}(t))$, corresponding to a random phase modulation.

As described in Chap. 1, when the system propagates over longer time scales, other effective phenomena can be observed. Over longer propagation scales $s > 1/(2\hat{\mathfrak{g}}_\varrho)$ such phase modulations oscillate fast enough to affect the coherence of the wave function. To understand how these phenomena take place a properly scaled Wigner transform and initial conditions are considered.

## 2.3 The Wigner transform

To highlight the loss of coherence effects on the wave function $\phi_\varepsilon$, we can consider an initial conditions $\phi_{0,\varepsilon}$, for (2.3), of the form:

$$ \phi_{0,\varepsilon}(\mathbf{x}) = \phi_0(\mathbf{x}) \exp(i\zeta \cdot \mathbf{x}/\varepsilon^{s-s_c}). \quad (2.9) $$

This initial condition represents a quasi-plane wave with direction $\zeta$, and slowly varying envelop $\phi_0$. The plane wave is slowly varying for $s_c = s$, and highly oscillating when $s_c < s$. In this latter case self-averaging effects take place and the loss of coherence mechanism will be deterministic. It is worth mentioning that the spatial frequency of the initial condition ($\sim \varepsilon^{-(s-s_c)}$) is low compared to the one of the random fluctuations ($\sim \varepsilon^{-s}$). For random media with rapidly decaying correlations, such low spatial frequency sources do not interact with the random medium, and provide trivial limits. However, they do interact with slowly decorrelating random media and provide nontrivial limits. Due to technical difficulties in the asymptotic analysis, it is convenient to consider the direction $\zeta$ as being randomly distributed, according to some probability law $\nu(d\zeta)$ over some space $U \subset \mathbb{R}^d$, in order to provide a rigorous derivation of the self-averaging effects.

The loss of coherence phenomena are analyzed through the following scaled Wigner transform

$$ \begin{aligned} W_\varepsilon(t, \mathbf{x}, \mathbf{k}) &:= \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d \times U} dy \nu(d\zeta) e^{i\mathbf{k} \cdot \mathbf{y}} \phi_\varepsilon \left( t, \mathbf{x} - \varepsilon^{s-s_c} \frac{\mathbf{y}}{2}, \zeta \right) \overline{\phi_\varepsilon \left( t, \mathbf{x} + \varepsilon^{s-s_c} \frac{\mathbf{y}}{2}, \zeta \right)} \\ &= \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d \times U} dy \nu(d\zeta) e^{i\mathbf{k} \cdot \mathbf{y}} \phi \left( \frac{t}{\varepsilon^s}, \frac{\mathbf{x}}{\varepsilon^s} - \frac{\mathbf{y}}{2\varepsilon^{s_c}}, \zeta \right) \overline{\phi \left( \frac{t}{\varepsilon^s}, \frac{\mathbf{x}}{\varepsilon^s} + \frac{\mathbf{y}}{2\varepsilon^{s_c}}, \zeta \right)}, \end{aligned} $$

The terminology loss of coherence refers to diffusion phenomena through the momentum variable $\mathbf{k}$. The Wigner transform can be seen as the Fourier transform of the two-point correlation function of the wave function, and where $s_c \in [0, s]$ corresponds to a spatial correlation parameter. This parameter captures the evolution of the degree of correlation of the wave function. We refer to Fig. 2.1 for some illustrations. Let us note that the cases $s_c < s$ study the local loss of coherence, while the case $s = s_c$ study the nonlocal loss of coherence. As we will see, for a given propagation scale parameter $s$ the loss of coherence of the wave field $\phi_\varepsilon$ can be observed at a particular spatial correlation parameter $s_c$. The larger the propagation scale parameter is the shorter the correlation scale parameter $s_c$ is. In other words, depending on the propagation time and distance we adjust the spatial correlation scale to exhibit the loss of coherence. It appears that as the wave function propagates over longer scales, the correlation scales affected by the random fluctuations become smaller.

The developments of this chapter are based on the following equation satisfied by the Wigner transform,

$$ \begin{aligned} \partial_t W_\varepsilon(t, \mathbf{x}, \mathbf{k}) + \varepsilon^{s_c} \mathbf{k} \cdot \nabla_\mathbf{x} W_\varepsilon(t, \mathbf{x}, \mathbf{k}) = & \\ \frac{\varepsilon^{(1-\varrho)/2-s}}{(2\pi)^d} \int_{\mathbb{R}^d} \hat{V} \left( \frac{t}{\varepsilon^{s+\varrho}}, d\mathbf{p} \right) e^{i\mathbf{p} \cdot \mathbf{x}/\varepsilon^s} \left( W_\varepsilon \left( t, \mathbf{x}, \mathbf{k} - \frac{\mathbf{p}}{2\varepsilon^{s_c}} \right) - W_\varepsilon \left( t, \mathbf{x}, \mathbf{k} + \frac{\mathbf{p}}{2\varepsilon^{s_c}} \right) \right). & \end{aligned} \quad (2.10) $$

If $V \equiv 0$, the Wigner transform is simply given by $W_\varepsilon(t, \mathbf{x}, \mathbf{k}) = W_\varepsilon(t=0, \mathbf{x}-\varepsilon^{s_c}t\mathbf{k}, \mathbf{k})$, the momentum of $W_\varepsilon$ is preserved during the propagation, that is there is no loss of coherence. For nonnull random fluctuations, the radiative transfer equation (2.10) describes the loss of coherence of the function $\phi_\varepsilon$ through the random integral operator. However, depending on the scaling under consideration, this operator can be of negligible effect, and no loss of coherence can be observed. For example, when considering the correlation scale $s_c = 0$, no loss of coherence can be observed before $s = 1$ (semi-classical limit). Also, the dispersion term $\mathbf{k} \cdot \nabla_\mathbf{x}$ is of order $\varepsilon^{s_c}$, so that dispersion is only capture for $s_c = 0$.
<figure><img src="image_4.png" /><figcaption>Fig. 2.1: Illustration of the initial condition (2.9). (a) and (b) represent low spatial frequency initial conditions. (a) represents the case $s_c = s$, and (b) represents the case $s_c < s$. (c) represents the case $s_c = 0$ and $s = 1$ corresponding to the radiative transfer regime.</figcaption></figure>
## 2.4 Loss of coherence

As described in Thm. 2.4.2, Thm. 2.4.1, and Thm. 2.5.1 below, at a given propagation scale $s$, the corresponding spatial correlation parameter is given by

$$s_c = \frac{1-s}{\theta}, \qquad (2.11)$$

with

$$\theta = 2(\alpha + \beta) - 1 \in (0, 1), \qquad (2.12)$$

so that for $s_c \le s$ we have

$$\frac{1}{2\mathfrak{H}_{\varrho}} \le \frac{1}{1+\theta} \le s.$$

As a result, no loss of coherence can be observed before the propagation scale $s = 1/(1+\theta)$, which is larger than the propagation scale to observe the onset of the phase modulation. This relation holds true for any spatial correlation parameter $s_c \le s$. Furthermore, Eq. (2.11) expresses how the spatial correlation scale decreases as the propagation scale increases.

The case $s = s_c = 1/(1+\theta)$ corresponds to the first effects happening after the onset of the random phase modulation, and they hold on the envelop $\phi_0$ itself and not at a local level. These effects can be described as follows.

**Theorem 2.4.1** For either $\varrho > 0$ or $\beta < 1/2$, and

$$s_c = s = \frac{1}{1+\theta},$$

the family $(W_{\varepsilon})_{\varepsilon}$ converges in law in $C([0, \infty), L^2(\mathbb{R}^{2d}))$, as $\varepsilon \to 0$, to a limit $W$ defined by

$$W(t, \mathbf{x}, \mathbf{k}) = \frac{1}{(2\pi)^d} \int d\mathbf{q} \hat{W}_0^{\mathbf{k}}(\mathbf{x}, \mathbf{q}) \exp \left( i\mathbf{k} \cdot \mathbf{q} + i \int \mathcal{B}(t, d\mathbf{p}) e^{i\mathbf{p} \cdot \mathbf{x}} \left( e^{-i\mathbf{q} \cdot \mathbf{p}/2} - e^{i\mathbf{q} \cdot \mathbf{p}/2} \right) \right), \quad (2.13)$$

where $\hat{W}_0^{\mathbf{k}}$ stands for the Fourier transform of

$$W_0(\mathbf{x}, \mathbf{k}) := W(t=0, \mathbf{x}, \mathbf{k}) = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^{2d}} d\mathbf{y} \nu(d\zeta) e^{i(\mathbf{k}-\zeta)\cdot\mathbf{y}} \phi_0(\mathbf{x}-\mathbf{y}/2) \overline{\phi_0(\mathbf{x}+\mathbf{y}/2)}$$

with respect to the $\mathbf{k}$-variable. Here, $(\mathcal{B}_t)_t$ is a Wiener process on the dual of the Hilbert space

$$\mathcal{H}_{\theta} := \left\{ \varphi : \varphi(\mathbf{p}) = \overline{\varphi(-\mathbf{p})} \quad \text{and} \quad \int \frac{d\mathbf{p}}{|\mathbf{p}|^{d+\theta}} |\varphi(\mathbf{p})|^2 < \infty \right\},$$

such that $\mathcal{B}^*(t, d\mathbf{p}) = \mathcal{B}(t, -d\mathbf{p})$, and with covariance function

$$\mathbb{E}[\mathcal{B}_t(\varphi)\mathcal{B}_s(\psi)] = t \wedge s \int \frac{d\mathbf{p}}{|\mathbf{p}|^{d+\theta}} \varphi(\mathbf{p}) \overline{\psi(\mathbf{p})} \quad \varphi, \psi \in \mathcal{H}_{\theta}.$$

Moreover, $W$ is the unique weak solution of the stochastic partial differential equation

$$dW(t, \mathbf{x}, \mathbf{k}) = -\sigma(\theta)(-\Delta_{\mathbf{k}})^{\theta/2} W(t, \mathbf{x}, \mathbf{k}) \\ + \frac{2ia(0)}{(2\pi)^d} \int d\mathcal{B}(t, d\mathbf{p}) e^{i\mathbf{x} \cdot \mathbf{p}} \left( W(t, \mathbf{x}, \mathbf{k} - \frac{\mathbf{p}}{2}) - W(t, \mathbf{x}, \mathbf{k} + \frac{\mathbf{p}}{2}) \right), \quad (2.14)$$

where

$$\sigma(\theta) = \frac{2a(0)\theta\Gamma(1-\theta)}{(2\pi)^d} \int_{\mathbb{S}^{d-1}} dS(\mathbf{u}) |\mathbf{e}_1 \cdot \mathbf{u}|^{\theta}, \quad (2.15)$$

with $\mathbf{e}_1 \in \mathbb{S}^{d-1}$, and $\Gamma$ the standard Gamma function.
The representation (2.13) is still random as it exhibits, in the Fourier domain, a random phase modulation. In this case the wave does not propagate enough to produce self-averaging phenomena as for the standard semi-classical limit $s = 1$ [11, 17, 18]. The case $\rho = 0$, $\beta = 1/2$, and $s = s_c$ is addressed by Thm. 2.2.1, since in this particular case $1/(1 + \theta) = 1/(2\zeta)$, and no loss of coherence takes place.

As described in the following result, for a longer propagation scale $s > 1/(1 + \theta)$ the random phase modulation in the Wigner transform starts to oscillate fast enough to average out and produces a homogenization phenomenon. It turns out that the deterministic limit corresponds to the expectation of the limiting Wigner transform obtained for $s_c = s = 1/(1 + \theta)$. Both Thm. 2.4.1 and the following result are proved through approximation-diffusion theorems based on a perturbed-test-function approach [139].

**Theorem 2.4.2** For any $s \in (1/(1 + \theta), 1)$, and

$$s_c = \frac{1-s}{\theta} < s, \qquad (2.16)$$

the family $(W_\varepsilon)_\varepsilon$ converges in probability in $C([0, +\infty), L_w^2(\mathbb{R}^{2d}))$, as $\varepsilon \to 0$, to the unique solution of the fractional diffusion equation

$$\partial_t W = -\sigma(\theta)(-\Delta_k)^{\theta/2} W, \qquad (2.17)$$

with

$$W_0(\mathbf{x}, \mathbf{k}) := W(t=0, \mathbf{x}, \mathbf{k}) = |\phi_0(\mathbf{x})|^2 \nu(\mathbf{k}), \qquad (2.18)$$

and $\sigma(\theta)$ defined by (2.15). Here, $L_w^2(\mathbb{R}^{2d})$ stands for $L^2(\mathbb{R}^{2d})$ equipped with the weak topology, $(-\Delta_k)^{\theta/2}$ is the fractional Laplacian with Hurst index $\theta \in (0, 1)$. Moreover, $W$ is given by

$$W(t, \mathbf{x}, \mathbf{k}) = \frac{1}{(2\pi)^d} \int d\mathbf{q} \exp(i\mathbf{k} \cdot \mathbf{q} - \sigma(\theta)|\mathbf{q}|^\theta t) \hat{W}_0^{\mathbf{k}}(\mathbf{x}, \mathbf{q}),$$

where $\hat{W}_0^{\mathbf{k}}$ stands for the Fourier transform of $W_0$ with respect to the $\mathbf{k}$-variable.

The observed fractional diffusion exhibits a damping term that obeys to a power law with exponent $\theta \in (0, 1)$. This term characterizes how fast the wave function $\phi_\varepsilon$ loses its coherence.

We conclude this section by providing a heuristic explanation for why such loss of coherence phenomena are observable under long-range correlations but not under short-range correlations. After homogenization of the stochastic dynamic, the diffusion operator w.r.t. the momentum variable is approximately given by

$$\varepsilon^{1-s} \int d\mathbf{p} \sigma(\mathbf{p}) \left( \varphi\left(\mathbf{k} + \frac{\mathbf{p}}{\varepsilon^{s_c}}\right) - \varphi(\mathbf{k}) \right).$$

In case of rapidly decaying correlations $\sigma(\mathbf{p}) \in L^1(\mathbb{R})$, unless $s = 1$, the term $\varepsilon^{1-s}$ cannot be compensated, and we cannot observe any significant interaction between the wave function and the random fluctuations. However, random fluctuations with long-range correlations allow the following asymptotic

$$\varepsilon^{1-s} \int d\mathbf{p} \frac{a(\mathbf{p})}{|\mathbf{p}|^{d+\theta}} \left( \varphi\left(\mathbf{k} + \frac{\mathbf{p}}{\varepsilon^{s_c}}\right) - \varphi(\mathbf{k}) \right) \sim \varepsilon^{1-s-\theta s_c} a(0) \int \frac{d\mathbf{p}}{|\mathbf{p}|^{d+\theta}} \left( \varphi(\mathbf{k} + \mathbf{p}) - \varphi(\mathbf{k}) \right),$$

and therefore significant momentum diffusion as soon as $s_c = (1-s)/\theta$. In other words, the asymptotic autosimilarity property of the correlation function enables significant interactions between the low-frequency structures of the wave function, characterized by $s_c$, with the random fluctuations.

## 2.5 Radiative transfer regime

The scaling $s = 1$ and $s_c = 0$ corresponds to the standard semi-classical limit for which the loss of coherence is described by a radiative transfer equation (RTE). This limiting RTE holds for both slowly- and rapidly- decaying correlations.
**Theorem 2.5.1** The family $(W_\epsilon)_\epsilon$ converges in probability in $C([0, +\infty), L^2_w(\mathbb{R}^{2d}))$, as $\epsilon \to 0$, to a limit denoted by $W$, which is the unique, uniformly bounded in $L^2(\mathbb{R}^{2d})$, classical solution of the RTE

$$ \partial_t W + \mathbf{k} \cdot \nabla_x W = \mathcal{L}W, \qquad (2.19) $$

with $W(t=0) = W_0 \in L^2(\mathbb{R}^{2d})$ defined by (2.18). Here, $\mathcal{L}$ is defined by

$$ \mathcal{L}\varphi(\mathbf{k}) := \int d\mathbf{p} \sigma(\mathbf{p} - \mathbf{k}, \frac{|\mathbf{k}|^2 - |\mathbf{p}|^2}{2}) (\varphi(\mathbf{p}) - \varphi(\mathbf{k})), $$

for $\varphi \in C^\infty(\mathbb{R}^d)$, where

$$ \sigma(\mathbf{p}, \omega) = \sigma(\mathbf{p}) := \frac{2a(\mathbf{p})}{(2\pi)^d |\mathbf{p}|^{d+\theta}} \quad \text{for} \quad \varrho > 0, \qquad (2.20) $$

and

$$ \sigma(\mathbf{p}, \omega) := \frac{2\mu a(\mathbf{p})}{(2\pi)^d |\mathbf{p}|^{d+2\alpha-2\beta-1} (\mu^2 |\mathbf{p}|^{4\beta} + \omega^2)} \quad \text{for} \quad \varrho = 0. \qquad (2.21) $$

Moreover, $W \in C^\infty((0, \infty) \times \mathbb{R}^d \times \mathbb{R}^d)$ for any initial condition $W_0 \in L^2(\mathbb{R}^{2d})$.

Originally, the case $\varrho = 0$ was not covered in [92]. In this case the difficulty comes from the identification of the accumulation points once the tightness of the family $(W_\epsilon)_\epsilon$ is obtained. Due to the singular kernel in $\mathcal{L}$, this can be achieved by proving the weak uniqueness of the limiting radiative transfer equation. To obtain this weak uniqueness, it is enough to prove the smoothness of classical solutions for this equation. While the regularity can be obtained with relatively straightforward algebra for $\varrho > 0$, the case $\varrho = 0$ requires more efforts as described in Sect. 4.1 of Chap. 4 [97,98].

Interestingly, the radiative transfer model remains well-defined regardless of whether the integral $\int d\mathbf{p} \sigma(\mathbf{p}, \omega)$ is finite or not. In other words, the RTE (2.19) is valid under both slowly and rapidly decaying correlations in time. However, although long-range correlations lead to regularizing effects, short-range correlations do not yield similar properties. This regularizing properties highlight an important qualitative distinction between short-range and long-range correlations.

For long-range correlations, although the behavior of the momentum in (2.19) does not exactly align with the diffusion mechanism obtained in Thm. 2.4.2, they are very closed to each other. The momentum diffusion provided in Thm. 2.4.2 is characterized by a fractional Laplacian, whereas in the radiative transfer regime the momentum diffusion is described using a nonlocal operator that is almost a fractional Laplacian. Nevertheless, both of these diffusion mechanisms fall under the category of anomalous diffusion, as they result in damping terms obeying to a power law with exponent $\theta \in (0, 1)$. From the radiative transfer model, one needs to consider a proper scaling regime to once again observe a momentum diffusion given by a fractional Laplacian [92].

## 2.6 Asymptotic preserving time-splitting schemes

Time-splitting schemes for the Schrödinger equation have been widely used for applications that range from quantum transport to wave propagation in random media [26, 125, 151]. In this section, we discuss how the aforementioned asymptotic regimes (random phase modulation, loss of coherence, and radiative transfer regime) can be captured by a time-splitting scheme. It is a natural question to wonder how can we select sufficiently large stepsizes to reduce computational expenses while still accurately recovering the desired limiting regimes.

In order to define the splitting scheme for (2.3), consider a fixed stepsize $h > 0$, $n \in \mathbb{N}^*$, and let

$$ \mathcal{A}_\epsilon = \frac{i\epsilon^s}{2} \Delta_x, \quad \text{and} \quad \mathcal{B}_{n,h}^\epsilon(\mathbf{x}) = \frac{1}{\epsilon^{s-(1-\varrho)/2}} \int_{(n-1)h}^{nh} V\left(\frac{s}{\epsilon^{s+\varrho}}, \frac{\mathbf{x}}{\epsilon^s}\right) ds. $$

The standard first-order semi-discrete split-step scheme then reads

$$ \Psi_n^\epsilon = e^{h\mathcal{A}_\epsilon} e^{-i\mathcal{B}_{n,h}^\epsilon} \Psi_{n-1}^\epsilon, \quad \text{for } n = 1, 2, \dots, \quad \text{with} \quad \Psi_{0,\epsilon}(\mathbf{x}) = \phi_{0,\epsilon}(\mathbf{x}). \qquad (2.22) $$

Following the approach of [157, 194] for instance, the order in time of this scheme can be improved by taking other combinations of exponentials of operators. Note that the spatial variable $\mathbf{x}$ is not
discretized as the semi-group $e^{hA_\epsilon}$ can be evaluated at low cost with high accuracy using spectral methods, such as Fast Fourier Transforms (FFT), provided the spatial stepsize $\Delta x$ satisfies $\Delta x \ll \epsilon^s$. Absorbing boundary conditions are necessary since the wave function solves the Schrödinger equation in the entire space, and boundary conditions compatible with FFT can be found in [100].

Splitting techniques are efficient numerical methods for solving Schrödinger equations under the semi-classical limit. This scheme is always stable since it preserves the $L^2$ norm w.r.t. the variable $x$. For short-range correlations, and under the transport scaling ($s = 1$ and $\varrho = 0$), it has been shown in [23] that the scheme accurately captures (in the sense of a convergence in probability) the radiative transfer regime if the stepsize $h$ is of the form $\epsilon h'$, where $h' \ll 1$ is small but independent of $\epsilon$. It worth mentioning that a larger stepsize would still lead to a radiative transfer regime, but with an incorrect scattering cross-section (the one given by (2.20) for $\varrho > 0$ in Thm. 2.5.1 instead of (2.21)). The constraint $h = \epsilon h'$ is already a great improvement on the consistency constraint to capture the statistical behavior of the wave function. In fact, for a given realization of the random fluctuation $V$, the local error estimate between the exact solution and its approximation after one single step is readily shown to be

$$ \|\phi_{\epsilon}(h) - \Psi_{1}^{\epsilon}\|_{L^2(\mathbb{R}^d)} \le \frac{h^2}{\epsilon} \|[\epsilon \Delta_{\mathbf{x}}, \sqrt{\epsilon}V]e^{-\frac{i\epsilon}{2}\Delta_{\mathbf{x}}\phi_{0,\epsilon}}\|_{L^2(\mathbb{R}^d)} + \frac{h^2}{\epsilon} \|\phi_{0,\epsilon}\|_{L^2(\mathbb{R}^d)} \sup_{t,\mathbf{x}} |V(t, \mathbf{x})|^2, $$

where $[\cdot, \cdot]$ denotes the commutator between two operators. Assuming $\phi_{0,\epsilon}$ and $V$ smooth enough, the local error is then controlled by $Ch^2\epsilon^{-3/2}$, and the global error by $h\epsilon^{-3/2}$, requiring $h \ll \epsilon^{3/2}$ for the scheme to be consistent in the pathwise sense. This condition is considerably more restrictive than $h = \epsilon h'$ required to capture the transport regime in probability (in a statistical sense).

In what follows, we discuss the constraints on the stepsize $h$ to capture the three asymptotic regimes described above. For most of our parameter configurations, all three regimes are captured for any $h > 0$, independently of $\epsilon$. Such a result is a consequence of several factors. First, we only need to approximate the behavior of the wave function in the statistical sense, not in the pathwise sense, as the above asymptotic regimes are derived in a statistical sense and not w.r.t a particular realization of the fluctuations. Second, the random medium fluctuations exhibit multiscale properties allowing an accurate description of the dynamics at larger scales in the statistical sense.

## 2.6.1 The phase modulation

To capture the behavior of (2.6), a discretized version of the wave function is considered:

$$ \zeta_n^\epsilon(\mathbf{k}) = \frac{1}{\epsilon^{d(s-s_c)}} \hat{\Psi}_n^\epsilon \left( \frac{\mathbf{k}}{\epsilon^{s-s_c}} \right) e^{inh \frac{|\mathbf{k}|^2}{2\epsilon^{s-2s_c}}}, $$

where $\hat{\Psi}_n^\epsilon(\mathbf{k})$ stands for the Fourier transform of $\Psi_n^\epsilon(\mathbf{x})$ defined by (2.22). We have the following result for

$$ s = 1/(2\mathfrak{H}_{\varrho}), $$

with $\mathfrak{H}_{\varrho}$ given by (2.5).

**Theorem 2.6.1** For any $\mathbf{k} \in \mathbb{R}^d$, $h > 0$, $n > 0$, and $s_c \in [0, s]$, the family $(\zeta_n^\epsilon(\mathbf{k}))_\epsilon$ converges in law, as $\epsilon \to 0$, to

$$ \zeta(nh, \mathbf{k}) = \hat{\phi}_0(\mathbf{k}) \exp \left( i \sqrt{D_{n,h}(\mathbf{k})} W_{\mathfrak{H}}(nh) \right). $$

Here, $W_{\mathfrak{H}}$ is a standard fractional Brownian motion with Hurst index $\mathfrak{H}$, and $D_{n,h}$ satisfies

$$ \lim_{h \to 0} D_{[T/h],h}(\mathbf{k}) = D(\alpha, \beta, \mathbf{k}), \qquad (2.23) $$

given by (2.7) in the case $\beta = 1/2$, $\varrho = 0$, and $s_c = s$, and $D_{n,h}$ is given by (2.8) otherwise.

This result shows that the proposed time-splitting scheme captures the correct behavior as $\epsilon \to 0$ for any stepsize $h > 0$, independently of $\epsilon$, when at least one of the three following conditions hold: $\beta < 1/2$, $\varrho > 0$, $s_c < s$. When $\beta = 1/2$, $\varrho = 0$, and $s_c = s$, the condition $h \ll 1$ is required but independently of $\epsilon$. As already mentioned, such a result holds only because the statistical behavior of the wave function is of interest.
## 2.6.2 Loss of coherence

Let us recall that after the onset the random phase modulation the loss of coherence phenomena are described through two different behaviors. When $s = s_c = 1/(\theta + 1)$ (for $\beta < 1/2$, or $\varrho > 0$) the Wigner transform converges in law to the solution of a stochastic fractional heat equation (see Thm. 2.4.1). When $0 < s_c < s$, with $s_c$ given by (2.16), homogenization takes place and the Wigner transform converges in probability to the solution of a fractional heat equation (see Thm. 2.4.2).

These two behaviors are captured using the following discretized version of the Wigner transform:

$$W_{n,h}^{\varepsilon}(\mathbf{x}, \mathbf{k}) = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d \times S} dy \nu(d\zeta) e^{i\mathbf{k} \cdot \mathbf{y}} \Psi_n^{\varepsilon}(\mathbf{x} - \varepsilon^{s-s_c} \frac{\mathbf{y}}{2}, \zeta) \overline{\Psi_n^{\varepsilon}(\mathbf{x} + \varepsilon^{s-s_c} \frac{\mathbf{y}}{2}, \zeta)},$$

where $\Psi_n^{\varepsilon}(\mathbf{x})$ is defined by (2.22), with initial condition (2.9).

**Theorem 2.6.2** For either $\varrho > 0$ or $\beta < 1/2$, and $s_c = s = 1/(1 + \theta)$, where $\theta \in (0, 1)$ is defined by (2.12), and for all $n > 0$ and $h > 0$, $(W_{n,h}^{\varepsilon})_{\varepsilon}$ converges in law on $L^2(\mathbb{R}^{2d})$, as $\varepsilon \to 0$, to a limit $W(nh, \mathbf{x}, \mathbf{k})$ defined by

$$W(t, \mathbf{x}, \mathbf{k}) = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} dq \widehat{W}_0^{\mathbf{k}}(\mathbf{x}, \mathbf{q}) \exp \left( i\mathbf{k} \cdot \mathbf{q} + i \int_{\mathbb{R}^d} \mathcal{B}(t, d\mathbf{p}) e^{i\mathbf{p} \cdot \mathbf{x}} (e^{-i\mathbf{q} \cdot \mathbf{p}/2} - e^{i\mathbf{q} \cdot \mathbf{p}/2}) \right), \quad (2.24)$$

which is the unique solution of the stochastic fractional heat equation (2.14). In (2.24), $\widehat{W}_0^{\mathbf{k}}$ stands for the Fourier transform of $W_0$ defined by (2.18) with respect to the $\mathbf{k}$-variable.

As already explained, the remaining case $s = s_c$, $\varrho = 0$ and $\beta = 1/2$ is included in Thm. 2.6.1. The splitting scheme captures again the correct asymptotics for any $h > 0$. In this critical case, the limiting Wigner transform is random because the wave function does not propagate enough to average out the stochastic effects. This self-averaging effect is observed for longer propagation scale $s$ and the splitting scheme is again able to capture this effect.

**Theorem 2.6.3** Let $s_c = (1-s)/\theta < s$, where $\theta$ is defined by (2.12). Then, for all $n > 0$, $h > 0$, $(W_{n,h}^{\varepsilon})_{\varepsilon}$ converges in probability in $L^2_w(\mathbb{R}^{2d})$, as $\varepsilon \to 0$, to a limit $W(nh, \mathbf{x}, \mathbf{k})$ defined by

$$W(t, \mathbf{x}, \mathbf{k}) = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} dq e^{i\mathbf{k} \cdot \mathbf{q} - \sigma(\theta) |\mathbf{q}|^\theta} \widehat{W}_0^{\mathbf{k}}(\mathbf{x}, \mathbf{q}),$$

which is the unique solution of (2.17) uniformly bounded in $L^2(\mathbb{R}^{2d})$.

In this regime as well, the splitting scheme captures the correct asymptotic behavior for any $h > 0$.

## 2.6.3 Radiative transfer regime

In this regime ($s = 1$ and $s_c = 0$), there is some requirements on the stepsize $h$. The following result guarantees that the time-splitting scheme (2.22), for a stepsize $h \ll 1$ independent of $\varepsilon$, leads to the correct limit for $\varrho > 0$.

**Theorem 2.6.4** For any $\varrho > 0$, $n > 0$, $h > 0$, the family $(W_{n,h}^{\varepsilon})_{\varepsilon}$ converges in probability on $L^2_w(\mathbb{R}^{2d})$, as $\varepsilon \to 0$, to a limit defined by

$$W_{n,h}(\mathbf{x}, \mathbf{k}) = \frac{1}{(2\pi)^{2d}} \int_{\mathbb{R}^d} dy dq e^{i(\mathbf{x} \cdot \mathbf{y} + \mathbf{k} \cdot \mathbf{q}) + \psi_{n,h,\theta}(\mathbf{y}, \mathbf{q})} \widehat{W}_0(\mathbf{y}, \mathbf{y} + t\mathbf{q}).$$

Here,

$$\psi_{n,h,\theta}(\mathbf{y}, \mathbf{q}) = h \sum_{j=1}^{n} \psi(\mathbf{y} + nh(1 - j/n + 1/n)\mathbf{q}),$$

with

$$\psi(\mathbf{q}) = \frac{2}{\mu(2\pi)^d} \int_{\mathbb{R}^d} d\mathbf{p} \frac{a(\mathbf{p})}{|\mathbf{p}|^{d+\theta}} (e^{i\mathbf{p} \cdot \mathbf{q}} - 1),$$

so that

$$\lim_{h \to 0} W_{[T/h],h}(\mathbf{x}, \mathbf{k}) = W(T, \mathbf{x}, \mathbf{k}),$$

which is the unique solution to (2.19) uniformly bounded in $L^2(\mathbb{R}^{2d})$.
When $\rho = 0$ the situation is different. The transport term and momentum variables interact, which leads to another expression of the collision operator (see (2.20) versus (2.21)). This situation corresponds to the one considered in [23] under random media with short-range correlations, and for which one needs $h = \varepsilon h'$, with $h' \ll 1$ independent of $\varepsilon$ to recover the correct scattering operator. Without this restriction the scattering kernel (2.20) is obtained instead of (2.21).

## 2.7 Perspectives

The Schrödinger equation exhibits intriguing effective multiscale properties when coupled with a slowly decorrelating random potential. Capturing numerically these effects is also an interesting question as these properties seem to influence the performance of the numerical schemes. The perspectives in these directions are twofolds. First, these time-splitting schemes have not been implemented so far. It would be interesting to investigate if all the scaling regimes, described theoretically in this chapter, can effectively be captured numerically, and evaluate how these schemes are sensitive to the time stepsize.

In [15] under short-range correlations, or [105] where the authors start directly from a white-noise model for the random fluctuations, the long-time behavior of the wave function itself, and not the Wigner transform, is investigated. The asymptotic wavefield can be described through a Gaussian random field with a correlation function related to the radiative transfer model (2.19). A full description of the correlation function holds significant importance in the context of wave propagation, particularly for imaging problems through strongly scattering media relying on incoherent recorded signals. While the quadratic nature of the Wigner transform does not capture the slow phase of the wavefield (as it corresponds to the energy density) the direct analysis of the wavefield enables to extract information about the phase itself, as illustrated in Chap. 1 through the slow angle variable. It is well known in the imaging community that information carried by the phase is highly relevant for inversion algorithms. The second perspective consists in describing the loss of coherence phenomena on the wavefield itself and investigate whether it exhibits Gaussian statistics or not.
Chapter 3

Wave propagation in random media with long-range correlations

Analysis of physical measurements reveals that propagation media may exhibit perturbations with long-range dependencies [49,60,117,172,188]. This has sparked interest in a mathematical understanding of how waves propagate through multiscale media. Multiscale random media with long-range correlations find applications in modeling the heterogeneous earth crust, the turbulent atmosphere, and also biological tissues for instance. While there exists a large literature on wave propagation in heterogeneous media which vary on a well-defined microscale [70], the study of wave behavior in multiscale random media remains an area with many open questions. In order to be efficient, imaging and communication algorithms necessitate a deep understanding on how waves are affected by the rough medium fluctuations. In view of its potential for applications, mathematical descriptions of wave propagation in multiscale random media with long-range correlations has attracted a lot of interest over the last decade [15, 84, 153, 154, 191].

The full scalar wave equation is considered in this chapter with a continuous multiscale medium characterized by a stochastic process with long-range correlations:

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

with null initial conditions, $p(t = 0) = \partial_t p(t = 0) = 0$, so that the system is supposed to be initially at rest. Here, $p$ represents the acoustic pressure, and $c$ the wave-speed profile. The coordinate $z$ represents the main propagation axis, and the coordinate $\mathbf{x} \in \mathcal{D}$ represents the transverse variable (see Fig. 3.1). In this chapter, $\mathcal{D} = (0, d)$ is considered in Sect. 3.1 to define a planar waveguide with

<figure><img src="image_5.png" /><figcaption>Fig. 3.1: Illustration of the wave propagation model. The source $F$ generates a wave with profile $f$, which propagates along the $z$-axis. The section $z \in [0, L_z]$ is randomly heterogeneous with long-range correlations. Our objective is to characterize the pulse as it emerges at the end of the random section ($z = L_z$).</figcaption></figure>
a bounded cross-section, and $\mathcal{D} = \mathbb{R}^2$ in Sect. 3.2 and 3.3 to model open propagation media. The forcing term $\mathbf{F}(t, z, \mathbf{x})$ is given by

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

where $\mathbf{e}_z$ is the unit vector pointing in the z-direction. This term models a source located in the plan $z = L_S$, emitting a wave of the form $f(t, \mathbf{x})$. The wave-speed profile is assumed to be of the form

$$\frac{1}{c^2(z, \mathbf{x})} := \frac{1}{c_0^2} (1 + \nu(z, \mathbf{x}) \mathbf{1}_{(0, L_z)}(z)),$$

where $\nu$ models the fluctuations of the propagation medium. The aim of this chapter is to present mathematical descriptions of the wavefield as it exits the random section at $z = L_z$ in three different scenarios, emphasizing the crucial role of the fluctuations with long-range correlations.

As described in Chap. 2 for the Schrödinger equation, random media with long-range correlations affect the wave propagation in a way that differs greatly from the corruption caused by media fluctuating on a well-defined microscale (short-range correlations) [70]. Wave propagation in random media with long-range correlations has already been considered in one-dimensional propagation media [84,153], or open media under the paraxial approximation [15,68,69]. These works have demonstrated that wave propagating in such random media exhibits stochastic effects at various propagation scales. Unlike for propagation media with short-range correlations, all the stochastic effects do not appear on the same propagation scale. Perturbations with long-range correlations first induce a random phase modulation on the waves driven by a single standard fractional Brownian motion, which does not depend on the frequency band [15,152]. For larger propagation distances, the random phase modulation starts to oscillate very fast up to produce diffusion phenomena and affects the energy propagation [84]. This picture follows the ones presented in Chap. 1 and 2, where the phase modulation is exhibited by a non-central limit theorem, and the diffusion phenomena by approximation-diffusion theorems.

This chapter is organized as follows. First, we consider in Sect. 3.1 the context of a waveguide supporting a discrete number of modes. This section considers $\mathcal{D} = (0, d)$ as transverse section and constitutes an intermediate situation between a one-dimensional propagation medium and the full ambient physical space (obtained for $\mathcal{D} = \mathbb{R}^2$). In this section, an analysis similar to the non-central limit theorem is developed to exhibit phase modulations driven by a unique fractional Brownian motion for each propagating modes in a high-frequency limit. In Sect. 3.2 a similar approach is developed for $\mathcal{D} = \mathbb{R}^2$, but this time two phenomena are captured at once in a single high-frequency limit: a non-central limit type result, as well as the paraxial approximation. This approximation consists in describing the wave propagation along a privileged axis, and has been extensively studied and used in applications (see [10, 24, 54, 201] for instance). Under suitable assumptions on physical parameters, this approximation can greatly simplify the description of propagation phenomena as well as their numerical simulations [151]. The paraxial approximation together with the fractional white-noise approximation lead to a model of fractional Itô-Schrödinger equation driven by a fractional noise, not necessarily Gaussian, to describe the impact of the medium fluctuations on the propagating waves. The evolution variable of this fractional Itô-Schrödinger equation corresponds to the one of the main propagation axis. Under the setting of this section, the backscattering effects are shown to be negligible, and the energy carried by the propagating wave is conserved. Finally, in Sect. 3.3, under the simpler setting of a randomly layered medium (that is $\nu(z, \mathbf{x}) = \nu(z)$), but still with $\mathcal{D} = \mathbb{R}^2$, the waves are allowed to propagate on a larger propagation distance than the one of the non-central limit theorem exhibiting the phase modulation. Typically, despite considering long-range correlations, we let the waves propagate over distances leading to diffusion effects under short-range correlations. At this scale, the random phase modulation oscillates very fast up to produce effective diffusion phenomena. The backscattering effects become now nonnegligible and the energy is now affected by the medium fluctuations. In this section as well, the diffusion phenomena and paraxial approximation are obtained under a single high-frequency limit. These phenomena can be described through a deterministic fractional paraxial equation, wherein the order of the fractional time derivative is directly determined by the decay rate of the correlation function of the random fluctuations.

This chapter covers the results obtained in [10, 12, 17] of the bibliography list on page 6.
## 3.1 Random waveguides

Considering $\mathcal{D} = (0, d)$, Eq. (3.1) describes scalar wave propagation in a planar waveguide with propagation axis given by the z-direction (see Fig. 3.1). To simplify the presentation, only Dirichlet boundary conditions are here considered to complete the system on $\partial \mathcal{D}$. A similar analysis can be applied to other boundary conditions. The wave equation (3.1) being linear in $p$, its solution can be expressed as a superposition of monochromatic waves through the Fourier transform:

$$ \hat{p}(\omega, \mathbf{x}, z) = \int p(t, \mathbf{x}, z)e^{i\omega t} dt \quad \text{and} \quad p(t, \mathbf{x}, z) = \frac{1}{2\pi} \int \hat{p}(\omega, \mathbf{x}, z)e^{-i\omega t} d\omega. $$

The wavefield $\hat{p}(\omega, \mathbf{x}, z)$ satisfies, in the randomly perturbed section $z \in (0, L_z)$, the following Helmholtz equation (time-harmonic wave equation)

$$ \partial_{zz}^2 \hat{p}(\omega, \mathbf{x}, z) + \partial_{xx}^2 \hat{p}(\omega, \mathbf{x}, z) + k_{\omega}^2 (1 + \varepsilon \nu(z, \mathbf{x})) \hat{p}(\omega, \mathbf{x}, z) = 0, $$

where $k_{\omega} = \omega/c_0$ is the wavenumber. Here, the medium fluctuations are assumed to be small, of order $\varepsilon \ll 1$ (weak scattering regime). To observe non-trivial cumulative stochastic effects on the propagating waves a propagation distance of order

$$ L_z \sim \frac{L}{\varepsilon^s} \gg 1 $$

is required. The central wavelength is of order one in this problem ($\lambda \sim 1$), that is small compared to the propagation distance ($\lambda \ll L_z$), to place ourselves in a high-frequency regime. The propagation scale parameter $s$ represents how long the propagation distance needs to be to observe the onset of effective random effects on the waves. Assuming

$$ \mathbb{E}[\nu(z + z', \mathbf{x})\nu(z', \mathbf{y})] \underset{|z| \to \infty}{\sim} \frac{R_0}{|z|^\gamma} R(\mathbf{x}, \mathbf{y}), \quad (3.2) $$

describing the long-range correlation property for the medium fluctuations, we need to take

$$ s = \frac{1}{\mathfrak{H}} < 2 \quad \text{with} \quad \mathfrak{H} = 1 - \frac{\gamma}{2} \in (1/2, 1). $$

With this choice, $L/\varepsilon^s \ll L/\varepsilon^2$, where the latter corresponds to the standard scaling to observe non-trivial stochastic effects under fluctuations with short-range correlations. While not being restrictive, in addition to the boundedness of $\nu$, some technical assumptions need to be set on $\nu$ to be compatible with a moment technique. Basically, $\nu$ is a nonlinear function (with Hermite rank 1) of a Gaussian filed satisfying a relation similar to (3.2). The precise description of these assumptions are omitted here, and we refer to [99] for the complete assumptions. The requirement on the Hermite rank could be relaxed to arbitrary rank, as proposed in [154], with an adaptation of the amplitude $\varepsilon$ of the fluctuations. Also, the width of the waveguide $d$ is assumed to be of the order of the central wavelength $\lambda$ (of order 1), so that only a discrete number of mode is supported by the waveguide. The wavefield $\hat{p}$ can then be expressed as

$$ \hat{p}(\omega, \mathbf{x}, z) = \hat{p}_a(\omega, \mathbf{x}, z) + \hat{p}_e(\omega, \mathbf{x}, z), $$

with

$$ \hat{p}_a(\omega, \mathbf{x}, z) := \sum_{j=1}^{N(\omega)} \frac{\hat{a}_j(\omega, z)}{\sqrt{\beta_j(\omega)}} e^{i\beta_j(\omega)z} \phi_j(\mathbf{x}) + \sum_{j=1}^{N(\omega)} \frac{\hat{b}_j(\omega, z)}{\sqrt{\beta_j(\omega)}} e^{-i\beta_j(\omega)z} \phi_j(\mathbf{x}), \quad (3.3) $$

and

$$ \hat{p}_e(\omega, \mathbf{x}, z) := \sum_{j \ge N(\omega)+1} \hat{p}_j(\omega, z) \phi_j(\mathbf{x}). $$

Here, $\hat{p}$ is expanded according to the spectral decomposition associated to the transverse Laplacian $-\partial_{xx}^2$ equipped with Dirichlet boundary conditions on $\partial \mathcal{D}$. Its spectrum is composed of a countable number of positive eigenvalues

$$ \lambda_j = \frac{j^2 \pi^2}{d^2}, \quad j \ge 1, $$
<figure><img src="image_6.png" /><figcaption>Fig. 3.2: Illustration of the right-going and left-going propagating mode amplitudes $\hat{a}_j(\omega, z)$ and $\hat{b}_j(\omega, z)$. The source generates the probing wave $\hat{a}_0(\omega)$, the reflected wave is $\hat{b}(\omega, 0)$ and the transmitted wave, which is our interest, is $\hat{a}(\omega, L/\epsilon^s)$. Note that there is no energy coming from the right-hand-side of the waveguide $z > L/\epsilon^s$.</figcaption></figure>

with eigenvectors

$$ \phi_j(\mathbf{x}) := \sqrt{\frac{2}{d}} \sin\left(\frac{j\pi\mathbf{x}}{d}\right) \quad j \ge 1, $$

providing an orthonormal basis of $L^2(\mathcal{D})$. The family $(\beta_j(\omega))_{j \ge 1}$ is defined as

$$ \beta_j(\omega) := \sqrt{k_{\omega}^2 - \lambda_j} \quad \text{for } j \in \{1, \dots, N(\omega)\}, $$

where $N(\omega)$ is the integer such that $\lambda_{N(\omega)} \le k_{\omega}^2 < \lambda_{N(\omega)+1}$, and represents the number of propagating modes. For our planar waveguide model, one has

$$ N(\omega) = \left[ \frac{\omega d}{\pi c} \right], $$

where $[\cdot]$ stands for the integer part. In (3.3), the terms $\hat{a}_j$ and $\hat{b}_j$ represent the amplitudes of the right- and left-propagating modes respectively, as they corresponds to the oscillatory terms $e^{\pm i\beta_j(\omega)z}$ w.r.t. the variable of the main propagation axis (see Fig. 3.2 for an illustration). Here, there is exactly $2N(\omega)$ propagating modes. The terms $\hat{p}_j$ in $\hat{p}_e$ represent the amplitudes of the evanescent modes and corresponds to non-oscillatory terms w.r.t. the z-variable. Although the evanescent modes are negligible after propagating over large propagation distances, their influences can be incorporated to the mode-coupling mechanism of the propagating modes $\hat{a}_j$ and $\hat{b}_j$ following the strategy of [79]:

$$ \frac{d}{dz} \begin{bmatrix} \hat{a}(\omega, z) \\ \hat{b}(\omega, z) \end{bmatrix} = \left[ \varepsilon \mathbf{H}(\omega, z) + \varepsilon^2 \mathbf{G}(\omega, z) \right] \begin{bmatrix} \hat{a}(\omega, z) \\ \hat{b}(\omega, z) \end{bmatrix}, \quad (3.4) $$

completed with the boundary conditions

$$ \hat{a}_j(\omega, 0) = \hat{a}_{j,0}(\omega) \quad \text{and} \quad \hat{b}_j(\omega, L/\epsilon^s) = 0. \quad (3.5) $$

The second relation means that no wave is coming from the right-hand side of the random section, and the $\hat{a}_{j,0}(\omega)$'s are the propagating mode amplitudes generated by the source term. The coupling matrices $\mathbf{H}$ and $\mathbf{G}$ are defined by:

$$ \mathbf{H}(\omega, z) := \begin{bmatrix} \mathbf{H}^a(\omega, z) & \overline{\mathbf{H}^b(\omega, z)} \\ \mathbf{H}^b(\omega, z) & \overline{\mathbf{H}^a(\omega, z)} \end{bmatrix} \quad \text{and} \quad \mathbf{G}(\omega, z) := \begin{bmatrix} \mathbf{G}^a(\omega, z) & \overline{\mathbf{G}^b(\omega, z)} \\ \mathbf{G}^b(\omega, z) & \overline{\mathbf{G}^a(\omega, z)} \end{bmatrix}, $$

with

$$ \mathbf{H}_{jl}^a(z) := \frac{ik_{\omega}^2}{2} C_{jl}(z) e^{i(\beta_l - \beta_j)z}, \quad \mathbf{H}_{jl}^b(z) := -\frac{ik_{\omega}^2}{2} C_{jl}(z) e^{i(\beta_l + \beta_j)z} $$

$$ \mathbf{G}_{jl}^a(z) := \frac{ik_{\omega}^4}{4} \sum_{l' \ge N+1} \int_{-\infty}^{\infty} C_{jl'}(z) C_{l'l}(z+u) e^{i\beta_{l'}(z+u) - i\beta_j z - \beta_{l'}|u|} du, $$

$$ \mathbf{G}_{jl}^b(z) := -\frac{ik_{\omega}^4}{4} \sum_{l' \ge N+1} \int_{-\infty}^{\infty} C_{jl'}(z) C_{l'l}(z+u) e^{i\beta_{l'}(z+u) + i\beta_j z - \beta_{l'}|u|} du, $$
and where the coupling coefficient $C_{jl}(z)$ are defined by

$$C_{jl}(z) := \frac{1}{\sqrt{\beta_j \beta_l}} \int_D \nu(z, \mathbf{x}) \phi_j(\mathbf{x}) \phi_l(\mathbf{x}) d\mathbf{x}.$$

The matrix $\mathbf{H}$ describes the coupling between the propagating modes ($\mathbf{H}^a$ describes the coupling between the forward-going modes and backward-going modes themselves, while $\mathbf{H}^b$ describes the coupling between the forward- and backward-going modes). These coupling coefficients involve directly the medium fluctuations $\nu$. Without these fluctuations, there is no mode coupling and the solution to (3.4) is straightforward.

To account for the cumulative stochastic effects over the large propagation distance $L/\epsilon^s$, the propagating mode amplitudes are rescaled as

$$\hat{a}^{\epsilon}(\omega, z) = \hat{a}(\omega, z/\epsilon^s) \quad \text{and} \quad \hat{b}^{\epsilon}(\omega, z) = \hat{b}(\omega, z/\epsilon^s).$$

The corresponding two-point boundary value problem (3.4) is not convenient to analyze directly. It is more convenient to recast this problem into a initial value problem by mean of a propagator matrix defined as the solution to

$$\frac{d}{dz} \mathbf{P}^{\epsilon}(\omega, z) = \left[ \frac{1}{\epsilon^{s-1}} \mathbf{H} \left( \omega, \frac{z}{\epsilon^s} \right) + \epsilon^{2-s} \mathbf{G} \left( \omega, \frac{z}{\epsilon^s} \right) \right] \mathbf{P}^{\epsilon}(\omega, z), \quad \text{with} \quad \mathbf{P}^{\epsilon}(\omega, 0) = \mathrm{Id}_{2N(\omega)},$$

so that for all $z \in [0, L]$,

$$\begin{bmatrix} \hat{a}^{\epsilon}(\omega, z) \\ \hat{b}^{\epsilon}(\omega, z) \end{bmatrix} = \mathbf{P}^{\epsilon}(\omega, z) \begin{bmatrix} \hat{a}^{\epsilon}(\omega, 0) \\ \hat{b}^{\epsilon}(\omega, 0) \end{bmatrix} \quad \text{and} \quad \begin{bmatrix} \hat{a}^{\epsilon}(\omega, L) \\ 0 \end{bmatrix} = \mathbf{P}^{\epsilon}(\omega, L) \begin{bmatrix} \hat{a}_0(\omega, 0) \\ \hat{b}^{\epsilon}(\omega, 0) \end{bmatrix},$$

according to (3.5). Due to the form of $\mathbf{H}(\omega, z)$ and $\mathbf{G}(\omega, z)$, the propagator can be expressed as

$$\mathbf{P}^{\epsilon}(\omega, z) = \begin{bmatrix} \mathbf{P}^{a, \epsilon}(\omega, z) & \mathbf{P}^{b, \epsilon}(\omega, z) \\ \mathbf{P}^{b, \epsilon}(\omega, z) & \mathbf{P}^{a, \epsilon}(\omega, z) \end{bmatrix},$$

where $\mathbf{P}^{a, \epsilon}$ describes the coupling mechanisms between the right-going modes (resp. left-going modes) with themselves, while $\mathbf{P}^{b, \epsilon}$ describes the coupling mechanisms between the right-going and left-going modes.

Due to the form of the propagator $\mathbf{P}^{\epsilon}$, one can restrict the asymptotic analysis to the one of the two blocks:

$$\mathcal{P}^{\epsilon}(\omega, z) = \begin{bmatrix} \mathbf{P}^{a, \epsilon}(\omega, z) \\ \mathbf{P}^{b, \epsilon}(\omega, z) \end{bmatrix} \in \mathcal{M}_{N(\omega)}(\mathbb{C}) \times \mathcal{M}_{N(\omega)}(\mathbb{C}),$$

where $\mathcal{M}_{N(\omega)}(\mathbb{C})$ stands for the set of $N(\omega) \times N(\omega)$ matrices with complex coefficients. The asymptotic behavior of the propagator $\mathcal{P}^{\epsilon}$ can be described as follows.

**Theorem 3.1.1** For all $z \in [0, L]$, the family $(\mathcal{P}^{\epsilon}(\omega, z))_{\epsilon}$ converges in law in $\mathcal{M}_N(\mathbb{C}) \times \mathcal{M}_N(\mathbb{C})$ to

$$\begin{bmatrix} D(\omega, z) \\ \mathbf{0} \end{bmatrix},$$

with

$$D(\omega, z) = \mathrm{diag}(e^{i\sigma_1 s_1(\omega)W_{s_1}(z)}, \dots, e^{i\sigma_N s_N(\omega)W_{s_N}(z)}),$$

where $W_{s_j}$ is a standard fractional Brownian motion with Hurst index $\mathfrak{H}$,

$$\sigma_{j,s_j}(\omega) = \frac{k_{\omega}^2}{2\beta_j(\omega)} \sqrt{\frac{R_j}{\mathfrak{H}(2\mathfrak{H} - 1)}}, \quad (3.6)$$

and

$$R_j = \iint_{\mathcal{D} \times \mathcal{D}} R(\mathbf{x}, \mathbf{x}') \phi_j^2(\mathbf{x}) \phi_j^2(\mathbf{x}') d\mathbf{x} d\mathbf{x}'.$$
The proof of Thm. 3.1.1 is based on a moment technique following the ideas of [15, 66]. Let us clarify the meaning of this result. First, Thm. 3.1.1 implies that the first significant stochastic effects affecting the wave propagation take place for $s = 1/\mathfrak{H} < 2$ which is in contrast with the classical results under short-range correlations with $s = 2$ (see [70, Chap. 20]). The second one concerns the convergence of $\mathbf{P}^{b,\varepsilon}(\omega, z)$ to the null matrix in probability, meaning that the coupling mechanisms between the right- and left-going modes are negligible in the high-frequency limit $\varepsilon \to 0$. In other words, the backscattering is negligible for $\varepsilon$ small, and the backscattered amplitudes $\hat{b}_j^\varepsilon(\omega, z)$ (and then the backscattered wavefield) tend to 0. Third, the convergence in distribution of $\mathbf{P}^{a,\varepsilon}(\omega, z)$ to a diagonal matrix means that the coupling mechanisms between two different right-going modes are also negligible in the limit $\varepsilon \to 0$. Finally, the propagating modes are only affected by mode- and frequency-dependent phase modulations, but driven by the same fractional Brownian motion, which does not depend on the frequency $\omega$. This effect of mode-dependent phase modulations driven by a single fractional Brownian motion has already been observed in [15, Thm 1.2] for the random Schrödinger equation with long-range correlations and stated in Thm 2.2.1 of Chap. 2.

A version of Thm. 3.1.1 involving multiple frequencies can be derived following the same strategy, and allows to describe the role of these stochastic effects in the time domain. Considering now a narrowband source term,

$$f^{\varepsilon}(t, \mathbf{x}) = f(\varepsilon^s t) e^{-i\omega_0 t} \Psi(\mathbf{x}),$$

with carrier frequency $\omega_0$ (of order 1) and bandwidth of order $\varepsilon^s$, we have the following result.

**Theorem 3.1.2** For all $j \in \{1, \dots, N(\omega_0)\}$, let us consider the projection

$$p_{j,pr}^{\varepsilon}(t, L) = e^{-i\beta_j(\omega_0)(L/\varepsilon^s - L_S)} e^{i\omega_0 t/\varepsilon^s} \langle p_{pr} \left( \frac{t}{\varepsilon^s}, \cdot, \frac{L}{\varepsilon^s} \right), \phi_j \rangle_{L^2(D)}.$$

The family $(p_{1,pr}^{\varepsilon}(\cdot, L), \dots, p_{N(\omega_0),pr}^{\varepsilon}(\cdot, L))_{\varepsilon}$ converges in law in $C([-T, T], \mathbb{C}^{N(\omega_0)})$, as $\varepsilon \to 0$, to $(p_{j,pr}^{0}(\cdot, L))_{j \in \{1, \dots, N(\omega_0)\}}$, where

$$p_{j,pr}^{0}(t, L) = \frac{1}{2} e^{i\sigma_{j,\mathfrak{H}}(\omega_0)W_{\mathfrak{H}}(L)} f(t - \beta'_{j}(\omega_0)L) \langle \phi_j, \Psi \rangle_{L^2(D)}.$$

Here, $W_{\mathfrak{H}}$ is a standard fractional Brownian motion with Hurst index $\mathfrak{H} \in (1/2, 1)$, and $\sigma_{j,\mathfrak{H}}(\omega_0)$ is defined by (3.6).

In this result the wavefield is projected over the propagating modes. The reason is that the pulsewidth of the source and the propagation distance are of the same order ($\sim 1/\varepsilon^s$), so that the propagating modes overlap, there is no modal dispersion in this situation. By projecting over the propagating modes, we can compensate the rapid phases $e^{i\beta_j(\omega)L/\varepsilon^s}$ and $e^{-i\omega_0 t/\varepsilon^s}$ (inherent of a high-frequency problem) mode by mode and pass to the limit in $\varepsilon$. Roughly speaking, the transmitted wave can be formally described as

$$p_{tr} \left( \frac{t}{\varepsilon^s}, \mathbf{x}, \frac{L}{\varepsilon^s} \right) \underset{\varepsilon \to 0}{\sim} \frac{e^{-i\omega_0 t/\varepsilon^s}}{2} \sum_{j=1}^{N(\omega_0)} e^{i\sigma_{j,\mathfrak{H}}(\omega_0)W_{\mathfrak{H}}(L)} f(t - \beta'_{j}(\omega_0)L) e^{i\beta_j(\omega_0)(L/\varepsilon^s - L_S)} \phi_j(\mathbf{x}) \langle \phi_j, \Psi \rangle_{L^2(D)}.$$

This transmitted wave is therefore a superposition of modes, each of them is centered around its arrival time

$$t_j = \beta'_{j}(\omega_0)L,$$

$\beta'_{j}(\omega)$ being the derivative of $\beta_j(\omega)$ w.r.t. $\omega$, but also modulated by a mode-dependent and frequency-dependent random phase. Once again the randomness comes from the same fractional Brownian motion for all the propagating modes.

## 3.2 Paraxial approximation and fractional Itô-Schrödinger equation

In this section, a full ambient physical space is considered, and two particular asymptotic behaviors are described: the paraxial approximation and the fractional white-noise limit. In contrast with
Chap. 2, where we initially examine the paraxial approximation formally and then proceed with the analysis of the stochastic aspects separately, both of the limits are here obtained under the same limit. Considering both of these aspects under a single limit is much more challenging mathematically than considering two distinct limits.

The first aspect is the paraxial (parabolic) approximation, which is valid when the wave predominantly propagates in a specific direction (here the z-direction) and is sufficiently collimated around it. In the frequency domain, the 3D Helmholtz equation is reduced to 2D Schrödinger equation where the variable along the primary propagation axis takes the role of the evolution variable. This approximation offers a considerable advantage as it transforms a boundary value problem into an evolution problem with lower spatial dimensions. In homogeneous media, deriving the paraxial wave equation is relatively straightforward and involves an asymptotic expansion of the principal symbol of the operator describing the propagation (here the one of the scalar wave equation). However, the situation becomes significantly more complex in heterogeneous media since the interaction with the medium generates some backscattering and a coupling with the evanescent modes. To justify the paraxial approximation under such conditions, one must rely on specific features of the medium, for instance small amplitudes of the fluctuations [24], or oscillatory behaviors [10, 83]. The second limit corresponds to a non-central limit type result arising from the long-range dependencies of the propagation medium. In case of short-range correlations, the random medium fluctuations are then asymptotically statistically equivalent to a white-noise w.r.t. the main propagation axis. The resulting limiting model corresponds to an Itô-Schrödinger equation, as derived in [83], and studied mathematically in detail in [58].

This section depicts the rigorous simultaneous derivation of the paraxial and the (fractional) white-noise approximations in the context of random media with long-range correlations w.r.t. the z-direction. Heuristically, the limiting classical white noise obtained in [83] is replaced by a fractional white noise, leading to a fractional Itô-Schrödinger equation. From the mathematical viewpoint, this is a significantly more difficult problem than the ones addressed in [10, 83]. As already discussed earlier in this manuscript, the martingale techniques used in [10, 83] and standard approximation-diffusion theorems do not apply, and we make use of moment techniques inspired from [15, 66] which are fairly involved analytically. The proof also involves a fine analysis of the backscattering and the coupling between the propagating and evanescent modes. Note as well that the existence theory for the fractional Itô-Schrödinger equation is not trivial, and leads to some additional difficulties in the asymptotic theory.

### 3.2.1 Paraxial scaling regime

To introduce an appropriate parameter scaling regime under which the paraxial approximation is valid, some assumptions are required on the source term and the random medium fluctuations. Going back to the wave equation (3.1), under a dimensionless setting, the source term is assumed to be of the form

$$F(t, z, \mathbf{x}) := f\left(\frac{t}{\lambda}, \frac{\mathbf{x}}{r_0}\right) \delta(z - L_S)\mathbf{e}_z,$$

modeling a source located in the plane $z = L_S < 0$, and emitting a pulse in the z-direction with profile $f$ (see Fig. 3.1), central wavelength $\lambda$, and transverse width $r_0$. The wave-speed profile is assumed to be of the form:

$$\frac{1}{c^2(z, \mathbf{x})} := \begin{cases} \frac{1}{c_0^2} \left( 1 + \sigma \nu \left( \frac{z}{l_c}, \frac{\mathbf{x}}{l_c} \right) \right) & \text{if } z \in [0, L], \\ \frac{1}{c_0^2} & \text{if } z \in (-\infty, 0) \cup (L, \infty), \end{cases}$$

for $\mathbf{x} \in \mathbb{R}^2$, and where $c_0$ stands for the background (constant) velocity. The correlation length of the fluctuations $l_c$ can be seen as their typical scale of variation, and $\sigma$ characterizes the order of magnitude of these fluctuations. To introduce a high-frequency regime, the central wavelength $\lambda$ is assumed to be small compared to the propagation distance $L_z$:

$$\varepsilon^2 := \frac{\lambda}{L_z} \ll 1.$$
Strong interactions between the wave and the medium fluctuations are obtained by assuming

$$r_0 \sim l_c, \quad \text{and} \quad \frac{\lambda}{l_c} \sim \varepsilon.$$

We also assume

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

so that the Rayleigh length is of order the propagation distance, which is crucial to derive the paraxial approximation. In fact, the Rayleigh length is defined as the distance from the beam waist to the place where its cross-section is doubled, and in homogeneous media it is of order $r_0^2/\lambda$. For simplicity in the presentation, we set

$$L_z = L \sim 1, \quad \lambda = \varepsilon^2, \quad l_c = \varepsilon, \quad \text{and} \quad r_0 = \varepsilon. \qquad (3.7)$$

Both the paraxial approximation and the non-central limit type result (fractional white-noise approximation) are obtained by sending $\varepsilon$ to 0.

Finally, the strength of the fluctuations is assumed to be small (weak scattering regime)

$$\sigma = \varepsilon^s, \quad \text{with} \quad s = 2 - \frac{\gamma}{2}. \qquad (3.8)$$

Note that this approach is equivalent to the ones introduced in the previous section and Chap. 2 where the parameter $s$ was related to the propagation distance [15]. The parameter $\gamma \in (0, 1)$ still describes the decay rate of the correlation function

$$\mathbb{E}[\nu(z + z', \mathbf{x})\nu(z', \mathbf{y})] \underset{|z| \to \infty}{\sim} \frac{R_0}{|z|^{\gamma}} R(\mathbf{x} - \mathbf{y}).$$

As for the previous section, while not being restrictive, similar technical assumptions need to be set on $\nu$ to be compatible with a moment technique. The precise assumptions are omitted here, and we refer to [95] for their complete descriptions.

### 3.2.2 Main result and properties

Under the scaling (3.7) and (3.8), the original wave problem (3.1) becomes

$$\Delta p_{\varepsilon} - \frac{1}{c_0^2} \left(1 + \varepsilon^s V\left(\frac{z}{\varepsilon}, \frac{\mathbf{x}}{\varepsilon}\right) \mathbf{1}_{(0,L)}(z)\right) \partial_{tt}^2 p_{\varepsilon} = f\left(\frac{t}{\varepsilon^2}, \frac{\mathbf{x}}{\varepsilon}\right) \delta'(z - L_S).$$

The main result of this section concerns the asymptotic description of the pulse front exiting from the random section at $z = L$. It is defined by

$$p_L^\varepsilon(t, \mathbf{x}) = p\left(\lambda t + \frac{L - L_S}{c_0}, L, r_0 \mathbf{x}\right) = p\left(\varepsilon^2 t + \frac{L - L_S}{c_0}, L, \varepsilon \mathbf{x}\right),$$

and it corresponds to the wave observed on a time window of order the pulse width centered around the expected arrival time, and at the transverse scale of the source profile.

**Theorem 3.2.1** The family $(p_L^\varepsilon)_\varepsilon$ converges in law in the space $C^0((-\infty, +\infty), L^2(\mathbb{R}^2)) \cap L^2((-\infty, +\infty) \times \mathbb{R}^2)$, as $\varepsilon \to 0$, to a limit given by

$$p_L^0(t, \mathbf{x}) = \frac{1}{2\pi} \int e^{-i\omega t} \Psi_\omega(L, \mathbf{x}) d\omega,$$

where $\Psi_\omega$ is the unique pathwise solution to the fractional Itô-Schrödinger equation

$$d\Psi_\omega(z, \mathbf{x}) = \frac{i}{2k_\omega} \Delta_\mathbf{x} \Psi_\omega(z, \mathbf{x}) + ik_\omega \Psi_\omega(z, \mathbf{x}) dW_{\mathfrak{H}}(z, \mathbf{x}) = 0. \qquad (3.9)$$

Here, $k_\omega = \omega/c_0$, $W_{\mathfrak{H}}$ is a random mixture of fractional Brownian fields with Hurst index $\mathfrak{H}$ and covariance operator given by

$$\mathbb{E}[W_{\mathfrak{H}}(z_1, \mathbf{x}) W_{\mathfrak{H}}(z_2, \mathbf{y})] = \frac{R_0(\mathbf{x} - \mathbf{y})}{2\mathfrak{H}(2\mathfrak{H} - 1)} (z_1^{2\mathfrak{H}} + z_2^{2\mathfrak{H}} - |z_1 - z_2|^{2\mathfrak{H}}),$$

for all $z_1, z_2 \in [0, +\infty)$ and $\mathbf{x}, \mathbf{y} \in \mathbb{R}^2$.
The notion of solution for this fractional Itô-Schrödinger equation is made precise in the next section. As for the waveguide, due to a negligible backscattering the energy carried by the pulse is conserved at the end of the random section, that is

$$ \|p_L^0\|_{L^2((-\infty, +\infty) \times \mathbb{R}^2)} = \frac{1}{2} \|f\|_{L^2((-\infty, +\infty) \times \mathbb{R}^2)}. $$

Moreover, all the moments

$$ \mathbb{E} \left[ \left\langle \prod_{j_1=1}^{M_1} \Psi_{\omega_{1,j_1}}(z) \prod_{j_2=1}^{M_2} \overline{\Psi_{\omega_{2,j_2}}(z)}, \varphi \right\rangle_{L^2(\mathbb{R}^{2(M_1+M_2)})} \right], $$

for all $(M_1, M_2) \in \mathbb{N}^2$, distinct frequencies $(\omega_{i,j_i})_{(i,j_i) \in \{1,2\} \times \{1,\dots,M_i\}}$, and test function $\varphi \in L^2(\mathbb{R}^{2(M_1+M_2)})$, can be explicitly computed through classical Born series (see [120, Sect. 17.2]). These moments are important in imaging applications where they help to quantify the stability of reconstructions with respect to changes in the random medium [19, 20, 39, 86].

### 3.2.3 Stochastic integral and fractional Itô-Schrödinger equation

The stochastic integral w.r.t. a fractional noise $W_β$ in the limiting fractional Itô-Schrödinger equation (3.9) is of pathwise type. It is defined here according to the work of Zähle [211] and the approach of [156, 168]. To defined this integral we introduce the Banach space

$$ W^{\alpha}(0, L, B) := \left\{ \psi \in C^0([0, L], B) \text{ such that } \|\psi\|_{\alpha, B} < +\infty \right\}, $$

with

$$ \|\phi\|_{\alpha, B} := \sup_{z \in [0, L]} \left[ \|\phi(z)\|_B + \int_0^z \frac{\|\phi(z) - \phi(u)\|_B}{(z-u)^{\alpha+1}} du \right], $$

and where $B$ is a given Banach space. For $\beta \in (0, 1)$, we denote by $C^\beta([0, L], B)$ the Banach space of $\beta$-Hölder functions on $[0, L]$ with values in $B$, and equipped with

$$ \|\phi\|_{\beta, C, B} := \sup_{z \in [0, L]} \|\phi(z)\|_B + \sup_{0 \le v < u \le L} \frac{\|\phi(u) - \phi(v)\|_B}{(u-v)^\beta}. $$

Note that for $\alpha \in (0, 1/2)$, we have

$$ \|\phi\|_{\alpha, B} \le \left(1 + \frac{L^{1-2\alpha}}{1-2\alpha}\right) \|\phi\|_{1-\alpha, C, B} \quad \text{so that} \quad C^{1-\alpha}([0, L], B) \subset W^{\alpha}(0, L, B). \quad (3.10) $$

Now, for $\alpha \in (0, 1)$ and $z \in (0, L)$, we introduce, for a real-valued function $f$, the so-called Weyl's derivative:

$$ D_{0+}^{\alpha} \phi(z) := \frac{1}{\Gamma(1-\alpha)} \left[ \frac{\phi(z)}{z^{\alpha}} + \alpha \int_{0}^{z} \frac{\phi(z) - \phi(u)}{(z-u)^{\alpha+1}} du \right], $$

$$ D_{L^{-}}^{\alpha} \phi(z) := \frac{(-1)^{\alpha}}{\Gamma(1-\alpha)} \left[ \frac{\phi(z)}{(L-z)^{\alpha}} + \alpha \int_{z}^{L} \frac{\phi(z) - \phi(u)}{(u-z)^{\alpha+1}} du \right], $$

whenever these quantities are well-defined, and where $\Gamma$ is the standard Gamma function. Following the idea of [211], the generalized Stieljes integral of a function $\phi \in C^\nu([0, L], \mathbb{R})$ with respect to $\psi \in C^u([0, L], \mathbb{R})$, with $\nu + \mu > 1$, $\nu > \alpha$, and $\mu > 1 - \alpha$ is defined by

$$ \int_{0}^{L} \phi d\psi := (-1)^{\alpha} \int_{0}^{L} D_{0+}^{\alpha} \phi(u) D_{L^{-}}^{1-\alpha} \psi_{L^{-}}(u) du, \quad (3.11) $$

where

$$ \psi_{L^{-}}(u) := \psi(u) - \psi(L^{-}). $$

This definition does not depend on $\alpha$, and we have

$$ \int_{0}^{z} \phi d\psi = \int_{0}^{L} \phi \mathbf{1}_{(0,z)} d\psi. $$
Moreover, according to [168], this integral can be extended to more general classes of function thanks to the inequality

$$ \left| \int_{0}^{L} \phi d\psi \right| \leq \|\phi\|_{\alpha,1} \Lambda_{\alpha}(\psi), $$

where

$$ \|\phi\|_{\alpha,1} := \int_{0}^{L} \left( \frac{|\phi(u)|}{u^{\alpha}} + \int_{0}^{u} \frac{|\phi(u) - \phi(v)|}{(u-v)^{\alpha+1}} dv \right) du, $$

and

$$ \Lambda_{\alpha}(\psi) := \frac{1}{\Gamma(1-\alpha)\Gamma(\alpha)} \sup_{0<u<z<L} |D_{z-}^{1-\alpha}\psi_{z-}(u)|. $$

Consequently, this integral is well-defined as soon as $\phi \in W^{\alpha}(0, L, \mathbb{R})$ and $\Lambda_{\alpha}(\psi) < \infty$.

As a result, for a random function $F \in W^{\alpha}(0, L, L^2(\mathbb{R}^2))$, the stochastic integral with respect to the fractional field $W_{\mathfrak{S}}$,

$$ \int_{0}^{z} F(u, \mathbf{x}) dW_{\mathfrak{S}}(u, \mathbf{x}), $$

is defined by (3.11) almost everywhere in $\mathbf{x}$ and $\mathbb{P}$-almost surely. In fact, for all $\alpha \in (1 - \mathfrak{S}, 1/2)$ and $z \in [0, L]$, we have

$$ \left\| \int_{0}^{z} F(u) dW_{\mathfrak{S}}(u) \right\|_{L^{2}(\mathbb{R}^{2})} \leq C \|F\|_{W^{\alpha}(0, L, L^{2}(\mathbb{R}^{2}))} \sup_{\mathbf{x} \in \mathbb{R}^{2}} \Lambda_{\alpha}(W_{H}(\mathbf{x})), $$

knowing that

$$ \mathbb{E} \left[ \sup_{\mathbf{x} \in \mathbb{R}^2} \Lambda_{\alpha}(W_{\mathfrak{S}}(\mathbf{x})) \right] < \infty. $$

**The fractional Itô-Schrödinger equation.** The definition of a solution for the fractional Itô-Schrödinger equation (3.9) is made precise in the definition below, but first let us introduce some additional notations. Let $k \in \mathbb{N}^*$, and $H^k(\mathbb{R}^2)$ be the $k$-th Sobolev space on $\mathbb{R}^2$. Consider also $W_k^\alpha(0, L) := W^\alpha(0, L, H^k(\mathbb{R}^2))$, equipped with the norm $\| \cdot \|_{\alpha, H^k(\mathbb{R}^2)}$, and the complete metric space

$$ W_{\infty}^{\alpha}(0, L) := \bigcap_{k \in \mathbb{N}^{*}} W_{k}^{\alpha}(0, L), $$

equipped with

$$ d_{\alpha, \infty}(\phi, \psi) := \sum_{k \ge 1} \frac{1}{2^k} (1 - \|\phi - \psi\|_{\alpha, H^k(\mathbb{R}^2)}). $$

**Definition 3.2.1** Let $\mathfrak{S} \in (1/2, 1)$ and $\alpha \in (1 - \mathfrak{S}, 1/2)$. We say that $\Psi_\omega \in W_\infty^\alpha(0, L)$ is a pathwise solution of (3.9) if, with probability one, for all $(z, \mathbf{x}) \in [0, L] \times \mathbb{R}^2$, we have

$$ \Psi_{\omega}(z, \mathbf{x}) = \Psi_{\omega}(0, \mathbf{x}) + \frac{i}{2k_{\omega}} \int_{0}^{z} \Delta_{\mathbf{x}} \Psi_{\omega}(u, \mathbf{x}) du + ik_{\omega} \int_{0}^{z} \Psi_{\omega}(u, \mathbf{x}) dW_{\mathfrak{S}}(u, \mathbf{x}). $$

In other words, a solution to (3.9) is a pointwise solution of this equation for almost all realizations. It is worth mentioning that a solution to (3.9) has automatically Hölder regularity

$$ \Psi_{\omega} \in C_{\infty}^{\mathfrak{S}-\theta}(0, L) := \bigcap_{k \in \mathbb{N}^*} C^{\mathfrak{S}-\theta}([0, L], H^k(\mathbb{R}^2)) \text{ with } \theta = \mathfrak{S} + \alpha - 1. $$

Here, $C_{\infty}^{\mathfrak{S}-\theta}(0, L)$ is a complete metric space equipped with

$$ d_{\mathfrak{S}-\theta, \mathcal{C}, \infty}(\phi, \psi) := \sum_{k \ge 1} \frac{1}{2^k} (1 - \|\phi - \psi\|_{\mathfrak{S}-\theta, \mathcal{C}, H^k(\mathbb{R}^2)}), $$

so that $C_{\infty}^{\mathfrak{S}-\theta}(0, L) \subset W_{\infty}^{\alpha}(0, L)$ according to (3.10). In the proof of Thm. 3.2.1, $\Psi_\omega(z=0)$ is assumed to belong to all the $H^k(\mathbb{R}^2)$'s in order to apply a fixed point type argument. This is the reason why intersections with all the $H^k(\mathbb{R}^2)$'s are involved in Def. 3.2.1.

Finally, note that the stochastic integral is here the fractional equivalent to the Itô-Stratonovich integral for a standard Brownian motion, and as such satisfies the classical integration by parts formula. This formally yields the conservation relation: for all $z \in [0, L]$,

$$ \|\Psi_{\omega}(z)\|_{L^2(\mathbb{R}^2)} = \|\Psi_{\omega}(0)\|_{L^2(\mathbb{R}^2)}. $$
## 3.3 Effective fractional paraxial wave equation

In various contexts such as geophysics, laser beam propagation through the atmosphere, or medical imaging for instance, frequency-dependent attenuations have been observed at rates proportional to

$$|\omega|^{\lambda} \quad \lambda \in (0, 2), \qquad (3.12)$$

for a given angular frequency $\omega$ [62, 77, 90, 101, 116, 130, 166, 190]. Depending on the field of application, this type of power-law attenuation can be referred to as anomalous diffusion, nonexponential relaxation, inelastic damping, hysteric damping, singular hereditary, or singular memory media [49]. The development of accurate wave propagation models exhibiting power-law attenuations is therefore of great importance for applications in imaging and inverse problems (see [162] for a survey focusing on medical applications and references therein). To capture such power-law decays, several models have been proposed involving fractional derivatives (see [46, 89, 110, 116, 162, 177, 198, 199] for instance). A particular attention has been paid to space-time models as they enable numerical simulations of a wide range of boundary value problems [109] and can be more straightforward to implement and less costly [210] compared to frequency-domain models.

An alternative approach, which has gained more attention recently, involves randomly perturbed propagation media with fractal correlation structures or long-range dependencies (see [84, 85, 191] and [116, Chap. 9]). This section aims at describing a mathematical derivation, from first principles of physics, of a paraxial wave equation exhibiting a power-law attenuation (3.12) with exponent $\lambda \in (1, 2]$, in the context of random propagation media with long-range correlations. Both aspects, paraxial approximation and random medium fluctuations with long-range correlations, are again considered under the same limit. The case of mixing random fluctuations is also discussed for comparison. The two previous sections, based on non-central limit type results and corresponding to [153, 154] for 1D propagation media, do not provide any power-law attenuation in frequency. To obtain such a property, we rather consider a central limit theorem scaling within the context of long-range correlations, which presents mathematical challenges. This approach has been explored for 1D propagation media [84, 85], but the method used to exhibit a power-law attenuation of the form (3.12) does not seem to apply in a 3D setting. For technical reasons similar to those of Chap. 1, the random fluctuations of the propagation medium are here driven by (1.6). Despite this restrictive assumption, compared to [84, 85], the approach we propose can also be applied to more general 3D settings with non-layered fluctuations, at the cost of additional technical difficulties. Nevertheless, more general fluctuation models with long-range correlations should not change the overall results as the asymptotic equation and scattering coefficients depend only on the correlation function of the fluctuation model, and not its precise definition.

As already pointed out, the scaling associated to the central limit theorem can be seen as propagating the non-central limit scaling over longer propagation distances. This latter scaling, which already leads to an effective phase modulation driven by a fractional Brownian motion, results in waves undergoing very fast oscillations over extended propagation distances, and then have to be treated appropriately to still exhibit effective nontrivial effects. In Chap. 1, we look at the energy of the Hamiltonian system (1.13) to exhibit the diffusion effects. For the random Schrödinger equation the Wigner transform is used to study the energy density propagation by analyzing the correlations of the wave function. The quadratic nature of the Wigner transform w.r.t. the wave function naturally cancels out some rapid phases, and provides an effective description of the energy propagation [92, 93]. For classical wave propagation problems, a similar approach consists in looking at the wave-front along a proper random characteristic time-frame. The rapid phases still produce some effects by averaging the stochasticity to obtain a deterministic spreading for the wave-front [84, 85].

### 3.3.1 The random fluctuations

As for Sect. 3.1 and 3.2, the scalar wave model (3.1) is considered, but here with $D = \mathbb{R}^2$. The main difference w.r.t. the two previous sections lies in the wave-speed profile, which has here the form

$$\frac{1}{c^2(z)} := \frac{1}{c_0^2} \left( 1 + \nu \left( \frac{z}{l_c} \right) \mathbf{1}_{(0,L)}(z) \right), \qquad (3.13)$$
representing a randomly layered propagation medium. For this section, the random fluctuations are assumed to be of the form

$$ \nu(z) := \Theta(\sigma V(z)). $$

Here, $\Theta$ is an odd smooth bounded function with

$$ \theta'_{0} := \Theta'(0) \neq 0 \quad \text{and} \quad \sup |\Theta| < 1, $$

and $V$ is the mean-zero stationary Gaussian random process given by (1.6) satisfying

$$ \mathbb{E}[V(z+z')V(z')] \underset{|z| \to +\infty}{\sim} \frac{R_0}{|z|^\gamma}. \qquad (3.14) $$

The function $\Theta$ plays no significant role here. Gaussian random processes being not bounded, this function just guarantees, for a modelization purpose, that $c^2(z)$ is actually always positive (see (3.13)). For $\sigma$ small enough, the Taylor expansion

$$ \Theta(\sigma V(z)) = \sigma V(z) \theta'_{0} + \mathcal{O}(\sigma^3) $$

($\Theta''(0) = 0$ since $\Theta$ is odd) indicates that the medium fluctuations are driven by $V$.

As in Sect. 1.3.2 of Chap. 1, for the sake of mathematical tractability, a specific form for $V$ is considered. More general fluctuation models have been considered in [153, 154], but the proposed method, based on the rough-path theory, does not seem to apply in the scaling regime described in the next section.

### 3.3.2 Scaling regime

The scaling regime considered in this section corresponds to the one proposed in [84], and is slightly different from the one of Sect. 3.1 and 3.2. The high-frequency regime is here described by

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

and both the correlation length and the central wavelength are of the same order

$$ l_c \sim \lambda, $$

providing a full interaction between the random fluctuations and the propagating waves. For the paraxial approximation to be valid, we assume that the Rayleigh length is of order the propagation distance,

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

Finally, the strength of the fluctuations is assumed to be small to place ourselves in a weak scattering regime,

$$ \sigma \ll 1. $$

To fix the ideas, we set

$$ L \sim 1, \quad \lambda = l_c = \varepsilon^2, \quad \text{and} \quad r_0 = \sigma = \varepsilon, $$

where the choice of $\sigma = \varepsilon$ corresponds to the scaling of the central-limit theorem. This choice allows us to derive a nontrivial limit for both short-range and long-range correlations.

Under this scaling the wave equation (3.1) reads now

$$ \Delta p_{\varepsilon} - \frac{1}{c_{\varepsilon}^2(z)} \partial_{tt}^2 p_{\varepsilon} = f\left(\frac{t}{\varepsilon^2}, \frac{\mathbf{x}}{\varepsilon}\right) \delta(z) \quad (t, \mathbf{x}, z) \in \mathbb{R} \times \mathbb{R}^2 \times \mathbb{R}, $$

with

$$ \frac{1}{c_{\varepsilon}^2(z)} = \frac{1}{c_0^2} \left( 1 + \nu_{\varepsilon} \left( \frac{z}{\varepsilon^2} \right) \mathbf{1}_{(0,L)}(z) \right), \quad \text{and} \quad \nu_{\varepsilon}(z) := \Theta(\varepsilon V(z)). $$
### 3.3.3 Main results

The main result of this section is obtained following the strategy of [84,85]. We introduce the random travel time

$$T_{\varepsilon}(L) := \frac{L}{c_0} + \frac{1}{2c_0} \int_{0}^{L} \nu_{\varepsilon}(z/\varepsilon^2) dz,$$

which corresponds to the expected travel time $L/c_0$ with a random correction, and the wave-front

$$p_{tr,\varepsilon}^{L}(s, \mathbf{y}) := p_{\varepsilon}(T_{\varepsilon}(L) + \varepsilon^2 s, \varepsilon \mathbf{y}, L) \quad (s, \mathbf{y}) \in \mathbb{R} \times \mathbb{R}^2. \quad (3.15)$$

This wave-front corresponds to the wave observed at the end of the random section ($z = L$), on a time window corresponding to the pulse width $\varepsilon^2$, and centered at the random travel time $T_{\varepsilon}(L)$.

Before stating the main result, we introduce some notations. First, the Fourier convention considered in this section is the following:

$$\hat{f}(\omega, \kappa) := \iint f(s, \mathbf{y}) e^{i\omega(s-\kappa \cdot \mathbf{y})} ds \, d\mathbf{y},$$

and

$$f(s, \mathbf{y}) := \frac{1}{(2\pi)^3} \iint \hat{f}(\omega, \kappa) e^{-i\omega(s-\kappa \cdot \mathbf{y})} \omega^2 d\omega \, d\kappa,$$

which is convenient to study space-time problems. Second, we introduce

$$S_0(\mathbb{R} \times \mathbb{R}^2) = \{ \psi \in S(\mathbb{R} \times \mathbb{R}^2) : \int \psi(s, \mathbf{y}) ds = 0, \quad \forall \mathbf{y} \in \mathbb{R}^2 \},$$

where $S(\mathbb{R} \times \mathbb{R}^2)$ stands for the Schwartz class, and $S'_{0,s,\mathbf{y}}(\mathbb{R} \times \mathbb{R}^2)$ the set of tempered distributions restricted to $S_0(\mathbb{R} \times \mathbb{R}^2)$ w.r.t. the variables $s$ and $\mathbf{y}$. This restriction to $S_0(\mathbb{R} \times \mathbb{R}^2)$ is required for the paraxial wave equation (3.18) to be well-posed. One can remark that the source term $f$ belongs to $S_0(\mathbb{R} \times \mathbb{R}^2)$ as soon as its Fourier transform does not support the frequency $\omega = 0$. Below, $C_z^0$ (reps. $C_z^1$) stands for the set of $C^0$-functions (resp. $C^1$-functions) w.r.t. the $z$-variable.

**Theorem 3.3.1** The family $(p_{tr,\varepsilon}^L)_{\varepsilon}$ converges in probability in $C(\mathbb{R} \times \mathbb{R}^2)$, as $\varepsilon \to 0$, to

$$p_{tr}^{L}(s, \mathbf{y}) = \frac{1}{2} \mathcal{K}(\cdot, \cdot, L) * f(s, \mathbf{y}) \quad (s, \mathbf{y}) \in \mathbb{R} \times \mathbb{R}^2, \quad (3.16)$$

where, in the Fourier domain,

$$\hat{\mathcal{K}}(\omega, \kappa, z) := e^{-\theta_0^2 \omega^2 (\Gamma_c(\omega) + i\Gamma_s(\omega))z/(8c_0^2)} e^{-i\omega c_0 |\kappa|^2 z/2}, \quad (3.17)$$

with

$$\Gamma_c(\omega) := 2 \int_0^\infty R(s) \cos \left( \frac{2\omega s}{c_0} \right) ds \quad \text{and} \quad \Gamma_s(\omega) := 2 \int_0^\infty R(s) \sin \left( \frac{2\omega s}{c_0} \right) ds.$$

Here, $R$ is the correlation function of the medium fluctuations (1.7). The convolution kernel $\mathcal{K}$ is the unique solution in $C_z^0([0, \infty), S'_{0,s,\mathbf{y}}(\mathbb{R} \times \mathbb{R}^2)) \cap C_z^1((0, \infty), S'_{0,s,\mathbf{y}}(\mathbb{R} \times \mathbb{R}^2))$ of the paraxial wave equation

$$\partial_{sz}^2 \mathcal{K} - \frac{c_0}{2} \Delta_{\mathbf{y}} \mathcal{K} - \mathcal{I}(\mathcal{K}) = 0, \quad (3.18)$$

with $\mathcal{K}(s, \mathbf{y}, z = 0) = \delta(s)\delta(\mathbf{y})$, and

$$\mathcal{I}(\phi)(s) := \frac{\theta_0^2}{8c_0^2} \int_{-\infty}^{s} R\left(\frac{c_0(s-\tau)}{2}\right) \partial_{sss}^3 \phi(\tau) d\tau \quad s \in \mathbb{R}. \quad (3.19)$$

The asymptotic transmitted wave-front $p_{tr}^L$, at the end of the random section ($z = L$), can be written in term of a convolution where $\mathcal{K}$ represents the pulse deformation. From the explicit formulation of $\mathcal{K}$ in the Fourier domain, the pulse shape is affected in a way which is consistent with the standard ODA theory [70, Chap. 8], even if we are not considering fluctuations with short-range correlations. Typically, according to this theory, the transmitted pulse exhibits a deterministic spreading characterized by a frequency-dependent attenuation and phase modulation. Here, these effects are given
<figure><img src="image_7.png" /><figcaption>Fig. 3.3: Illustrations of the profile $p_{tr}^L(s, y_1, y_2 = 0)$ in the homogeneous case (blue lines) and for two values of $\beta$ ($\beta = 1/2$ for the orange curves and $\beta = 1/6$ for the green curves). We take $\alpha = 1/4$, $\mu = 2$, $a(p) = \mathbf{1}_{(-10,10)}(p)$, $L = 5$, $c_0 = \theta'_0 = 1$ and a source profile given by $\hat{f}(\omega, \kappa) = 2\omega^2 e^{-\omega^2(1+\kappa^2)}$ centered at $s = 0$ in the time domain for simplicity.</figcaption></figure>

through $\omega^2\Gamma_c(\omega)$ (which is positive thanks to the Bochner theorem) and $\omega^2\Gamma_s(\omega)$ respectively. These two terms are well defined even for slowly decaying correlations due to the oscillatory functions. We refer to Fig. 3.3 for illustrations regarding the influence of the kernel $\mathcal{K}$ on the pulse spreading.

The standard ODA theory also provides a random time-shift driven by a standard Brownian motion, meaning that the transmitted pulse exhibits a random arrival time at $z = L$ of order the pulse width. In the context of long-range correlations the situation is more delicate. The aforementioned random time-shift is already compensated in Thm. 3.3.1 by considering the random travel time $T_\epsilon(L)$ in (3.15). Due to the asymptotic behavior of $T_\epsilon(L)$, which is described in the following result, this time-shift allows to remove pathological behaviors when studying the asymptotic of the transmitted wave front.

**Proposition 3.3.1 Let us defined**

$$W_{\epsilon}(L) := \frac{1}{\sigma_{\epsilon}} \left( T_{\epsilon}(L) - \frac{L}{c_{0}} \right),$$

where

$$\sigma_{\epsilon} := \begin{cases} \epsilon^{1+\gamma} & \text{if } \gamma \in (0, 1), \\ \epsilon^2 |\ln(\epsilon)|^{1/2} & \text{if } \gamma = 1, \\ \epsilon^2 & \text{if } \gamma > 1. \end{cases}$$

The family $(W_\epsilon(L))_\epsilon$ converges in law to a limit $W_0(L)$, as $\epsilon \to 0$, where:

* for $\gamma \in (0, 1)$, $W_0$ is a fractional Brownian motion with Hurst index

$$\mathfrak{H} = 1 - \frac{\gamma}{2} \in (1/2, 1),$$

and

$$\mathbb{E}[W_0(L)^2] = L^{2\mathfrak{H}} \frac{\theta_0'^2 R_0}{\mathfrak{H}(2\mathfrak{H} - 1)};$$
• for $\gamma \ge 1$, $W_0$ is a Brownian motion such that

$$ \mathbb{E}[W_0(L)^2] = L \theta_0^{\prime 2} \Gamma_0, \quad \text{with} \quad \Gamma_0 := \begin{cases} \frac{2a(0)}{\mu \beta} & \text{if } \gamma = 1, \\ \Gamma_c(0) & \text{if } \gamma > 1. \end{cases} $$

For rapidly decaying correlations ($\gamma > 1$), $T_\varepsilon(L)$ can be approximated by a Gaussian variable with mean $L/c_0$ and a standard deviation of order the pulse width, which is consistent with the standard ODA theory. In the case of long-range correlations, $T_\varepsilon(L)$ can also be approximated by a Gaussian variable but with a standard deviation very large w.r.t. the pulse width. This is the reason why we compensate this term in (3.15), as it leads to blowing up terms in the derivation of the pulse spreading.

Regarding the backscattered signal at $z = 0$, it can be shown that

$$ p_{bk,\varepsilon}(s, \mathbf{y}) := p_{\varepsilon}(\varepsilon^2 s, \varepsilon \mathbf{y}, 0) $$

converges in probability to 0, in $C(\mathbb{R} \times \mathbb{R}^2)$, as $\varepsilon \to 0$. This is consistent with [70, Chap. 9], in which the authors show that the backscatter wave is made of a small amplitude incoherent signal that can be described through a random field.

In case of long-range correlations, the operator $\mathcal{I}$ can be approximated by a Weyl fractional derivative whose order depends directly on the decay rate $\gamma$ of the correlation function $R$ at infinity. The Weyl fractional derivative is given for $\gamma \in (0, 1)$ by

$$ D^{\gamma} \phi(s) := \frac{\gamma}{\Gamma(1-\gamma)} \int_{-\infty}^{s} \frac{\phi(s) - \phi(\tau)}{(s-\tau)^{1+\gamma}} d\tau \quad s \in \mathbb{R}, $$

whenever this quantity is well-defined, and $\Gamma$ stands for the Gamma function. For instance, $\phi$ can be a bounded $\gamma'$-Hölder function with $\gamma < \gamma'$. However, for $C^1$-functions with fast enough decay at $-\infty$, the Weyl derivative can be rewritten as

$$ D^{\gamma} \phi(s) = \frac{1}{\Gamma(1-\gamma)} \int_{-\infty}^{s} \frac{\phi'(\tau)}{(s-\tau)^{\gamma}} d\tau. $$

To define higher order derivatives, one can just set

$$ D^{j+\gamma} \phi(s) := D^{\gamma} \phi^{(j)}(s) = \frac{1}{\Gamma(1-\gamma)} \int_{-\infty}^{s} \frac{\phi^{(j+1)}(\tau)}{(s-\tau)^{\gamma}} d\tau \quad j \in \mathbb{N}, $$

assuming $\phi$ smooth enough, with enough decay at $-\infty$ of its derivatives $\phi^{(j)}$. These latter requirements hold true for the kernel $\mathcal{K}$ as soon as $z > 0$ thanks to the damping term $\omega^2 \Gamma_c(\omega)$ in (3.17). Therefore, (3.14) and (3.19) suggests that the integro-differential operator $\mathcal{I}$ in (3.18) can be approximated as follows,

$$ \mathcal{I}(\mathcal{K}) \propto D_s^{2+\gamma} \mathcal{K}. $$

Below, we emphasize that the fractional derivative $D^{2+\gamma}$ acts on the $s$-variable with the notation $D_s^{2+\gamma}$. This observation can be derived by considering the scaled correlation function

$$ \sigma(l_0) R(z/l_0), $$

where $l_0$ will be sent to 0, and

$$ \sigma(l_0) := \begin{cases} \frac{1}{l_0^{\gamma}} & \text{if } \gamma \in (0, 1), \\ \frac{1}{l_0 |\ln(l_0)|} & \text{if } \gamma = 1, \\ \frac{1}{l_0} & \text{if } \gamma > 1. \end{cases} $$

In other words, the correlation length $l_0$ is assumed to be small compared to the pulse duration, and we can defined the transmitted wave-front $p_{tr,l_0}^L$ from (3.16) accordingly.
**Theorem 3.3.2** The family $(p_{tr,l_0}^L)_{l_0}$ converges in $C(\mathbb{R} \times \mathbb{R}^2)$, as $l_0 \to 0$, to

$$p_{tr,0}^{L}(s, \mathbf{y}) := \frac{1}{2} \mathcal{K}_{0}(\cdot, \cdot, L) * \Psi(s, \mathbf{y}) \quad (s, \mathbf{y}) \in \mathbb{R} \times \mathbb{R}^{2}.$$

Here, $\mathcal{K}_0$ is defined in the Fourier domain by

$$\hat{\mathcal{K}}_0(\omega, \kappa, z) := \begin{cases} e^{-\theta_0^{\prime 2} \omega^2 \Gamma_0 z / (8c_0^2)} e^{-i\omega c_0 |\kappa|^2 z / 2} & \text{if } \gamma \ge 1, \\ e^{-\theta_0^{\prime 2} R_0 |\omega|^{1+\gamma} (\Gamma_{c,0}(\omega) + i\Gamma_{s,0}(\omega))z / (8c_0^2)} e^{-i\omega c_0 |\kappa|^2 z / 2} & \text{if } \gamma \in (0, 1), \end{cases} \quad (3.20)$$

with

$$\Gamma_{c,0}(\omega) = \Gamma(1 - \gamma) \cos\left(\frac{(1 - \gamma)\pi}{2}\right) \left(\frac{2}{c_0}\right)^{\gamma - 1},$$

$$\Gamma_{s,0}(\omega) = \Gamma(1 - \gamma) \sin\left(\frac{(1 - \gamma)\pi}{2}\right) \left(\frac{2}{c_0}\right)^{\gamma - 1} \text{sign}(\omega),$$

and

$$\Gamma_0 := \begin{cases} \frac{2a(0)}{\mu\beta} & \text{if } \gamma = 1, \\ \Gamma_c(0) & \text{if } \gamma > 1. \end{cases}$$

Moreover, $\mathcal{K}_0$ is the unique solution in $C_z^0([0, \infty), S_{0,s,\mathbf{y}}'(\mathbb{R} \times \mathbb{R}^2)) \cap C_z^1((0, \infty), S_{0,s,\mathbf{y}}'(\mathbb{R} \times \mathbb{R}^2))$ of the paraxial wave equation

$$\partial_{sz}^2 \mathcal{K}_0 - \frac{c_0}{2} \Delta_{\mathbf{y}} \mathcal{K}_0 - \mathcal{I}_0(\mathcal{K}_0) = 0,$$

with $\mathcal{K}_0(s, \mathbf{y}, z = 0) = \delta(s)\delta(\mathbf{y})$, and

$$\mathcal{I}_0(\phi) := \begin{cases} \frac{\theta_0^{\prime 2} R_0}{8c_0^2} \partial_{sss}^3 \phi & \text{if } \gamma \ge 1, \\ \frac{\theta_0^{\prime 2} R_0 \Gamma(1-\gamma)}{2^{3-\gamma} c_0^{1+\gamma}} D_s^{2+\gamma} \phi & \text{if } \gamma \in (0, 1). \end{cases}$$

The difference between the case $\gamma \ge 1$ and $\gamma \in (0, 1)$ can easily be highlighted from this result. In the former case, $\mathcal{I}_0$ is a classical third order differential operator, while for $\gamma \in (0, 1)$, it is a fractional derivative of order $2 + \gamma \in (2, 3)$. Moreover, even if the case $\gamma = 1$ corresponds to slowly decaying correlations, the kernel $\mathcal{K}_0$ behaves as for short-range correlations $\gamma > 1$. The case $\gamma = 1$ plays somehow the role of a continuity point w.r.t. the order of derivation in $\mathcal{I}_0$. For long-range correlations, one can observe in (3.20) the frequency-dependent attenuation given by the power law

$$|\omega|^{1+\gamma} \quad \gamma \in (0, 1],$$

with exponent depending directly on the decay rate of the correlation function of the medium fluctuations (3.14).

## 3.4 Perspectives

Wave propagation in randomly perturbed media with long-range correlations exhibits interesting multiscale phenomena that differs from the standard ODA theory. While being mathematically challenging, the main perspective in this direction concerns the derivation of effective diffusion effects with random fluctuations involving the lateral variable. Although the context of a waveguide does not bring much more difficulty compared to the last section of this chapter, extending this type of result in the context of Sect. 3.2 appears to be a challenging problem. A second aspect that needs to be investigated, as in Chap. 1 and 2, concerns the restriction over the noise model. For instance, in
the context of Sect. 3.3, we are not able to capture power-law attenuations with exponent between 0 and 1. Such a case can be captured under short-range correlations with the additional anti-persistence assumption [84],

$$ \int_{-\infty}^{\infty} R(z) dz = 0, $$

but our noise model (1.6) is not compatible with this assumption.

A deeper analysis of the fractional Itô-Schrödinger equation could be of interest for practical uses. In fact, due to the nature of the stochastic integral, it is not straightforward to provide a simple formulation for basic quantities such as the coherent wave (the expectation of the solution of the fractional Itô-Schrödinger equation). Even if our approach provides an explicit expression for all the order moments, the resulting formulas are quite complex. This is in contrast with the standard Itô-Schrödinger equation involving a classical Itô-Stratonovich integral, and the resulting martingale property, for which the coherent wave can easily be evaluated. As a perspective, our fractional Itô-Schrödinger equation could be investigated under the so-called scintillation or spot dancing regimes [58,86]. These asymptotics could lead to simplified expression of the resulting moments and quantities of interest.
Chapter 4

Radiative transfer with nonintegrable
singular scattering kernels

Radiative transfer models have been used for more than a century to describe wave energy propagation
through complex/random media [48, 124], as well as neutron transport [142, 143, 192], heat transfer
[207], and are still active research subjects in astrophysics, geophysics, and optical tomography for
instance [150, 167, 175, 176].

This chapter discusses interesting effects and numerical difficulties for radiative transfer models
involving a singular scattering kernel:

$$
\begin{equation}
\begin{cases}
\partial_t u + \hat{\mathbf{k}} \cdot \nabla_x u = \mathcal{Q} u \\
u(t=0, \mathbf{x}, \hat{\mathbf{k}}) = u_0(\mathbf{x}, \hat{\mathbf{k}})
\end{cases}
\qquad (t, \mathbf{x}, \hat{\mathbf{k}}) \in (0, \infty) \times \mathbb{R}^d \times \mathbb{S}^{d-1}, \tag{4.1}
\end{equation}
$$

where $\mathbb{S}^{d-1}$ denotes the unit sphere in $\mathbb{R}^d$, and $u$ is the wave energy density in the context of wave propagation or the particle distribution function in the context of neutronics. The scattering operator $\mathcal{Q}$ has the standard form

$$
(\mathcal{Q}u)(\hat{\mathbf{k}}) = \int_{\mathbb{S}^{d-1}} \Phi(|\hat{\mathbf{p}} - \hat{\mathbf{k}}|) (u(\hat{\mathbf{p}}) - u(\hat{\mathbf{k}})) \sigma(d\hat{\mathbf{p}}),
$$

where $\sigma(d\hat{\mathbf{p}})$ stands for the surface measure on $\mathbb{S}^{d-1}$. Here, the scattering kernel $\Phi$ is assumed to be
of the form

$$
\Phi(|\hat{\mathbf{p}} - \hat{\mathbf{k}}|) := \frac{a(|\hat{\mathbf{p}} - \hat{\mathbf{k}}|)}{|\hat{\mathbf{p}} - \hat{\mathbf{k}}|^{d-1+\alpha}} = F(\hat{\mathbf{k}} \cdot \hat{\mathbf{p}}),
$$

with

$$
F(s) := \frac{a(\sqrt{2(1-s)})}{(2(1-s))^{(d-1+\alpha)/2}}, \quad s \in [-1, 1], \quad \alpha \in (0, 2),
$$

and $a$ is a nonnegative function. The difficulty in this situation resides in the vanishing of the mean free time $t_0$ associated with $\mathcal{Q}$:

$$
\frac{1}{t_0} = \int_{\mathbb{S}^{d-1}} \Phi(|\hat{\mathbf{k}} - \hat{\mathbf{p}}|) \sigma(d\hat{\mathbf{p}}) = \infty. \quad (4.2)
$$

In the context of gas dynamics and the Boltzmann equation, a singular collision operator refers to
the so-called non cut-off case [3,206], where the particle interactions have long-range properties. Such
a scenario also arises in the context of highly peaked-forward light scattering in biological tissues, in
turbulent atmosphere, or more generally in the context of wave propagation in random media with
long-range correlations.

Here are some examples for $d=3$. When the momentum variable is restricted to the unit sphere,
radiative transfer equations (RTEs) can be derived from high-frequency wave propagation in random