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

HAL
open science

# Stochastic Dynamics of Incoherent Branched Flow

Josselin Garnier, Antonio Picozzi, Theo Torres

► To cite this version:

Josselin Garnier, Antonio Picozzi, Theo Torres. Stochastic Dynamics of Incoherent Branched Flow. Physical Review Letters, 2025, 134, pp.223803. 10.48550/arXiv.2502.07028 . hal-05053067

HAL Id: hal-05053067

https://hal.science/hal-05053067v1

Submitted on 5 May 2025

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
Stochastic Dynamics of Incoherent Branched Flow

Josselin Garnier¹, Antonio Picozzi², and Theo Torres²

¹ CMAP, CNRS, Ecole polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France and

² Laboratoire Interdisciplinaire Carnot de Bourgogne, CNRS, Université Bourgogne Europe, Dijon, France

Waves propagating through weakly disordered smooth linear media undergo a universal phenomenon called branched flow. Branched flow has been observed and studied experimentally in various systems by considering coherent waves. Recent experiments have reported the observation of optical branched flow by using an incoherent light source, thus revealing the key role of coherent phase-sensitive effects in the development of incoherent branched flow. By considering the paraxial wave equation as a generic representative model, we elaborate a stochastic theory of both coherent and incoherent branched flow. We derive closed-form equations that determine the evolution of the intensity correlation function, as well as the value and the propagation distance of the maximum of the scintillation index, which characterize the dynamical formation of incoherent branched flow. We report accurate numerical simulations that are found in quantitative agreement with the theory without free parameters. Our theory highlights the important impact of coherence and interference on branched flow, thereby providing a framework for exploring branched flow in nonlinear media, in relation with the formation of freak waves in oceans.

PACS numbers: 42.25.Dd, 05.40.-a.

**Introduction.** Waves passing through a weakly disordered smooth medium with a correlation radius larger than the wavelength form long, narrow filaments called branches. Instead of random speckle patterns, the disordered potential focuses waves into branches that split, creating a tree-like structure, known as branched flow (BF). Originally observed in electrons [1–5] and microwave cavities [6, 7], BFs have then been anticipated to occur with vastly different wavelength scales [8]. They may serve as a catalyst for the emergence of extreme nonlinear events [9–13], and freak waves on the ocean [14–19]. BFs have also been suggested to occur for sound waves [20], ultra-relativistic electrons in graphene [21], flexural waves in elastic plates [22], while they can act as a conduit for energy transmission in scattering media [23]. BFs have been extended to random potentials in space and time [24], to periodic potentials [25, 26], and even to active random walks [27]. More recently, BFs have been observed experimentally with optical waves propagating in soap films [28], and the control of light BF through weakly disordered media has become an important challenge [23, 29–31].

The formation of BFs can be explained using geometrical optics, where local maxima of the random refraction index act as lenses, creating caustics and high wave intensities, as originally described in [32, 33]. Numerical simulations show that the scintillation index (i.e., the relative variance of the intensity fluctuations) can exceed one in such cases. Recent studies have used geometrical optics or diffraction integrals in the framework of catastrophe optics [34, 35] to derive the scaling behavior of BF dynamics [36, 37] and extreme waves [38, 39]. Actually, except for some particular theoretical studies [40, 41], BFs have been essentially treated in the framework of ray caustics, then disregarding coherence or interference effects [8]. Along this way, experiments have been carried out essentially with coherent waves, such as coherent electron waves [1], coherent microwave [6, 7], or with coherent laser light [28]. On the other hand, in recent experiments, optical BFs have been studied by using incoherent light sources [42], revealing intriguing properties

about the role of coherence in the formation and the evolution of BFs, such as coherent interference between the different wave fronts and the sensitivity of BFs to the coherence of the waves.

Our aim in this Letter is to elaborate a stochastic formulation of BFs by considering an initial random wave function propagating in a random potential. Using the paraxial wave (Schrödinger) equation as a representative model, we show that interference effects deeply modify the statistical properties of BFs. Employing multiscale and stochastic calculus, we derive closed-form equations that give the evolution of the intensity correlation function. In particular we describe the evolution of the scintillation index that characterizes the dynamical formation of incoherent BFs. We determine that the scintillation index is a function of two dimensionless parameters that we identify and that involve the statistics of the medium and of the initial field. The theory of the stochastic dynamics of BFs is validated by accurate numerical simulations, which are found in quantitative agreement with the theory, without using any adjustable parameter.

**Model.** We consider the two-dimensional paraxial wave equation [43, 44]:

$$i \partial_z \psi_z = -\alpha \partial_x^2 \psi_z + V(z, x) \psi_z, \quad z > 0, x \in \mathbb{R}, \quad (1)$$

starting from $\psi_{z=0}(x) = \psi_o(x)$, where $\psi_o$ is a coherent or partially coherent field and $V$ is a smooth and slowly varying potential, which we assume to be a random process. We will denote by $\mathbb{E}[\cdot]$ the expectation with respect to the distribution of this random process.

We present our work in optics as a concrete example, but the paraxial wave Eq.(1) is widespread in physics, making the processes discussed herein broadly applicable to various systems. In optics, the parameter $\alpha$ and the potential $V$ are related to the index of refraction $n$ as follows: $\alpha = 1/(2k_o n_o)$, $V(z, x) = k_o(n_o^2 - n^2(z, x))/(2n_o)$, where $k_o$ is the wavenumber in free space, $n_o$ is the homogeneous background index of refraction, and $n(z, x)$ is the spatially dependent index of refraction of the medium.
We will consider two different types of initial field.

1. We will first consider the coherent case in which the initial field is a plane wave: $\psi_o(x) = 1$. The measured intensity is $|\psi_z(x)|^2$, the mean intensity is $\mathbb{E}[|\psi_z(x)|^2]$, and the scintillation index (i.e., the relative variance of the intensity) is

$$S_z(x) = \frac{\mathbb{E}[|\psi_z(x)|^4] - \mathbb{E}[|\psi_z(x)|^2]^2}{\mathbb{E}[|\psi_z(x)|^2]^2}. \quad (2)$$

2. We will then consider in detail the situation in which the initial field is a coherent or partially coherent speckled field. We will consider the two following situations:

(c) $\psi_o$ is a coherent speckled field, which will be modeled as a stationary random field with Gaussian statistics and correlation radius $\rho_o$ (the width of the field correlation function).

(pc) $\psi_o$ is a partially coherent speckled field, which will be modeled as a time-dependent random field with Gaussian statistics and correlation radius $\rho_o$.

Such fields can be generated by passing a time-harmonic plane wave through a static (case (c)) or rotating (case (pc)) diffuser, which features a random arrangement of scattering centers. This experimental setup has attracted considerable interest due to its ability to mimic the properties of a thermal light source [45, 46], with the added advantage of controlling the spatial and temporal coherence properties from the degree of roughness of the diffuser and its rotation speed. It has been used to investigate speckle phenomena [46], ghost imaging [47, 48] and incoherent BFs [42].

We will denote by $\langle \cdot \rangle$ the expectation with respect to the distribution of the initial field. In situation (c), the measured intensity is $|\psi_z(x)|^2$, the mean intensity is $\mathbb{E}[\langle |\psi_z(x)|^2 \rangle]$ and the scintillation index is

$$S_z^{(c)}(x) = \frac{\mathbb{E}[\langle |\psi_z(x)|^4 \rangle] - \mathbb{E}[\langle |\psi_z(x)|^2 \rangle]^2}{\mathbb{E}[\langle |\psi_z(x)|^2 \rangle]^2}. \quad (3)$$

In situation (pc), assuming that the response time of the photodetector is larger than the coherence time of the field, the measured intensity is $\langle |\psi_z(x)|^2 \rangle$ (the averaging $\langle \cdot \rangle$ is experimentally carried out by time averaging by the detector over the multiple initial conditions generated by the rotating diffuser), the mean intensity is $\mathbb{E}[\langle |\psi_z(x)|^2 \rangle]$, and the scintillation index is

$$S_z^{(\text{pc})}(x) = \frac{\mathbb{E}[\langle |\psi_z(x)|^2 \rangle^2] - \mathbb{E}[\langle |\psi_z(x)|^2 \rangle]^2}{\mathbb{E}[\langle |\psi_z(x)|^2 \rangle]^2}. \quad (4)$$

*Coherent initial plane wave.* In this paragraph we assume a regime in which: i) the wavelength $\lambda = 2\pi/k_o$ is much smaller than the correlation radius $\ell_c$ of the index of refraction of the medium, ii) the variance $\sigma_n^2$ of the index of refraction is small (hence, the variance $\sigma^2 = 4\pi^2\sigma_n^2/\lambda^2$ of the random potential satisfies $\sigma^2 \ll 1/\lambda^2$); iii) the propagation distance is large enough so that the evolution of the variance of the intensity is of order one.

This situation has been intensively studied [43, 49]. By a multiscale analysis closed-form equations can be derived for the field and intensity correlation functions [50, 51]. These equations depend on the medium statistics via the integrated

medium correlation function $\gamma$ defined by

$$\gamma(x) = \int_{\mathbb{R}} \mathbb{E}[V(0,0)V(z,x)]dz, \quad (5)$$

which can be written in the form: $\gamma(x) = \sigma^2 \ell_c \tilde{\gamma}(x/\ell_c)$. As a particular example of smooth random medium, we can consider a medium with Gaussian correlation function $\mathbb{E}[V(0,0)V(z,x)] = \sigma^2 \exp(-\ell_c^2/(x^2+z^2))$, so that $\tilde{\gamma}(\tilde{x}) = \sqrt{\pi} \exp(-\tilde{x}^2)$. The field correlation function is $\mathbb{E}[\psi_z(x + \frac{y}{2})\psi_z(x - \frac{y}{2})] = \exp[z(\gamma(y) - \gamma(0))]$ [50]. This shows that the mean intensity is constant in $z$ and $x$ and that the correlation radius of the field decays as $1/\sqrt{z}$ [52].

We introduce two relevant parameters that will play a key role: $X_c = \sigma^{2/3}\ell_c/\alpha^{1/3}$ (which is dimensionless) and $z_c = \ell_c/(2\sigma^{2/3}\alpha^{2/3})$ (which is homogeneous to a length), that will be shown to correspond to the propagation distance at which the scintillation index reaches a maximum for large values of $X_c$. They can also be expressed as $X_c = (2^{4/3}\pi n_o^{1/3})\sigma_n^{2/3}\ell_c/\lambda$ and $z_c = (2^{-1/3}n_o^{2/3})\ell_c/\sigma_n^{2/3}$. From [50] we find that the scintillation index does not depend on $x$:

$$S_z = \tilde{D}_{z/z_c}(0,0) - 1, \quad (6)$$

where $\tilde{D}_{\tilde{z}}(\tilde{x}, \tilde{y})$ satisfies

$$\partial_{\tilde{z}} \tilde{D}_{\tilde{z}} = i X_c^{-1} \partial_{\tilde{x}\tilde{y}}^2 \tilde{D}_{\tilde{z}} + \frac{1}{2} X_c^2 \tilde{U}(\tilde{x}, \tilde{y}) \tilde{D}_{\tilde{z}}, \quad (7)$$

starting from $\tilde{D}_{\tilde{z}=0}(\tilde{x}, \tilde{y}) = 1$, with $\tilde{U}(\tilde{x}, \tilde{y}) = 2\tilde{\gamma}(\tilde{x}) + 2\tilde{\gamma}(\tilde{y}) - \tilde{\gamma}(\tilde{x} + \tilde{y}) - \tilde{\gamma}(\tilde{x} - \tilde{y}) - 2\tilde{\gamma}(0)$ [here $\tilde{x} = x/\ell_c$ and $\tilde{y} = y/\ell_c$]. This shows that the scintillation index is a function of $\tilde{z} = z/z_c$ and $X_c$ only. Eq. (7) can be solved by the split-step Fourier method [53]. Moreover, by expanding the solution for small $\tilde{z}$, we get $S_z \approx [\partial_x^4 \tilde{\gamma}(0)/6](z/z_c)^3$ at leading order [with $\partial_x^4 \tilde{\gamma}(0) = 12\sqrt{\pi}$ for a medium with Gaussian correlation].

Here are the main results.

Firstly, the scintillation index $S_z$ is close to 0 when $z$ is small (i.e., smaller than $z_c$) and it first increases cubically with $z$. Secondly, when $X_c$ is below a threshold value $X_c^{(t)}$ ($X_c^{(t)}$ is between one and three for a medium with Gaussian correlation function, see Fig. 1(b)), the scintillation index is monotonously increasing towards its limit value 1 when $z \to +\infty$. Thirdly, when $X_c$ is above the threshold value $X_c^{(t)}$, the scintillation index reaches a maximal value larger than one at finite propagation distance. The maximal value $\max_z S_z$ depends only on $X_c$, but the distance $z$ at which the maximum of the scintillation index is reached depends on $X_c$ and $z_c$ (it was predicted to be proportional to $\ell_c/\sigma_n^{2/3}$ in previous works [32, 33]). It then relaxes to its limit value 1 when $z \to +\infty$, where the wavefield acquires Gaussian statistics for very large propagation distances [38, 51, 54–56]. Finally, it is quite surprising to note that when $X_c$ is larger than $X_c^{(t)}$, then the scintillation index may present two maxima, one global and one local [see the small bump around $z/z_c = 2.5$ in Fig. 1(b)].

We remark that the distance $\ell_c/\sigma_n^{2/3}$ is also the typical spatial scale of evolution of the number of branches that a Hamiltonian flow develops in a random potential [37]. As
<figure><img src="image_3.png" /><figcaption>FIG. 1: Coherent initial plane wave with a medium with Gaussian correlation: (a) Numerical simulation of Eq.(1) showing the evolution of $|\psi_z(x)|^2$ starting from $\psi_o(x) = 1$. Parameters: $\ell_c/\lambda = 100$, $\sigma^2\lambda^2 = 10^{-4}$ ($X_c \approx 12.4$). (b) Scintillation index $S_z$ versus $z/z_c$ for different values of $X_c$: the black dashed lines report the theory, Eq.(6); the dotted line is the small $z$ prediction $S_z \approx 2\sqrt{\pi}(z/z_c)^3$; the colored lines are the results of the numerical simulations, averaged over 1000 independent realizations of the disordered potential. Parameters: from the bottom, $\ell_c/\lambda = 10, 25, 50, 75$, with $\sigma^2\lambda^2 = 10^{-4}$ for all curves, except for the top yellow curve ($X_c \approx 12.4$) where $\ell_c/\lambda = 50$, $\sigma^2\lambda^2 = 8 \times 10^{-4}$.</figcaption></figure>

discussed in [57, App. A], a ray theory of branched flow can predict the numbers and positions of local intensity maxima, however it cannot predict the values of the maxima that result from interference effects and that depend on the coherence properties of the initial field (see also [41]).

We have also computed the intensity correlation function

$$C_z^I(x) = \frac{\mathbb{E}[|\psi_z(y + \frac{x}{2})|^2 |\psi_z(y - \frac{x}{2})|^2] - \mathbb{E}[|\psi_z(y)|^2]^2}{\mathbb{E}[|\psi_z(y)|^2]^2}, \quad (8)$$

which is independent of $y$ and is given by

$$C_z^I(x) = \tilde{D}_{z/z_c}(x/\ell_c, 0) - 1, \quad (9)$$

where $\tilde{D}$ is solution of Eq.(7). The intensity correlation function is plotted in Fig. 2. Of course one has $C_{z=0}^I(x) = 0$, $C_z^I(0) = S_z$, $C_z^I(x) \to 0$ as $x \to +\infty$, and $\int C_z^I(x)dx = 0$ (this can be interpreted as an energy conservation relation).

We have tested the validity of the theoretical predictions by direct numerical simulations of the paraxial wave Eq.(1) (see [57, App. H]). The results for the evolution of the scintillation index in Fig. 1 and the intensity correlation function in Fig. 2 show excellent quantitative agreements, even though the separation of scales is not strong in the simulations.

**Incoherent initial wave: Scaling regime.** From now on we address the situation in which the initial field $\psi_o(x)$ is a speckled field. In addition to the assumptions (i) to (iii) considered above for the initial plane-wave case, we assume that $\psi_o(x)$ has Gaussian statistics, with correlation radius $\rho_o$ larger than the wavelength $\lambda$ and smaller than the correlation radius $\ell_c$ of the index of refraction, $\ell_c/\lambda \gg \rho_o/\lambda \gg 1$.

We carry out a multiscale analysis in which a dimensionless scale parameter $\varepsilon$ encapsulates the four assumptions listed above. Accordingly, we denote by $\varepsilon \sim \lambda/\ell_c$, the order of magnitude of the ratio of the wavelength over the correlation radius of the index of refraction. We assume that the typical

<figure><img src="image_4.png" /><figcaption>FIG. 2: Coherent initial plane wave with a medium with Gaussian correlation: (a) Theoretical intensity correlation function $C_z^I(x)$ from Eq.(9). (b) Comparison of $C_z^I(x)$ from Eq.(9) (black dashed lines), with the numerical simulations of Eq.(1) (colored lines), for different propagation lengths $z/z_c$. An average over 1000 simulations with different realizations of the random potential $V(z,x)$ has been carried out. Parameters: $\ell_c/\lambda = 100$, $\sigma^2\lambda^2 = 10^{-4}$ ($X_c \approx 12.4$).</figcaption></figure>

amplitude of the fluctuations of the index of refraction is $\varepsilon^c$, with $c > 0$. If we consider that the reference length is the correlation radius of the index of refraction, we can write $\alpha^c = 1/(2k_o n_o) = \varepsilon\alpha$ and $V^c = k_o(n_o^2 - n^2)/(2n_o) = \varepsilon^{c-1}V$, and the scaled paraxial wave equation has the form

$$i \partial_z \psi_z^c = - \varepsilon \alpha \partial_x^2 \psi_z^c + \varepsilon^{c-1} V(z, x) \psi_z^c, \quad z > 0, x \in \mathbb{R}, \quad (10)$$

starting from $\psi_{z=0}^c(x) = \psi_o^c(x)$. The initial field $\psi_o^c$ has a correlation radius of the order of $\varepsilon^d$ (relative to the correlation radius of the index of refraction) for some $d \in (0, 1)$, which means that it is larger than the wavelength (because $d < 1$) and smaller than the correlation radius of the index of refraction (because $d > 0$). The correlation function of the initial field is, therefore, of the form

$$\langle \psi_o^c(x + \varepsilon^d \frac{y}{2}) \bar{\psi}_o^c(x - \varepsilon^d \frac{y}{2}) \rangle = C_o(y). \quad (11)$$

For the numerical simulations we consider the Gaussian model [58, 59] in which $C_o(y) = \exp(-y^2/(4\rho_o^2))$. Finally, we consider the Wigner transform after a propagation distance of order $\varepsilon^{-b}$ (relative to the correlation radius of the index of refraction):

$$W_z^c(x, k) = \int_{\mathbb{R}} \langle \psi_{z=0}^c(x + \varepsilon^d \frac{y}{2}) \bar{\psi}_{z=0}^c(x - \varepsilon^d \frac{y}{2}) \rangle e^{-iky} dy. \quad (12)$$

In the scaling regime $d \in (1/5, 1)$, $b = 1 - d$, $c = 3(1 - d)/2$, we get from (10) (see [57, App. B]) that it satisfies the scaled Vlasov-type equation

$$\partial_z W_z^c + \partial_k \omega_k \partial_x W_z^c - \frac{1}{\varepsilon^{b/2}} \partial_x V(\frac{z}{\varepsilon^b}, x) \partial_k W_z^c = 0, \quad (13)$$

with the initial condition $W_{z=0}^c(x, k) = W_o(k) = \int_{\mathbb{R}} C_o(y) e^{-iky} dy$ and with $\omega_k = \alpha k^2$ (see [57, App. C]). Note that the scaling of the potential in (13) is appropriate for the use of limit theorems for random differential equations [60, Chapter 6] and we will carry out such a multiscale analysis.

Before going to the multiscale analysis, we remark that the solution of the Vlasov equation (13) can be
<figure><img src="image_5.png" /><figcaption>FIG. 3: Incoherent initial wave with a medium with Gaussian correlation: (a) Numerical simulation of Eq.(1) showing the evolution of $|\psi_z(x)|^2$ starting from a coherent speckle field [situation (c)], with $\rho_o/\lambda = 10$, $\ell_c/\lambda = 100$, $\sigma^2\lambda^2 = 10^{-4}$. (b-d) Evolution of $S_z^{(c)}$ versus $z/z_c$, by varying different parameters: the black dashed lines report the theory, Eq.(16); the colored lines are the results of the numerical simulations, averaged over 1000 independent realizations of the disordered potential and of the initial random field. Parameters: (b) $\rho_o/\lambda = 10$, $\sigma^2\lambda^2 = 10^{-4}$; (c) $\ell_c/\lambda = 100$, $\sigma^2\lambda^2 = 10^{-4}$; (d) $\rho_o/\lambda = 10$, $\ell_c/\lambda = 100$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 4: Incoherent initial wave with a medium with Gaussian correlation: Theoretical intensity correlation function $C_z^I(x)$: from Eq.(18) for a coherent speckle field (a) [situation (c)], and from Eq.(19) for a partially coherent speckle field (b) [situation (pc)]. Corresponding comparison of the theoretical correlation function $C_z^I(x)$ (black dashed lines), with the numerical simulations of Eq.(1) (colored lines), for different propagation lengths $z/z_c$. In situation (c), an average over 1000 realizations of $V(z, x)$ and of $\psi_o(x)$, has been considered. In situation (pc) an average over 300 realizations of $V(z, x)$, each with 400 realizations of $\psi_o(x)$. Parameters: $\ell_c/\lambda = 100$, $\rho_o/\lambda = 10$, $\sigma^2\lambda^2 = 10^{-4}$.</figcaption></figure>

expressed in terms of the solutions of random ordinary differential equations. Indeed, using the characteristic method, we have $W_z^\varepsilon(X_z^\varepsilon(x, k), K_z^\varepsilon(x, k)) = W_o(k)$, where $(X_z^\varepsilon(x, k), K_z^\varepsilon(x, k))$ satisfies the ray equations

$$ \frac{dX_z^\varepsilon}{dz} = 2\alpha K_z^\varepsilon, \quad \frac{dK_z^\varepsilon}{dz} = -\frac{1}{\varepsilon^{b/2}} \partial_x V\left(\frac{z}{\varepsilon^b}, X_z^\varepsilon\right), \quad (14) $$

starting from $X_{z=0}^\varepsilon(x, k) = x$, $K_{z=0}^\varepsilon(x, k) = k$. A key result (proved in [57, App. C]) that makes it possible to study the Wigner transform is the following one: For any $X, K \in \mathbb{R}$,

$$ W_z^\varepsilon(X, K) = \int_{\mathbb{R}^2} W_o(k) \delta(X_z^\varepsilon(x, k) - X) \delta(K_z^\varepsilon(x, k) - K) dx dk. \quad (15) $$

By taking an expectation (with respect to the distribution of the random medium), one can see that the mean Wigner transform involves the probability density function (pdf) of $(X_z^\varepsilon(x, k), K_z^\varepsilon(x, k))$. Higher-order moments of the Wigner transform involve multivariate pdf. Those pdf are computed in [57, App. D], and they give the following results.

**Mean Wigner transform.** From (15) we get the expression of the mean Wigner transform in the regime $\varepsilon \to 0$, which in turn gives the expression of the field correlation function $E[(\psi_z(x + \frac{y}{2})\bar{\psi}_z(x - \frac{y}{2}))] = C_o(y) \exp(-\gamma_2 zy^2/2)$, where $\gamma_2 = -\partial_x^2\gamma(0)$. This shows that the mean intensity is constant in $z$ and $x$ and that the correlation radius of the beam decays as $1/\sqrt{z}$ just as in the case of an initial coherent plane wave.

**Scintillation index.** We write the correlation function of the initial field in the dimensionless form $C_o(y) =$

$\tilde{C}_o(y/\rho_o)$, where $\rho_o$ is the correlation radius of the initial field. We introduce the relevant dimensionless parameter $X_o = \sigma^{2/3}\rho_o/\alpha^{1/3}$. We get that in the situation (pc) and (c) the scintillation index does not depend on $x$ (see [57, App. E]):

$$ S_z^{(\text{pc})} = \tilde{\Pi}_{z/z_c}(0,0) - 1, \quad S_z^{(c)} = 2\tilde{\Pi}_{z/z_c}(0,0) - 1, \quad (16) $$

where $\tilde{\Pi}_z(\tilde{x}, \tilde{y})$ is the solution to

$$ \partial_{\tilde{z}} \tilde{\Pi}_{\tilde{z}} = i \partial_{\tilde{x}} \partial_{\tilde{y}} \tilde{\Pi}_{\tilde{z}} - \frac{1}{2} (\tilde{\Gamma}(0) - \tilde{\Gamma}(\tilde{x})) \tilde{y}^2 \tilde{\Pi}_{\tilde{z}}, \quad (17) $$

starting from $\tilde{\Pi}_{\tilde{z}=0}(\tilde{x}, \tilde{y}) = \tilde{\pi}_o(\tilde{y}/X_o)$. Here $\tilde{\Gamma}(\tilde{x}) = -\partial_{\tilde{x}}^2\tilde{\gamma}(\tilde{x})$, $\tilde{\pi}_o(\tilde{y}) = |\tilde{C}_o(\tilde{y})|^2/\tilde{C}_o(0)^2$, while Eq.(17) can be solved by a split-step Fourier method [53]. By expanding the solution of (17) for small $\tilde{z}$, we get $S_{\tilde{z}}^{(\text{pc})} \simeq [\partial_{\tilde{x}}^4\tilde{\gamma}(0)/6](z/z_c)^3$, and $S_{\tilde{z}}^{(c)} \simeq 1 + [\partial_{\tilde{x}}^4\tilde{\gamma}(0)/3](z/z_c)^3$ at leading order, see [57, App. F] [with $\partial_{\tilde{x}}^4\tilde{\gamma}(0) = 12\sqrt{\pi}$ for a medium with Gaussian correlation]. This shows that the early dynamics of the scintillation index in situation (pc) does not depend on the correlation radius $\rho_o$ nor on the correlation function of the initial field, and is equivalent to the behavior valid for an initial plane wave. The scintillation index first grows cubically and then reaches a maximum value, which depends on $X_o = \sigma^{2/3}\rho_o/\alpha^{1/3}$ only (it increases with $X_o$). It is interesting to note that $X_o$ (and hence the maximal scintillation indices) depends on $\rho_o$ but not on $\ell_c$, while the distance $z$ at which the maximum of the scintillation index is reached depends on $\ell_c$ through $z_c = \ell_c/(2\sigma^{2/3}\alpha^{2/3})$. Remember
that when the initial field is a plane wave, the maximum of the scintillation index depends only on $X_c = \sigma^{2/3}\ell_c/\alpha^{1/3}$. This shows that speckled beams experience reduced intensity growth compared to plane waves: the smaller the correlation radius of the initial beam, the lower the maximal intensity reached by the beam as it propagates.

The theory (obtained in the limit $\epsilon \to 0$) is compared to simulations of Eq.(1). The intensity evolution in Fig. 3(a) exhibits distinct qualitative features with respect to the coherent excitation in Fig. 1(a). As observed experimentally [42], in the coherent case, each branch is accompanied by sidelobes arising from interference effects, which tend to disappear when the initial condition is incoherent. In the simulations we study the impact of $\rho_o$, $\sigma^2$ and $\ell_c$, on the scintillation index. Fig. 3 shows an excellent quantitative agreement, in spite of the rather limited separation of scales of the parameters. The simulations also confirm that the maximum of $S_z$ does not depend on $\ell_c$, see Fig. 3(b).

**Intensity correlation function.** Our theoretical approach can be exploited to compute an explicit form for the fourth-order moments of the field. In particular, the intensity correlation function in situation (c) is

$$C_z^{\mathcal{I},(c)}(x) = \tilde{\Pi}_{z/z_c}(x/\ell_c, 0) + \tilde{\Pi}_{z/z_c}(x/\ell_c, X_o x/\rho_o) - 1. \quad (18)$$

This expresses a two-scale behavior: At the small scale $x \sim \rho_o$, the intensity correlation function decays rapidly, while its behaviour exhibits complex variations at the large scale $x \sim \ell_c$. Similarly, the intensity correlation function in situation (pc) is given by

$$C_z^{\mathcal{I},(\text{pc})}(x) = \tilde{\Pi}_{z/z_c}(x/\ell_c, 0) - 1. \quad (19)$$

Note that it is equal to $C_z^{\mathcal{I},(c)}(x)$ when $x$ is of the order of $\ell_c$ because $\ell_c \gg \rho_o$ and $\tilde{\Pi}_{\tilde{z}}(\tilde{x}, \tilde{y}) \to 0$ as $\tilde{y} \to +\infty$. The intensity correlation functions are plotted in Fig. 4. The two-scale behavior of the intensity correlation function in situation (c) is clearly visible: the limit for large $x/\rho_o$ of the intensity correlation function is the initial value for $C_z^{\mathcal{I},(\text{pc})}$ (or $C_z^{\mathcal{I},(c)}$) when $x$ is of the order of $\ell_c$. This behavior can be seen in the numerical simulations as well in Fig. 4, and it satisfies the energy conservation relation $\int C_z^{\mathcal{I},(\text{pc})}(x)dx = 0$.

**Perspectives.** We have reported a general stochastic theory of BFs by considering both coherent and incoherent initial waves. Optics naturally provides an ideal framework for experimentally testing and observing the theoretical predictions. The results presented in the two-dimensional framework can be extended to the three-dimensional one (see [57, App. G]). Our work paves the way for a systematic approach to studying coherent phase-sensitive effects in BFs, focusing on the linear propagation regime. As shown in prior studies, linear BFs can trigger extreme nonlinear events [9–13, 61–63], with intensity peaks significantly influenced by nonlinearity [9, 64]. Our stochastic framework offers a basis for developing a theoretical model of nonlinear branched flow.

**Acknowledgements.** Funding was provided by Agence Nationale de la Recherche (Grants No. ANR-23-CE30-0021, ANR-21-ESRE-0040) and by Agence de l'Innovation de Défense (AID) via Centre Interdisciplinaire d'Études pour la Défense et la Sécurité (project PRODIPO). Calculations were performed using HPC resources from DNUM-CCUB (Université de Bourgogne Europe).

[1] M.A. Topinka, B.J. LeRoy, R.M. Westervelt, S.E.J. Shaw, R. Fleischmann, E.J. Heller, K.D. Maranowski, and A.C. Gossard, Coherent branched flow in a two-dimensional electron gas, Nature 410, 183–186 (2001).

[2] K.E. Aidala, et al. Imaging magnetic focusing of coherent electron waves, Nat. Phys. 3, 464–468 (2007).

[3] M.P. Jura, M.A. Topinka, L. Urban, A. Yazdani, H. Shtrikman, L.N. Pfeiffer, K.W. West, and D. Goldhaber-Gordon, Unexpected features of branched flow through high-mobility two-dimensional electron gases, Nat. Phys. 3, 841–845 (2007).

[4] D. Maryenko, et al., How branching can change the conductance of ballistic semiconductor devices, Phys. Rev. B 85, 195329 (2012).

[5] B. Liu and E.J. Heller, Stability of branched flow from a quantum point contact, Phys. Rev. Lett. 111, 236804 (2013).

[6] R. Höhmann, U. Kuhl, H.-J. Stöckmann, L. Kaplan, and E.J. Heller, Freak waves in the linear regime: a microwave study, Phys. Rev. Lett. 104, 093901 (2010).

[7] S. Barkhofen, J.J. Metzger, R. Fleischmann, U. Kuhl, and H.-J. Stöckmann, Experimental observation of a fundamental length scale of waves in random media, Phys. Rev. Lett. 111, 183902 (2013).

[8] E.J. Heller, R. Fleischmann, and T. Kramer, Branched flow, Physics Today 74(12), 44–51 (2021).

[9] G. Green and R. Fleischmann, Branched flow and caustics in nonlinear waves, New J. Phys. 21, 083020 (2019).

[10] Z.-Y. Sun and X. Yu, Nonlinear Schrödinger waves in a disordered potential: Branched flow, spectrum diffusion, and rogue waves, Chaos 32, 023108 (2022).

[11] K. Jiang, T.W. Huang, R. Li, M.Y. Yu, H.B. Zhuo, S.Z. Wu, C.T. Zhou, and S.C. Ruan, Branching of high-current relativistic electron beam in porous materials, Phys. Rev. Lett. 130, 185001 (2023).

[12] M. Mattheakis and G.P. Tsironis, Quodons in mica. In Extreme Waves and Branched Flows in Optical Media 425–454 (Springer International Publishing, 2015).

[13] M. Mattheakis, I.J. Pitsios, G.P. Tsironis, and S. Tzortzakis, Extreme events in complex linear and nonlinear photonic media. Chaos Solitons Fractals 84, 73–80 (2016).

[14] B.S. White and B. Fornberg, On the chance of freak waves at sea, J. Fluid Mech. 355, 113–138 (1998).

[15] M.V. Berry, Tsunami asymptotics, New J. Phys. 7, 129 (2005).

[16] M.V. Berry, Focused tsunami waves, Proc. R. Soc. A 463, 3055–3071 (2007).

[17] E.J. Heller, L. Kaplan, and A. Dahlen, Refraction of a Gaussian Seaway, J. Geophys. Res. 113, C09023 (2008).

[18] L.H. Ying, Z. Zhuang, E.J. Heller, and L. Kaplan, Linear and nonlinear rogue wave statistics in the presence of random currents. Nonlinearity 24, R67 (2011).

[19] H. Degueldre, J.J. Metzger, T. Geisel, and R. Fleischmann, Random focusing of tsunami waves. Nat. Phys. 12, 259–262 (2016).

[20] M.A. Wolfson and S. Tomsovic, On the stability of long-range
sound propagation through a structured ocean. J. Acoust. Soc. Am. **109**, 2693–2703 (2001).

[21] M. Mattheakis, G.P. Tsironis, and E. Kaxiras, Emergence and dynamical properties of stochastic branching in the electronic flows of disordered Dirac solids. EPL **122**, 27003 (2018).

[22] K. Jose, N. Ferguson, and A. Bhaskar, Branched flows of flexural waves in non-uniform elastic plates Comm. Phys. **5**, 152 (2023).

[23] A. Brandstötter, A. Girschik, P. Ambichl, and S. Rotter, Shaping the branched flow of light through disordered media, Proc. Natl. Acad. Sci. USA **116**(27), 13260–13265 (2019).

[24] J. St'avina and P. Bokes, Quantum and classical branching flow in space and time, Phys. Rev. A **106**, 052215 (2022).

[25] A. Daza, E.J. Heller, A.M. Graf, and E. Räsänen, Propagation of waves in high Brillouin zones: Chaotic branched flow and stable superwires, Proc. Natl. Acad. Sci. USA **118**(40), e2110285118 (2021).

[26] A. Wagemakers, A. Hartikainen, A. Daza, E. Räsänen, and M.A.F. Sanjuán, Chaotic dynamics creates and destroys branched flow, Phys. Rev. E **111**, 014214 (2025).

[27] K.H. Mok and R. Fleischmann, Branched flows in active random walks and the formation of ant trail patterns, Phys. Rev. Res. **5**, 043299 (2023).

[28] A. Patsyk, U. Sivan, M. Segev, and M.A. Bandres, Observation of branched flow of light, Nature **583**, 60–65 (2020).

[29] S. Chang, K.-H. Wu, S. Liu, Z.-K. Lin, J. Wu, S. Ge, L.-J. Chen, P. Chen, W. Hu, Y. Xu, H. Chen, D. He, D.-Q. Yang, J.-H. Jiang, Y. Lu, and J. Chen, Electrical tuning of branched flow of light, Nature Comm. **15**, 197 (2024).

[30] S. Rotter and S. Gigan, Light fields in complex media: Meso-scopic scattering meets wave control, Rev. Mod. Phys. **89**, 015005 (2017).

[31] H. Cao, A.P. Mosk, and S. Rotter, Shaping the propagation of light in complex media, Nature Physics **18**, 994–1007 (2022).

[32] V. Kulkarny and B.S. White, Focusing of waves in turbulent inhomogeneous media, Phys. Fluids **25**, 1770–1784 (1982).

[33] D.I. Zwillinger and B.S. White, Propagation of initially plane waves in the region of random caustics, Wave Motion **7**, 207–227 (1985).

[34] M.V. Berry and C. Upstill, IV catastrophe optics: morphologies of caustics and their diffraction patterns, Progress in Optics **18**, 257–346 (1980).

[35] J.F. Nye, Natural Focusing and Fine Structure of Light: Caustics and Wave Dislocations, CRC Press, 1999.

[36] L. Kaplan, Statistics of branched flow in a weak correlated random potential, Phys. Rev. Lett. **89**, 184103 (2002).

[37] J.J. Metzger, R. Fleischmann, and T. Geisel, Universal statistics of branched flows, Phys. Rev. Lett. **105**, 020601 (2010).

[38] J.J. Metzger, R. Fleischmann, and T. Geisel, Statistics of extreme waves in random media, Phys. Rev. Lett. **112**, 203903 (2014).

[39] M. Pradas, A. Pumir, and M. Wilkinson, Uniformity transition for ray intensities in random media, J. Phys. A **51**, 155002 (2018).

[40] J.J. Metzger, R. Fleischmann, and T. Geisel, Intensity fluctuations of waves in random media: What is the semiclassical limit?, Phys. Rev. Lett. **111**, 013901 (2013).

[41] M.V. Berry, Elementary branching: Waves, rays, decoherence, J. Opt. **22**, 115608 (2020).

[42] A. Patsyk, Y. Sharabi, U. Sivan, and M. Segev, Incoherent branched flow of light, Phys. Rev. X **12**, 021007 (2022).

[43] L.C. Andrews and R.L. Phillips, Laser Beam Propagation through Random Media, SPIE-International Society for Optical Engineering, 2005.

[44] V.I. Tatarskii, Waves Propagation in a Turbulent Medium, McGraw-Hill, 1961.

[45] B. Crosignani, B. Daino, and P. Di Porto, Light scattering by a rotating disk, J. Appl. Phys. **42**, 399–403 (1971).

[46] J.W. Goodman, Speckle Phenomena in Optics: Theory and Applications, 2nd ed. (SPIE Press, 2020).

[47] B.I. Erkmen and J.H. Shapiro, Unified theory of ghost imaging with Gaussian-state light, Phys. Rev. A **77**, 043809 (2008).

[48] J.H. Shapiro and R.W. Boyd, The physics of ghost imaging, Quantum Inf. Process. **11**, 949–993 (2012).

[49] A. Ishimaru, Wave Propagation and Scattering in Random Media, Academic Press, 1978.

[50] J. Garnier and K. Sølna, Scintillation in the white-noise paraxial regime, Communications in Partial Differential Equations **39**, 626–650 (2014).

[51] J. Garnier and K. Sølna, Fourth-order moments analysis for partially coherent electromagnetic beams in random media, Waves in Random and Complex Media **33**, 1346–1365 (2023).

[52] This results from the fact that, for z large, the correlation function has the form $\exp(-\gamma_2 z y^2/2)$, where $\gamma_2 = -\partial_x^2 \gamma(0) = \sigma^2 \ell_c^{-1} \tilde{\gamma}_2$, $\tilde{\gamma}_2 = -\partial_z^2 \gamma(0)$, is always positive because it is proportional to the power spectral density of the stationary process $\partial_x V(z, 0)$ at frequency 0.

[53] G. Strang, On the construction and comparison of difference schemes, SIAM J. Num. Anal. **5**, 506–517 (1968).

[54] J. Garnier and K. Sølna, Fourth-moment analysis for beam propagation in the white-noise paraxial regime, Archive on Rational Mechanics and Analysis **220**, 37–81 (2016).

[55] G. Bal and A. Nair, Complex Gaussianity of long-distance random wave processes, arXiv:2402.17107.

[56] G. Bal and A. Nair, Long distance propagation of light in random media with partially coherent sources, arXiv:2406.05252.

[57] See Supplementary Material for the theoretical developments and details about the numerical simulations.

[58] J.T. Foley and M.S. Zubairy, The directionality of Gaussian Schell-model beams, Optics Communications **26**, 297–300 (1978).

[59] A.T. Friberg and R.J. Sudol, Propagation parameters of Gaussian Schell-model beams, Optics Communications **41**, 383–387 (1982).

[60] J.-P. Fouque, J. Garnier, G. Papanicolaou, and K. Sølna, Wave Propagation and Time Reversal in Randomly Layered Media, Springer, 2007.

[61] J.M. Dudley, F. Dias, M. Erkintalo, and G. Genty, Instabilities, breathers and rogue waves in optics, Nat. Photon. **8**, 755–764 (2014).

[62] A. Safari, R. Fickler, M.J. Padgett, and R.W. Boyd, Generation of caustics and rogue waves from nonlinear instability, Phys. Rev. Lett. **119**, 203901 (2017).

[63] J.M. Dudley, G. Genty, A. Mussot, A. Chabchoub, and F. Dias, Rogue waves and analogies in optics and oceanography, Nature Reviews Physics **1**, 675 (2019).

[64] K. Jiang, T.W. Huang, C.N. Wu, M. Y. Yu, H. Zhang, S.Z. Wu, H.B. Zhuo, A. Pukhov, C.T. Zhou, and S.C. Ruan, Non-linear branched flow of intense laser light in randomly uneven media, Matter Radiat. Extremes **8**, 024402 (2023).
# Supplementary Material on the Article

## Appendix A: Ray theory of branched flow

One of the main tools used to describe and understand the properties of branched flow has been the ray tracing method. Rays are constructed as the characteristic curves of the eikonal equation obtained by considering rapidly oscillating solution of the paraxial Eq.(1) (main text) [S1]

$$ \psi_z(x) = A(z, x)e^{iS(z, x)/\delta}, \quad (A1) $$

where $\delta$ is an order counting parameter that expresses the assumption that the local phase $S(z, x)$ of the wave varies rapidly in comparison with its amplitude $A(z, x)$. Expanding the wave equation in a hierarchical fashion in powers of $\delta^{-1}$ yields a system of equations for the local phase and amplitude as asymptotic series in $\delta$. The leading-order equation is the eikonal equation:

$$ -\partial_z S_0 = \alpha (\partial_x S_0)^2 + V(z, x), \quad (A2) $$

with $S_0$ the first term in the asymptotic series of $S$. The eikonal equation can be interpreted as a Hamilton-Jacobi equation for a Hamiltonian $H(z, x, k_z, k_x)$ that determines the ray trajectories. The Hamiltonian is obtained by making the substitution $(\partial_z S_0, \partial_x S_0) \to (k_z, k_x)$:

$$ H = k_z + \alpha k_x^2 + V(z, x). \quad (A3) $$

The condition $H = 0$ is the local dispersion relation.

The rays are parametrized curves, $(x(s), z(s))$, that are the solutions of Hamilton's equations:

$$ \dot{x} = \partial_{k_x} H, \quad \dot{z} = \partial_{k_z} H, \quad \dot{k}_x = -\partial_x H, \quad \dot{k}_z = -\partial_z H, \quad (A4) $$

where the dot represents a derivative with respect to the ray parameter $s$. Since the Hamiltonian is linear in $k_z$, we have $\dot{z} = 1$ and therefore we can choose the $z$ coordinate to parametrize the rays. Plugging the expression for the Hamiltonian, the equations to be solved are

$$ \dot{x} = 2\alpha k_x, \quad \dot{k}_x = -\partial_x V, \quad \dot{k}_z = -\partial_z V, \quad (A5) $$

together with the initial conditions $(x(z = 0), k_x(z = 0), k_z(z = 0))$. From the definition of the momentum, we have $k_x(z = 0) = \partial_x S(z = 0, x)$ and $k_z$ is then found by solving the condition $H = 0$.

In the case of a coherent plane wave, $S(z = 0, x)$ is constant and therefore $k_x(z = 0) = 0$ for all rays (all the rays are parallel to each other). In the case of an incoherent initial speckle field, $k_x \neq 0$ and the launched rays are not parallel. Figures 5-6 report the ray dynamics and the solution of the paraxial Eq.(1) (main text) for the same realization of the random potential and for different initial conditions: coherent plane wave in Fig. 5, and incoherent speckled field in Fig. 6 (panels (a) and (b) report two different realizations of the speckle field). We can clearly see the formation of causatics which are associated with increases in wave intensities. This illustrates the fact that the rays equations can predict

<figure><img src="image_7.png" /><figcaption>FIG. 5: Coherent initial plane wave. Example of a congruence of rays for a single realization of the random potential. Panel (a) shows a realization of the random potential $V(z, x)$ together with the rays associated with an initial coherent plane wave at $z = 0$. Panel (b) shows the BF obtained by solving the paraxial Eq.(1) (main text) with the potential shown in (a), superimposed with the congruence of rays. We can see the formation of caustics which are associated with increases in wave intensity. Parameters are the same as in Fig. 1(a) (main text).</figcaption></figure>

the positions of the maximal intensities, but the values of the maxima result from interference effects that depend on the coherence properties of the initial field. This was discussed in [S2] by using a simple random phase model. The determination of the statistics of the values of the intensity maxima require a detailed multiscale analysis as carried out in the main text.

## Appendix B: Scaled Vlasov equation

We consider the Wigner transform $W_z^\epsilon(x, k)$ defined by Eq.(12) (main text). From Eq.(10) (main text), it satisfies the scaled Vlasov-type equation

$$ \partial_z W_z^\epsilon + \epsilon^{1-d-b} \partial_k \omega_k \partial_x W_z^\epsilon + \epsilon^{c-b-1} i \int_{\mathbb{R}} \left[ \mathcal{V}(\frac{z}{\epsilon^b}, x + \epsilon^d \frac{y}{2}) - \mathcal{V}(\frac{z}{\epsilon^b}, x - \epsilon^d \frac{y}{2}) \right] \times \langle \psi_{\frac{z}{\epsilon^b}}^\epsilon (x + \epsilon^d \frac{y}{2}) \overline{\psi_{\frac{z}{\epsilon^b}}^\epsilon (x - \epsilon^d \frac{y}{2})} \rangle \exp(-iky) dy = 0, $$
<figure><img src="image_8.png" /><figcaption>FIG. 6: Incoherent initial wave. (a) Evolution of the BF obtained by solving the paraxial Eq.(1) (main text), superimposed with the congruence of rays. (b) Same as in (a), except that a different realization of the initial incoherent speckle field is considered. The realization of the random potential V(z, x) in (a) and (b) is the same as in Fig. 5. As in the case of the initial plane wave, we can see the formation of caustics despite the fact that the rays are not all launched parallel to each other due to the initial phase variations. Parameters are the same as in Fig. 3(a) (main text)</figcaption></figure>

with $\omega_k = \alpha k^2$, which gives after expansion of the last term of the left-hand side

$$ \partial_z W_z^\varepsilon + \varepsilon^{1-d-b} \partial_k \omega_k \partial_x W_z^\varepsilon - \varepsilon^{c+d-b-1} \partial_x V(\frac{z}{\varepsilon^b}, x) \partial_k W_z^\varepsilon = O(\varepsilon^{c+3d-b-1}), \quad (B1) $$

with the initial condition $W_{z=0}^{\varepsilon}(x, k) = W_o(k)$. In the scaling regime $d \in (1/5, 1)$, $b = 1-d$, $c = 3(1-d)/2$, we have $1-d-b=0$, $c+d-b-1=-b/2$, and $c+3d-b-1=(5d-1)/2 > 0$, so that we can neglect the remainder in (B1) and we get that $W_z^\varepsilon$ satisfies (13) (main text) with the initial condition $W_{z=0}^\varepsilon(x, k) = W_o(k)$.

Note that the scaling regime addressed here is different from the one used to derive the paraxial white-noise (or Itô-Schrödinger) model [S3, S4, S5, S6]. The paraxial white-noise model is valid when $d=0$, $b=1$, $c=3/2$, that is to say, when the wavelength is much smaller than the correlation radius of the medium, which is itself of the same order as the correlation radius of the initial field.

## Appendix C: Proof of Equation (15) (main text)

We have

$$ W_z^{\varepsilon}(X, K) = \int_{\mathbb{R}^2} W_z^{\varepsilon}(x', k') \delta(x' - X) \delta(k' - K) dx' dk'. $$

We make the change of variables $(x', k') \mapsto (x, k)$ with $x' = X_z^\varepsilon(x, k)$, $k' = K_z^\varepsilon(x, k)$:

$$ W_z^{\varepsilon}(X, K) = \int_{\mathbb{R}^2} W_z^{\varepsilon}(X_z^{\varepsilon}(x, k), K_z^{\varepsilon}(x, k)) \delta(X_z^{\varepsilon}(x, k) - X) \times \delta(K_z^{\varepsilon}(x, k) - K) |\det \mathbf{J}_z^{\varepsilon}(x, k)| dx dk, $$

where $\mathbf{J}_z^\varepsilon(x, k)$ is the Jacobian

$$ \mathbf{J}_z^{\varepsilon}(x, k) = \begin{pmatrix} \frac{\partial X_z^{\varepsilon}}{\partial x} & \frac{\partial X_z^{\varepsilon}}{\partial k} \\ \frac{\partial K_z^{\varepsilon}}{\partial x} & \frac{\partial K_z^{\varepsilon}}{\partial k} \end{pmatrix}. $$

On the one hand we have $W_z^\varepsilon(X_z^\varepsilon(x, k), K_z^\varepsilon(x, k)) = W_o(k)$ and on the other hand we can compute

$$ \begin{aligned} \frac{d}{dz} \frac{\partial X_z^\varepsilon}{\partial x} &= 2\alpha \frac{\partial K_z^\varepsilon}{\partial x}, & \frac{\partial X_z^\varepsilon}{\partial x} \Big|_{z=0} (x, k) &= 1, \\ \frac{d}{dz} \frac{\partial K_z^\varepsilon}{\partial x} &= -\frac{1}{\varepsilon^{b/2}} \partial_x^2 V(\frac{z}{\varepsilon^b}, X_z^\varepsilon) \frac{\partial X_z^\varepsilon}{\partial x}, & \frac{\partial K_z^\varepsilon}{\partial x} \Big|_{z=0} (x, k) &= 0, \\ \frac{d}{dz} \frac{\partial X_z^\varepsilon}{\partial k} &= 2\alpha \frac{\partial K_z^\varepsilon}{\partial k}, & \frac{\partial X_z^\varepsilon}{\partial k} \Big|_{z=0} (x, k) &= 0, \\ \frac{d}{dz} \frac{\partial K_z^\varepsilon}{\partial k} &= -\frac{1}{\varepsilon^{b/2}} \partial_x^2 V(\frac{z}{\varepsilon^b}, X_z^\varepsilon) \frac{\partial X_z^\varepsilon}{\partial k}, & \frac{\partial K_z^\varepsilon}{\partial k} \Big|_{z=0} (x, k) &= 1, \end{aligned} $$

which gives

$$ \begin{aligned} \frac{d}{dz} \mathrm{Det} \mathbf{J}_z^\varepsilon &= \left( \frac{d}{dz} \frac{\partial X_z^\varepsilon}{\partial x} \right) \frac{\partial K_z^\varepsilon}{\partial k} + \frac{\partial X_z^\varepsilon}{\partial x} \left( \frac{d}{dz} \frac{\partial K_z^\varepsilon}{\partial k} \right) \\ &\quad - \left( \frac{d}{dz} \frac{\partial K_z^\varepsilon}{\partial x} \right) \frac{\partial X_z^\varepsilon}{\partial k} - \frac{\partial K_z^\varepsilon}{\partial x} \left( \frac{d}{dz} \frac{\partial X_z^\varepsilon}{\partial k} \right) \\ &= 2\alpha \frac{\partial K_z^\varepsilon}{\partial x} \frac{\partial K_z^\varepsilon}{\partial k} - \frac{1}{\varepsilon^{b/2}} \partial_x^2 \mathcal{V}(\frac{z}{\varepsilon^b}, X_z^\varepsilon) \frac{\partial X_z^\varepsilon}{\partial x} \frac{\partial X_z^\varepsilon}{\partial k} \\ &\quad + \frac{1}{\varepsilon^{b/2}} \partial_x^2 \mathcal{V}(\frac{z}{\varepsilon^b}, X_z^\varepsilon) \frac{\partial X_z^\varepsilon}{\partial x} \frac{\partial X_z^\varepsilon}{\partial k} - 2\alpha \frac{\partial K_z^\varepsilon}{\partial x} \frac{\partial K_z^\varepsilon}{\partial k} \\ &= 0, \end{aligned} $$

hence $\mathrm{Det} \mathbf{J}_z^\varepsilon = \mathrm{Det} \mathbf{J}_{z=0}^\varepsilon = \mathrm{Det} \mathbf{I} = 1$. This gives the desired result Eq.(15) (main text).

## Appendix D: Diffusion approximation theory

This section contains the technical results that are needed to characterize the statistics of the wave field, in particular the width of the envelope, the correlation radius of the field and the scintillation index. The potential V is a smooth, stationary, random process with mean zero and integrable covariance function. Applying diffusion-approximation theory [S7, Chapter 6], we can show from (14) (main text) that, for any integer n, for any $x_1, \dots, x_n \in \mathbb{R}$, for any $k_1, \dots, k_n \in \mathbb{R}$, the $\mathbb{R}^{2n}$-valued process $(X_z^\varepsilon(x_j, k_j), K_z^\varepsilon(x_j, k_j))_{j=1}^n$ converges
in distribution as $\varepsilon \to 0$ to the Markov diffusion process $(X_z(x_j, k_j), K_z(x_j, k_j))_{j=1}^n$ with the infinitesimal generator

$$ \mathcal{L}^{(n)} = \sum_{j=1}^{n} 2\alpha K_j \frac{\partial}{\partial X_j} + \frac{1}{2} \sum_{j,j'=1}^{n} \Gamma(X_j - X_{j'}) \frac{\partial^2}{\partial K_j \partial K_{j'}} \quad (\text{D1}) $$

where

$$ \Gamma(x) = \int_{-\infty}^{\infty} \mathbb{E}[\partial_x V(0,0) \partial_x V(z,x)] dz. \quad (\text{D2}) $$

As a particular example of smooth random medium, we can consider a potential $V$ with Gaussian correlation function, variance $\sigma^2$ and correlation radius $\ell_c$. We then have

$$ \Gamma(x) = 2\sqrt{\pi}\sigma^2 \ell_c^{-1} \left(1 - \frac{2x^2}{\ell_c^2}\right) \exp\left(-\frac{x^2}{\ell_c^2}\right). \quad (\text{D3}) $$

**Application** $n=1$. Let $x, k \in \mathbb{R}$. The pdf $p_z^{(1)}(X, K; x, k)$ of $(X_z(x, k), K_z(x, k))$ satisfies the Fokker-Planck equation

$$ \partial_z p_z^{(1)} = (\mathcal{L}^{(1)})^* p_z^{(1)}, \quad (\text{D4}) $$

starting from $p_{z=0}^{(1)}(X, K; x, k) = \delta(X - x)\delta(K - k)$, where

$$ \mathcal{L}^{(1)} = 2\alpha K \partial_X + \frac{\Gamma(0)}{2} \partial_K^2, $$

and $(\mathcal{L}^{(1)})^*$ is the adjoint of $\mathcal{L}^{(1)}$. Eq. (D4) has the form

$$ \partial_z p_z^{(1)} = -2\alpha K \partial_X p_z^{(1)} + \frac{\Gamma(0)}{2} \partial_K^2 p_z^{(1)}. \quad (\text{D5}) $$

It is possible to solve this equation (by taking a Fourier transform in $(X, K)$) and we get the expression of the pdf of the limit process $(X_z(x, k), K_z(x, k))$:

$$ \begin{aligned} p_z^{(1)}(X, K; x, k) = & \frac{1}{\sqrt{2\pi\Gamma(0)z}} \exp\left(-\frac{(K-k)^2}{2\Gamma(0)z}\right) \frac{1}{\sqrt{2\pi\Gamma(0)\frac{z^3}{3}}} \\ & \times \exp\left(-\frac{3(X-x-\alpha(K+k)z)^2}{2\Gamma(0)z^3}\right). \end{aligned} \quad (\text{D6}) $$

**Application** $n=2$. Let $x_1, x_2, k_1, k_2 \in \mathbb{R}$. The pdf $p_z^{(1)}(X_1, X_2, K_1, K_2; x_1, x_2, k_1, k_2)$ of $(X_z(x_j, k_j), K_z(x_j, k_j))_{j=1}^2$ satisfies the Fokker-Planck equation

$$ \partial_z p_z^{(2)} = (\mathcal{L}^{(2)})^* p_z^{(2)}, \quad (\text{D7}) $$

starting from $p_{z=0}^{(2)}(X_1, X_2, K_1, K_2; x_1, x_2, k_1, k_2) = \delta(X_1 - x_1)\delta(X_2 - x_2)\delta(K_1 - k_1)\delta(K_2 - k_2)$, where $\mathcal{L}^{(2)}$ is the infinitesimal generator of $(X_z(x_j, k_j), K_z(x_j, k_j))_{j=1}^2$:

$$ \begin{aligned} \mathcal{L}^{(2)} = & 2\alpha K_1 \frac{\partial}{\partial X_1} + 2\alpha K_2 \frac{\partial}{\partial X_2} + \frac{1}{2}\Gamma(0) \left( \frac{\partial^2}{\partial K_1^2} + \frac{\partial^2}{\partial K_2^2} \right) \\ & + \Gamma(X_1 - X_2) \frac{\partial^2}{\partial K_1 \partial K_2}. \end{aligned} \quad (\text{D8}) $$

We introduce

$$ R = \frac{X_1 + X_2}{2}, \quad Q = X_1 - X_2, \quad (\text{D9}) $$

$$ U = \frac{K_1 + K_2}{2}, \quad V = K_1 - K_2, \quad (\text{D10}) $$

where $X_j = X_z(x_j, k_j)$, $K_j = K_z(x_j, k_j)$, $j = 1, 2$. The infinitesimal generator of the process $(R_z, Q_z, U_z, V_z)$ is

$$ \mathcal{L} = 2\alpha U \partial_R + 2\alpha V \partial_Q + \frac{1}{4}(\Gamma(0) + \Gamma(Q))\partial_U^2 + (\Gamma(0) - \Gamma(Q))\partial_V^2. \quad (\text{D11}) $$

In particular, the process $(Q_z, V_z)$ is Markov with generator

$$ \mathcal{L} = 2\alpha V \partial_Q + (\Gamma(0) - \Gamma(Q)) \partial_V^2. \quad (\text{D12}) $$

## Appendix E: Expression of the scintillation index for incoherent initial conditions

From Eq. (15) (main text) we get the expression of the second-order moment of the Wigner transform in the limit $\varepsilon \to 0$:

$$ \begin{aligned} & \lim_{\varepsilon \to 0} \mathbb{E}[W_z^{\varepsilon}(X_1, K_1) W_z^{\varepsilon}(X_2, K_2)] \\ & = \int_{\mathbb{R}^4} \mathcal{W}_o(k_1) \mathcal{W}_o(k_2) \\ & \quad \times p_z^{(2)}(X_1, K_1, X_2, K_2; x_1, k_1, x_2, k_2) dx_1 dk_1 dx_2 dk_2, \end{aligned} \quad (\text{E1}) $$

where $p_z^{(2)}$ is the solution of the Fokker-Planck equation (D7). The second-order moment of the intensity in situation (pc) is

$$ \begin{aligned} & \mathbb{E}\left[\left|\psi_{\frac{\varepsilon}{\varepsilon_b}}^{\varepsilon}(X)\right|^2\right] \\ & = \frac{1}{(2\pi)^2} \int_{\mathbb{R}^2} \mathbb{E}[W_z^{\varepsilon}(X, K_1) W_z^{\varepsilon}(X, K_2)] dK_1 dK_2. \end{aligned} \quad (\text{E2}) $$

The second-order moment of the intensity in situation (c) is

$$ \begin{aligned} & \mathbb{E}\left[\left|\psi_{\frac{\varepsilon}{\varepsilon_b}}^{\varepsilon}(X)\right|^4\right] \\ & = \frac{2}{(2\pi)^2} \int_{\mathbb{R}^2} \mathbb{E}[W_z^{\varepsilon}(X, K_1) W_z^{\varepsilon}(X, K_2)] dK_1 dK_2, \end{aligned} \quad (\text{E3}) $$

where we have used Isserlis' theorem [S8]

$$ \langle \psi_o^\epsilon(x) \bar{\psi}_o^\epsilon(y) \psi_o^\epsilon(x') \bar{\psi}_o^\epsilon(y') \rangle = \langle \psi_o^\epsilon(x) \bar{\psi}_o^\epsilon(y) \rangle \langle \psi_o^\epsilon(x') \bar{\psi}_o^\epsilon(y') \rangle \\ + \langle \psi_o^\epsilon(x) \bar{\psi}_o^\epsilon(y') \rangle \langle \psi_o^\epsilon(x') \bar{\psi}_o^\epsilon(y) \rangle. $$

By Eq. (E1) and the change of variables (D9-D10) we then get:

$$ \begin{aligned} & \lim_{\epsilon \to 0} \mathbb{E}\left[\left|\psi_{\frac{\epsilon}{\epsilon_b}}^{\epsilon}(R)\right|^2\right] \\ & = \frac{1}{(2\pi)^2} \int_{\mathbb{R}^6} \mathcal{W}_o\left(u + \frac{v}{2}\right) \mathcal{W}_o\left(u - \frac{v}{2}\right) \\ & \quad \times p_z(R, 0, U, V | r, q, u, v) dU dV dr dq dudv \\ & = \frac{1}{(2\pi)^2} \int_{\mathbb{R}} \left[ \int_{\mathbb{R}} \mathcal{W}_o\left(u + \frac{v}{2}\right) \mathcal{W}_o\left(u - \frac{v}{2}\right) du \right] \\ & \quad \times \left[ \int_{\mathbb{R}^2} p_z(0, V | q, v) dq dV \right] dv, \end{aligned} \quad (\text{E4}) $$
which does not depend on $R$. By (D12), in the last line $p_z(Q, V|q, v)$ is the pdf solution of

$$ \partial_z p_z = -2\alpha V \partial_Q p_z + (\Gamma(0) - \Gamma(Q)) \partial_V^2 p_z, \quad (E5) $$

starting from $p_{z=0}(Q, V|q, v) = \delta(Q - q)\delta(V - v)$. Eq. (E4) can be rewritten as

$$ \lim_{\epsilon \to 0} \mathbb{E} \left[ \left\langle |\psi_{\tilde{z}}^{\epsilon}(R)|^2 \right\rangle^2 \right] = \Pi_z(0, 0) \quad (E6) $$

in terms of the function $\Pi_z$ defined by

$$ \Pi_z(Q, S) = \int_{\mathbb{R}^3} p_z(Q, V|q, v) \pi_o(v) e^{iSv} dq dvdV, \quad (E7) $$

with $\pi_o(v) = \frac{1}{(2\pi)^2} \int_{\mathbb{R}} W_o(u + \frac{v}{2})W_o(u - \frac{v}{2})du$. The function $\Pi_z$ is the solution of

$$ \partial_z \Pi_z = 2i\alpha \partial_Q \partial_S \Pi_z - (\Gamma(0) - \Gamma(Q)) S^2 \Pi_z, \quad (E8) $$

starting from $\Pi_{z=0}(Q, S) = \int \pi_o(v) e^{iSv} dv = |C_o(S)|^2$. This gives Eqs. (16-17) (main text).

## Appendix F: Proof of the small z-expansion

Let $\tilde{\Pi}_{\tilde{z}}$ be the solution of (17) (main text). We consider the functions

$$ \tilde{M}_{j, \tilde{z}}(\tilde{x}) = (-i)^j \partial_{\tilde{y}}^j \tilde{\Pi}_{\tilde{z}}(\tilde{x}, \tilde{y}) \bigg|_{\tilde{y}=0}. $$

They satisfy the equations

$$ \partial_{\tilde{z}} \tilde{M}_{0, \tilde{z}} = -\partial_{\tilde{x}} \tilde{M}_{1, \tilde{z}}, $$

$$ \partial_{\tilde{z}} \tilde{M}_{1, \tilde{z}} = -\partial_{\tilde{x}} \tilde{M}_{2, \tilde{z}}, $$

$$ \partial_{\tilde{z}} \tilde{M}_{2, \tilde{z}} = -\partial_{\tilde{x}} \tilde{M}_{3, \tilde{z}} + (\tilde{\Gamma}(0) - \tilde{\Gamma}(\tilde{x})) \tilde{M}_{0, \tilde{z}}, $$

starting from $\tilde{M}_{j, \tilde{z}=0}(\tilde{x}) = \tilde{M}_{j,o} := (-iX_o)^j \tilde{\pi}_o^{(j)}(0)$. For small $\tilde{z}$ and using the fact that $\tilde{M}_{j,o}$ does not depend on $\tilde{x}$, we get successively:

$$ \tilde{M}_{2, \tilde{z}}(\tilde{x}) = \tilde{M}_{2,o} + (\tilde{\Gamma}(0) - \tilde{\Gamma}(\tilde{x})) \tilde{M}_{0,o} \tilde{z} + o(\tilde{z}), $$

$$ \tilde{M}_{1, \tilde{z}}(\tilde{x}) = \tilde{M}_{1,o} + \frac{1}{2} \partial_{\tilde{x}} \tilde{\Gamma}(\tilde{x}) \tilde{M}_{0,o} \tilde{z}^2 + o(\tilde{z}^2), $$

$$ \tilde{M}_{0, \tilde{z}}(\tilde{x}) = \tilde{M}_{0,o} - \frac{1}{6} \partial_{\tilde{x}}^2 \tilde{\Gamma}(\tilde{x}) \tilde{M}_{0,o} \tilde{z}^3 + o(\tilde{z}^3). $$

By Eq. (16) (main text) this gives the desired result for the small z-expansions of $S_z^{(c)}$ and $S_z^{(pc)}$ since $\tilde{M}_{0,o} = 1$. More specifically, we get $S_z^{(pc)} = \frac{\tilde{\gamma}_4}{6} \frac{\tilde{z}^3}{\tilde{z}_c^3} + o(\frac{\tilde{z}^3}{\tilde{z}_c^3})$, $S_z^{(c)} = 1 + \frac{\tilde{\gamma}_4}{3} \frac{\tilde{z}^3}{\tilde{z}_c^3} + o(\frac{\tilde{z}^3}{\tilde{z}_c^3})$, with $\tilde{\gamma}_4 = -\partial_{\tilde{x}}^2 \tilde{\Gamma}(0) = \partial_{\tilde{x}}^4 \tilde{\gamma}(0)$, where $\gamma(x) = \int_{\mathbb{R}} \mathbb{E}[V(0,0)V(z,x)]dz = \sigma^2 \ell_c \tilde{\gamma}(x/\ell_c)$. For a medium with Gaussian correlation function, $\tilde{\gamma}(\tilde{x}) = \sqrt{\pi} \exp(-\tilde{x}^2)$ and $\tilde{\gamma}_4 = 12\sqrt{\pi}$.

## Appendix G: Extension to the three-dimensional case

The results described in this paper can be readily extended to the three-dimensional paraxial wave equation:

$$ i \partial_z \psi_z = -\alpha (\partial_{x_1}^2 + \partial_{x_2}^2) \psi_z + V(z, \mathbf{x}) \psi, \quad (G1) $$

for $z > 0$, $\mathbf{x} = (x_1, x_2) \in \mathbb{R}^2$. As an illustration, let us assume that that the initial field is Gaussian with Wigner transform independent of $\mathbf{x}$:

$$ \int_{\mathbb{R}^2} \langle \psi_o(\mathbf{x} + \frac{\mathbf{y}}{2}) \overline{\psi_o(\mathbf{x} - \frac{\mathbf{y}}{2})} \rangle e^{-i\mathbf{k} \cdot \mathbf{y}} d\mathbf{y} = \tilde{W}_o(\rho_o \mathbf{k}), \quad (G2) $$

and that the medium fluctuations have Gaussian covariance function:

$$ \mathbb{E}[V(0, \mathbf{0})V(z, \mathbf{x})] = \sigma^2 \exp \left( - \frac{|\mathbf{x}|^2 + z^2}{\ell_c^2} \right). \quad (G3) $$

Under such circumstances, in the situation (pc) the scintillation index defined by (4) (main text) has the form

$$ S_z^{(\text{pc})} = \tilde{\Pi}_{z/z_c}(\mathbf{0}, \mathbf{0}) - 1, \quad (G4) $$

while in the situation (c) the scintillation index defined by (3) (main text) has the form

$$ S_z^{(c)} = 2\tilde{\Pi}_{z/z_c}(\mathbf{0}, \mathbf{0}) - 1. \quad (G5) $$

The scintillation index in situations (c) and (pc) depends on the function $\tilde{\Pi}_{\tilde{z}}$ that is the solution of:

$$ \partial_{\tilde{z}} \tilde{\Pi}_{\tilde{z}} = i \nabla_{\tilde{x}} \cdot \nabla_{\tilde{y}} \tilde{\Pi}_{\tilde{z}} - \frac{1}{2} \sum_{j,l=1}^{2} (\tilde{\Gamma}_{jl}(\mathbf{0}) - \tilde{\Gamma}_{jl}(\tilde{x})) \tilde{y}_j \tilde{y}_l \tilde{\Pi}_{\tilde{z}}, \quad (G6) $$

starting from $\tilde{\Pi}_{\tilde{z}=0}(\tilde{x}, \tilde{y}) = |\tilde{C}_o(\tilde{y}/X_o)|^2/\tilde{C}_o(\mathbf{0})^2$, where $\tilde{C}_o$ is the inverse Fourier transform of $\tilde{W}_o$ and

$$ \tilde{\Gamma}(\tilde{x}) = 2\sqrt{\pi} \left( \mathbf{I} - 2 \begin{pmatrix} \tilde{x}_1^2 & \tilde{x}_1 \tilde{x}_2 \\ \tilde{x}_1 \tilde{x}_2 & \tilde{x}_2^2 \end{pmatrix} \right) \exp \left( - |\tilde{x}|^2 \right). \quad (G7) $$

## Appendix H: Normalization and simulations

We performed numerical simulations of the paraxial wave equation Eq. (1) (main text) by normalizing the spatial variables with respect to the wavelength $\lambda$:

$$ i \partial_{z'} \psi_{z'}(x') = -\alpha' \partial_{x'}^2 \psi_{z'} + V'(z', x') \psi_{z'}, \quad (H1) $$

where $x' = x/\lambda$, $z' = z/\lambda$, $V' = \lambda V = \pi(n_o^2 - n^2(z', x'))/n_o$, and $\alpha' = \alpha/\lambda = 1/(4\pi n_o) \approx 0.053$ with a reference refractive index of $n_o = 1.5$. Accordingly, the normalized initial correlation length is $\rho'_o = \rho_o/\lambda$, and the normalized variance of the random potential is $\sigma'^2 = \mathbb{E}[V'^2] = \lambda^2 \sigma^2$. Note that the relevant parameters are invariant with respect to the normalization, $X_c = X'_c = \sigma'^2/3 \ell_c'/\alpha'^{1/3}$, $X_o = X'_o = \sigma'^2/3 \rho_o'/\alpha'^{1/3}$, and $z_c/\lambda = z'_c = 1/(2\sigma'^2/3\alpha'^{2/3})$.

The normalized paraxial Eq. (H1) is solved using a pseudospectral split-step method, with a frequency cutoff of the
spectral grid $k'_c = 2\pi$ (i.e., $k_c = 2\pi/\lambda = k_o$ refers to the
light wavenumber in dimensional units), so that the spatial
discretization is $dx' = 1/2$ (i.e., $dx = \lambda/2$ in dimensional
units). In all simulations, the size of the spatial window,
$T_{x'}$, is chosen to be much larger than $\ell'_c$. Typically, we
take $T_{x'}/\ell'_c \simeq 40$. Each realization of the random processes
$V'(x', z')$ and $\psi_{z'=0}(x')$ are defined in the spectral domain
using Gaussian correlation functions characterized by $\ell'_c$ and
$\rho'_o$ respectively. The results presented in the main part of
the text are the results of the numerical simulations, aver-

aged over 1000 realizations in cases: 1) of an initial plane
wave and 2) of a coherent speckled field, corresponding to
situation (c). In the case of a partially coherent speckled
initial field, situation (pc), we perform 300 realizations of
the potential $V'(x', z')$. For each of those realizations, we
perform an average over 400 realizations of the initial field
$\psi_{z'=0}(x')$. The different realizations are performed in par-
allel using HPC resources from DNUM CCUB (Centre de
Calcul de l'Université de Bourgogne).

[S1] J.L. Synge, The Hamiltonian method and its application to water waves, Proceedings of the Royal Irish Academy. Section A: Mathematical and Physical Sciences **63**, 1–34 (1963).

[S2] M.V. Berry, Elementary branching: Waves, rays, decoherence, J. Opt. **22**, 115608 (2020).

[S3] L.C. Andrews and R.L. Phillips, Laser Beam Propagation through Random Media, SPIE-International Society for Optical Engineering, 2005.

[S4] J. Garnier and K. Sølna, Coupled paraxial wave equations in random media in the white-noise regime, Ann. Appl. Probab. **19**, 318–346 (2009).

[S5] J. Garnier and K. Sølna, Scaling limits for wave pulse trans-

mission and reflection operators, Wave Motion **46**, 122–143 (2009).

[S6] J. Garnier and K. Sølna, Fourth-moment analysis for beam propagation in the white-noise paraxial regime, Archive on Rational Mechanics and Analysis **220**, 37–81 (2016).

[S7] J.-P. Fouque, J. Garnier, G. Papanicolaou, and K. Sølna, Wave Propagation and Time Reversal in Randomly Layered Media, Springer, 2007.

[S8] L. Isserlis, On a formula for the product moment coefficient of any order of a normal frequency distribution in any number of variables, Biometrika **12**, 134–139 (1918).