WAVE PROPAGATION IN RANDOM MEDIA: BEYOND GAUSSIAN STATISTICS

. In this paper we review some aspects of wave propagation in random media. In the physics literature the picture seems simple: for large propagation distances, the wavefield has Gaussian statistics, mean zero, and second-order moments determined by radiative transfer theory. The results for the first two moments can be proved under general circumstances by multiscale analysis. The Gaussian conjecture for the statistical distribution of the wavefield can be proved in some propagation regimes, such as the white-noise paraxial regime that we address in the first part of this review. It may, however, be wrong in other regimes, such as in randomly perturbed open waveguides, that we address in the second part of this review. In the third and last part, we reconcile the two results by showing that the Gaussian conjecture is restored in randomly perturbed open waveguides in the high-frequency regime, when the number of propagating modes increases.


Introduction
In many wave propagation scenarios the medium is not spatially homogeneous, but may vary on a spatial scale that is small compared to the typical propagation distance.This is the case for wave propagation through the turbulent atmosphere, the Earth's crust, the ocean, and complex biological tissue for instance.If one aims to use transmitted or reflected waves for communication or imaging purposes it is important to characterize how such heterogeneities affect the waves.
Motivated by the situation described above we consider wave propagation through time-independent complex media with a spatially varying index of refraction that is unknown to us, so we model it as a realization of a random process.When the index of refraction is a random process, the wavefield is itself a random process

Multiscale analysis
In its most common form, the analysis of wave propagation in random media consists in studying the wavefield solution of the scalar time-harmonic wave or Helmholtz equation with a randomly heterogeneous index of refraction.Even though the scalar wave equation is simple and linear, the relation between the statistics of the index of refraction and the statistics of the wavefield is highly nontrivial.In order to simplify and understand this relation, one can carry out a multiscale analysis that transforms the random Helmholtz equation into a mathematically tractable yet physically relevant problem.This analysis is based on a separation of scales technique and limit theorems (homogenization, diffusion-approximation, ...), in the framework set forth by G. Papanicolaou and coauthors [2].The wave propagation problem can, indeed, be characterized by several length scales: the typical wavelength (which depends on the source), the correlation radius of the medium, the typical propagation distance.The bandwidth of the source and the relative amplitude of the medium fluctuations can also play a role.Different scaling regimes (corresponding to different physical configurations) can be analyzed when certain ratios between these length scales go to zero or infinity.They lead to tractable and relatively easy to interpret asymptotic results.Typically the solution of the random Helmholtz equation can be shown to converge to the solution of a deterministic partial differential equation or of a stochastic partial differential equation driven by a Brownian field.Stochastic calculus can then be used to compute quantities of interest.
In the random travel time model, which is a special high-frequency regime in which the wavelength is much smaller than the correlation radius of the medium, the fluctuations of the medium affect only the phase of the wave, which satisfies a random eikonal equation [8,56].In the random paraxial model in which the wavelength is smaller than the correlation radius of the medium, backscattering can be neglected but there is significant lateral scattering as the wave advances over long propagation distances and the wavefield satisfies a random Schrödinger-type equation [55,58].This is the regime that we address in Section 2. In the randomly layered regime, in which the medium is only varying along the longitudinal direction (along the propagation axis), there is significant backscattering and the plane wave mode amplitudes satisfy a system of ordinary random differential equations [22].In the radiative transfer regime, in which the wavelength is of the same order as the correlation radius of the medium, the angular distribution of the mean wave energy satisfies a transport equation [51,52].In the waveguide regime, in which the wave is -partially-trapped transversely by a special profile of the index of refraction, the wavefield can be decomposed as a superposition of wave modes and the wave mode amplitudes satisfy coupled stochastic differential equations [22,26].This is the regime that we address in Section 3.
In this review paper we first consider a scaling regime corresponding to long-range high-frequency beam propagation and small-scale medium fluctuations giving negligible backscattering.This is the so-called whitenoise paraxial regime, that gives rise to the Itô-Schrödinger model, which is presented in Section 2.1.This model is a simplification of the random Helmholtz equation since it corresponds to an evolution problem, but yet in the regime that we consider it describes the propagated field in a weak sense in that it gives the correct statistical structure of the wavefield.The Itô-Schrödinger model can be derived rigorously from the random Helmholtz equation by a separation of scales technique in the high-frequency regime (see [3] in the case of a randomly layered medium and [29][30][31] in the case of a three-dimensional random medium).It models many situations, for instance laser beam propagation [53], underwater acoustics [55], or migration problems in geophysics [12].The Itô-Schrödinger model makes it possible to use Itô's stochastic calculus, which in turn enables the closure of the hierarchy of moment equations [23,38].The equations for the first-order and second-order moments are easy to solve.The equation for the fourth-order moments is difficult and only approximations or numerical solutions have been available for a long time (see [20,57,59] and [38,Sec. 20.18]).In a special scaling regime, it is, however, possible to derive expressions for the fourth-order moments as presented in Section 2.2.
In the second part of this review paper we address wave propagation in a randomly perturbed open waveguide.This is a physically relevant model as it can be encountered in optics (dielectric slab waveguide [47]) or in underwater acoutics (Pekeris waveguide [61]) for instance.The wavefield can then be decomposed as the superposition of propagating, radiating, and evanescent modes, and the mode amplitudes satisfy ordinary differential equations that are coupled by the perturbations of the index of refraction or the waveguide geometry (Section 3.2).When the wavelength is smaller than or of the order of the diameter of the waveguide and the correlation length of the perturbation, which are themselves much smaller than the propagation distance along the waveguide axis, then it can be shown that the mode amplitudes satisfy stochastic differential equations driven by correlated Brownian motions as shown in Section 3.3.Using this model it is possible to write closed equations for the moments of the wave mode amplitudes.The analysis of these equations reveal non-classical forms of equipartition of energy as shown in Section 3.5 and exponential growth of the relative variances of the mode amplitudes as shown in Section 3.6.

The Gaussian conjecture
Star scintillation is a well-known paradigm, related to the observation that the irradiance of a star fluctuates due to interaction of the light with the turbulent atmosphere.Experimental observations indicate that the statistical distribution of the irradiance is exponential, with the irradiance being the square magnitude of the complex wavefield.In the physics literature it is a well-accepted conjecture that the statistics of the complex wavefield becomes circularly symmetric complex Gaussian when the wave propagates through the turbulent atmosphere [60,62], so that the irradiance is the sum of the squares of two independent real Gaussian random variables, which has chi-square distribution with two degrees of freedom, that is an exponential distribution.The mathematical proof of this conjecture has been obtained in randomly layered media [22,Chapter 9] but is still incomplete in three-dimensional random media [23,32,33,41].On the one hand, in Section 2.3 we report results for the fourth-order moments that are consistent with the Gaussian conjecture in the paraxial white noise regime.On the other hand, in Sections 3.6-3.7 we study the fluctuations of the mode amplitudes of a wave propagating in a randomly perturbed open waveguide and we show that the relative variances of the mode powers grow exponentially with the propagation distance.This demonstrates that the Gaussian conjecture is wrong in this regime and that the picture is not simple.
Certain functionals of the wavefield carry information about the medium and can be characterized in some specific regimes [4,5,19,50].For instance, the Wigner distribution (a transform in time-frequency analysis) is known to be convenient to study the solution of Schrödinger equation [34,51].An important issue is the socalled statistical stability property: we look for functionals that become deterministic in the considered scaling regime and that depend only on the statistics of the random medium and not on the particular realization.The conjecture is that this can happen for well chosen functionals in the limit of rapid decorrelation of the medium fluctuations.In [42,43] the authors consider a situation with rapidly fluctuating random medium fluctuations in which the Wigner distribution is statistically stable.As shown in [4], however, the statistical stability also depends on the initial data and can be lost for very rough initial data even with a high lateral diversity as considered there.In Section 2.4 we present a detailed and quantitative analysis of the stability of the Wigner distribution in the white-noise paraxial regime.We derive an explicit expression of the coefficient of variation of the smoothed Wigner distribution as a function of the smoothing parameters, in the general situation in which the standard deviation can be of the same order as the mean.This is a realistic scenario, which is not too deep into a statistical stabilization situation, but which gives partly coherent but fluctuating wave functionals.These results make it possible to quantify such fluctuations and how their magnitudes can be controlled by optimal smoothing of the Wigner distribution.

Wave propagation in an open random medium
In this section we address wave propagation in a d + 1-dimensional random medium and we describe how to derive the mathematically tractable Itô-Schrödinger model from the wave equation in a random medium.We consider the d + 1-dimensional scalar wave equation: Here the source emits a time-harmonic signal with frequency ω and it is localized in the plane z = 0: and the speed of propagation is spatially heterogeneous where µ is a zero-mean random process with stationary and ergodic properties.

The white-noise paraxial model
The time-harmonic field p such that p(t, ⃗ x) = p(⃗ x)e −iωt is solution of the random Helmholtz equation where ∆ ⃗ x = ∆ ⊥ + ∂ 2 z .Note that, when µ ≡ 0 and f s ≡ 1, then the solution is the plane wave p(⃗ x) = ico 2ω e i ωz co .When the initial source f s is spatially limited (belongs to L 2 (R d )) and µ is a stationary and ergodic process, the function φ (slowly-varying envelope of a plane wave going along the z-axis) defined by satisfies the modified Helmholtz equation In the white-noise paraxial regime "λ ≪ ℓ c , r o ≪ z" (which means, the wavelength λ = 2πc o /ω is much smaller than the correlation radius ℓ c of the medium and the radius r o of the source, which are themselves much smaller than the propagation distance z) the forward-scattering approximation in direction z is valid and φ satisfies the Itô-Schrödinger equation [30] where • stands for the Stratonovich integral, B(x, z) is a Brownian field, that is a Gaussian process with mean zero and covariance Remark.Existence, uniqueness and continuity of solutions of the Itô-Schrödinger model (6) are established in [16].The proof of the convergence of the solution to (5) to the solution to (6) is in [30].A first proof in the case of layered media (i.e. when µ ≡ µ(z)) can be found in [3].More precisely, the white-noise paraxial regime "λ ≪ ℓ c , r o ≪ z" corresponds to the scaled regime where ε is a small dimensionless parameter.The scaled version of (5) for φε (x, z) = φ( x ε 2 , z) is then When ε → 0 the second-order derivative in z is negligible (which proves the forward scattering approximation used in the physics literature) and the potential behaves like a white noise in z.The convergence in ) of the solution of the scaled version (8) to the solution of the Itô-Schrödinger equation (6) can then be established in distribution when ε → 0. This result requires mixing properties for the random process µ.If, however, the process µ has long-range properties, in the sense that as |z| → +∞ and α ∈ (0, 1), then, under appropriate technical and scaling assumptions [37], the limiting equation is the fractional Itô-Schrödinger model (6) in which B is a fractional Brownian field with Hurst index H = 1 − α/2 ∈ (1/2, 1), i.e. a Gaussian field with mean zero and covariance In this case the stochastic integral in (6) can be understood as a generalized Stieltjes integral (as H > 1/2).

Statistics of the wavefield
In this section we describe how to compute the moments of the wavefield.By Itô's formula and (6), the coherent (or mean) wave satisfies the Schrödinger equation with homogeneous damping (for z > 0): and therefore E φ(x, z) = φhom (x, z) exp(−z/ℓ sca ), where φhom is the solution in the homogeneous medium.The coherent wave amplitude decays exponentially with the propagation distance and the characteristic decay length is the scattering mean free path ℓ sca : This result shows that any coherent imaging or communication method fails in random media when the propagation distance is larger than the scattering mean free path [27].But the wave energy is not lost or dissipated, it is converted into incoherent (zero-mean) fluctuations that can be studied.
The mean Wigner distribution defined by is the local Fourier transform of the second-order moment of the wavefield (the bar stands for complex conjugation).It is the angularly-resolved mean wave energy density.By Itô's formula and (6), it solves the radiative transfer equation starting from W m (x, ξ, z = 0) = W s (x, ξ), the Wigner distribution of the source or initial field f s .Here is the Fourier transform of γ and determines the scattering cross section of the radiative transfer equation.This result shows that the fields observed at nearby points are correlated and their correlations contain information about the medium.Accordingly, one should use local cross correlations for imaging and communication in random media [8,10].
In order to quantify the stability of correlation-based imaging methods and to validate the Gaussian conjecture, one needs to evaluate variances of empirical correlations, which involves the fourth-order moment: By Itô's formula and (6), it satisfies the Schrödinger-type equation with the generalized potential These moment equations have been known for a long time [38].Recently it was shown [33] that in the regime "λ ≪ ℓ c ≪ r o ≪ z" the fourth-order moment can be expressed explicitly in terms of the function γ.These results can be used to address a wide range of applications in imaging and communication.They also clarify the Gaussian conjecture used in the physics literature as we will see below.

The scintillation index
In this section we study the intensity fluctuations of the wavefield solution of (6) and characterize the scintillation index which quantifies the relative intensity fluctuations.It is a fundamental quantity associated for instance with light propagation through the atmosphere [38].It is defined as the square coefficient of variation of the intensity [38,Eq. (20.151)]: Note that, if φ had (complex) Gaussian statistics, then the scintillation index would be equal to one.
When the spatial profile of the source (or the initial beam) has a Gaussian profile, and when "λ ≪ ℓ c ≪ r o ≪ z", the behavior of the scintillation index can be described as follows [33].
Proposition 2.1.Let us consider the following form of the covariance function of the medium fluctuations: with γ(0) = 1 and the width of the function γ is of order one.In the regime "λ ≪ ℓ c ≪ r o ≪ z" the scintillation index ( 16) has the following expression: In fact this result follows from the complete expressions of the second moment of the intensity and the second moment of the field that are given in [33].The scintillation index at the beam center x = 0 is a function of z/ℓ sca and z/ℓ dif only, where ℓ dif = ωr o ℓ c /c o is the typical propagation distance for which diffractive effects are of order one, as shown in [30,Eq. (4.4)].It is interesting to note that, even if the propagation distance is larger than the scattering mean free path, the scintillation index can be smaller than one if ℓ dif is small compared to ℓ sca .
In order to get more explicit expressions that facilitate interpretation of the results let us assume that γ(x) is smooth and can be expanded as When scattering is strong in the sense that the propagation distance is larger than the scattering mean free path z ≫ ℓ sca , the expressions of the second moments of the field and of the intensity can be simplified: where the beam radius R z is and the correlation radius of the beam ρ z is Note that, in this regime, the fourth-order moments satisfy the Isserlis formula (i.e. they can be expressed in terms of sums of products of second-order moments), and therefore the scintillation index S(x, z) is equal to one.This observation is consistent with the physical conjecture that, in the strongly scattering regime z/ℓ sca ≫ 1, the wavefield should have zero-mean complex circularly symmetric Gaussian statistics, and therefore the intensity should have exponential (or Rayleigh) distribution [20,38].

Fluctuations of the Wigner distribution
In this section we give an explicit characterization of the signal-to-noise ratio of the Wigner distribution.The Wigner distribution of the wavefield is defined by It can be interpreted as the angularly-resolved wave energy density (note, however, that it is real-valued but not always non-negative valued).Its expectation W m defined by (11) satisfies the radiative transfer equation (12).It is known that the Wigner distribution is not statistically stable, and that it is necessary to smooth it (that is to say, to convolve it with a smoothing kernel) to get a quantity that can be measured in a statistically stable way (that is to say, the smoothed Wigner distribution for one typical realization is approximately equal to its expected value) [4,50].Our goal in this section is to quantify this statistical stability.
Let us consider two positive parameters r s and ξ s and define the smoothed Wigner distribution: If we denote by ρ z the correlation radius of the field (given by (21) in the strongly scattering regime), we may anticipate that r s and 1/ξ s should be of the order of ρ z to ensure averaging.The coefficient of variation C s of the smoothed Wigner distribution, which characterizes its statistical stability, is defined by: An exact expression of the coefficient of variation of the smoothed Wigner distribution can be derived in the regime "λ ≪ ℓ c ≪ r o ≪ z" [33].It involves four-dimensional integrals and it is, therefore, difficult to derive its qualitative behavior.This expression becomes simple in the strongly scattering regime z ≫ ℓ sca .We then get the following expression for the coefficient of variation [33].
Proposition 2.2.In the regime "λ ≪ ℓ c ≪ r o ≪ z", if additionally z ≫ ℓ sca and γ can be expanded as (19), then the coefficient of variation of the smoothed Wigner distribution (23) satisfies: where ρ z is the correlation radius (21).
Note that the expression (25) of the coefficient of variation is independent of x and ξ.It is a formula that is simple enough to help determining the smoothing parameters ξ s and r s that are needed to reach a given value for the coefficient of variation: • For 2ξ s r s < 1 (resp.> 1) we have C s (x, ξ, z) > 1 (resp.< 1); in other words, the smoothed Wigner transform can be considered as statistically stable when 2ξ s r s > 1.
The critical value r s = 1/(2ξ s ) is indeed special.In this case, the smoothed Wigner distribution ( 23) can be written as the double convolution of the Wigner distribution W of the random field φ(•, z) (defined by ( 22)) with the Wigner distribution of the Gaussian state φg since we have , and therefore for r s = 1/(2ξ s ).It is known that the convolution of a Wigner distribution with a kernel that is itself the Wigner distribution of a function (such as W g ) is nonnegative real valued (the smoothed Wigner distribution obtained with the Gaussian W g is called Husimi function) [11,45].This can be shown easily in our case as the smoothed Wigner distribution can be written as for r s = 1/(2ξ s ).From this representation formula of W s valid for r s = 1/(2ξ s ), we can see that it is the square modulus of a linear functional of φ(•, z).The physical conjecture that φ(•, z) has circularly symmetric complex Gaussian statistics in strongly scattering media then predicts that W s (x, ξ, z) should have an exponential distribution, because the sum of the squares of two independent real-valued Gaussian random variables has an exponential distribution.This is indeed consistent with Proposition 2.2 that shows that C s = 1 for r s = 1/(2ξ s ).
If r s > 1/(2ξ s ), the smoothed Wigner distribution ( 23) can be expressed as: where the function Ψ is defined by From this representation formula for W s valid for r s > 1/(2ξ s ), we can see that it is nonnegative valued and that it is a local average of ( 29), which has a unit coefficient of variation in the strongly scattering regime.That is why the coefficient of variation of the smoothed Wigner distribution is smaller than one when r s > 1/(2ξ s ).

Wave propagation in a random waveguide
In this section we address wave propagation in a two-dimensional randomly perturbed open waveguide.Our model consists of a two-dimensional waveguide with range axis denoted by z ∈ R and transverse coordinate denoted by x ∈ R (see Figure 1).A point-like source at a fixed position (x, z) = (x s , 0) transmits a time-harmonic signal.The time-harmonic wavefield p(x, z) satisfies the Helmholtz equation: where k is the homogeneous wavenumber and n(x, z) is the index of refraction at position (x, z).
In the case of ideal (unperturbed) waveguides, the index of refraction is range-independent and equal to where n > 1 is the relative index of the core and d > 0 is its diameter.We are interested in randomly perturbed waveguides.We can address two types of random waveguides.Type I perturbation: in the first type, the index of refraction within the core region x ∈ (−d/2, d/2) is randomly perturbed [7,13,14,35,40]: The fluctuations are modeled by the zero-mean, bounded, stationary in z random process ν(x, z) with smooth covariance function It satisfies strong mixing conditions in z as defined for example in [48, section 2].The typical amplitude of the fluctuations of index of refraction is assumed to be much smaller than 1 and it is modeled by the small and positive dimensionless parameter ε.
Type II perturbation: in the second type (see Figure 1), the boundaries of the core are randomly perturbed [1,9,36,46,47]: + (z) and z ∈ (0, where The fluctuations are modeled by the zero-mean, bounded, independent and identically distributed stationary random processes ν 1 and ν 2 with smooth covariance function They satisfy strong mixing conditions.The typical amplitude of the fluctuations of the boundaries is assumed to be much smaller than the core diameter d and it is modeled in (37) by the small and positive dimensionless parameter ε.
We study the wavefield at z > 0, satisfying p(x, z) ∈ C 0 (0, +∞), H 2 (R) ∩ C 2 (0, +∞), L 2 (R) , and to set radiation conditions, we have assumed that the random fluctuations are supported in the range interval (0, L (ε) ).We will see that net scattering effect of these fluctuations becomes of order one at range distances of order ε −2 , so we consider the interesting case L (ε) = L/ε 2 .

Homogeneous waveguide
In this section, we consider an index of refraction of the form (33), which is stepwise constant.There is no fluctuation of the medium along the z-axis.The analysis of the perfect waveguide is classical [44,61].The Helmholtz operator has a spectrum of the form where the N = [ √ n 2 − 1kd] modal wavenumbers β j are positive and The generalized eigenfunctions ϕ t,γ , t ∈ {e, o}, associated to the spectral parameter γ in the continuous spectrum (−∞, k 2 ) and the eigenfunctions ϕ j , j = 0, . . ., N − 1, associated to the discrete spectrum, are given in [25].The generalized eigenfunctions ϕ e,γ are even and ϕ o,γ are odd.The eigenfunctions ϕ j are even for even j and odd for odd j.Any function can be expanded on the complete set of the eigenfunctions of the Helmholtz operator.In particular, any solution of the Helmholtz equation in homogeneous medium can be expanded as The modes for j = 0, . . ., N − 1 are guided, the modes for γ ∈ (0, k 2 ) are radiating, the modes for γ ∈ (−∞, 0) are evanescent.Indeed, the complex mode amplitudes satisfy for any z ̸ = 0. Therefore, if the source is of the form (32), we have for z > 0: where the mode amplitudes are constant and determined by the source:

Random waveguide
We consider the two types of random perturbations described in Section 3. In both cases we can write where the perturbation is of the form for type I perturbations, and for type II perturbations.The solution of the perturbed Helmholtz equation (32) can be expanded as ( 40) and the complex mode amplitudes satisfy the coupled equations for z ∈ (0, L (ε) ): for j = 0, . . ., N − 1, for γ ∈ (−∞, k 2 ) and t ∈ {e, o}, with and (•, •) L 2 stands for the standard scalar product in L 2 (R).These equations are obtained by substituting the ansatz ( 40) into (32) and by projecting onto the eigenmodes.
From the definitions ( 46) or (47) of V (ε) (x, z) and the Taylor expansions of the eigenfunctions ϕ j (x) and ϕ t,γ (x) around x = ±d/2, we can obtain power series (in ε) expressions of the coefficients C (ε) j,l : and similarly for C t,γ,l , and C (ε) t,γ,t ′ ,γ ′ .We finally introduce the generalized forward-going and backward-going mode amplitudes: for t ∈ {e, o}, which are defined such that and We can substitute (56-57) into (48)(49) in order to obtain the first-order system of coupled random differential equations satisfied by the mode amplitudes ( 55): with similar equations for b j and b t,γ .This system is complemented with the boundary conditions at z = 0 and z = L (ε) : , where a j,s and a t,γ,s are defined by (44)(45).The evanescent mode amplitudes pt,γ , t ∈ {e, o}, γ ∈ (−∞, 0), satisfy (49).

The effective Markovian dynamics for the mode amplitudes
We rename the complex mode amplitudes in the long-range scaling as We can follow the lines of [35] to get the following results.
(1) In the regime ε ≪ 1 the evanescent mode amplitudes, that satisfy (49), can be expressed to leading order in closed forms as functions of the guided and radiating mode amplitudes (60)(61).Indeed it is possible to invert the operator ∂ 2 z + γ in (49) for γ < 0 by using the Green's function that satisfies the radiation condition and to obtain: for z > 0, γ < 0 and t ∈ {e, o}.Here we recognize ) Under the assumption that the power spectral density R I (κ, x, x ′ ) for type-I perturbations (or R II (κ) for type-II perturbations) has compact support or fast decay, the forward-scattering approximation can be proved, i.e. the coupling between forward-going and backward-going mode amplitudes is negligible, so that we have (3) The forward-going guided mode amplitudes (a ε j ) N −1 j=0 and radiating mode amplitudes (a ε t,γ ) γ∈(0,k 2 ),t∈{e,o} then satisfy a closed linear system of the form with initial conditions for a ε at z = 0.Here F, resp.G, is an operator with zero mean, resp.non-zero mean, and ergodic properties inherited from those of the processes ν.
We can finally apply a diffusion approximation theorem to establish the following result (see [25,35] for the full statement, [40] for a first version in which the contributions of the evanescent modes is neglected, which means that the operator L 3 is missing in the expression of the generator L, and [25] for a detailed analysis).
Proposition 3.1.The random process ), the space of continuous functions from [0, L] to C N × L 2 ((0, k 2 )) 2 , to the Markov process (a j (z)) N −1 j=0 , (a t,γ (z)) γ∈(0,k 2 ),t∈{e,o} with infinitesimal generator L.Here C N × L 2 ((0, k 2 )) 2 is equipped with the weak topology and the infinitesimal generator has the form L = L 1 + L 2 + L 3 , where L j , 1 ≤ j ≤ 3, are the differential operators: In these definitions we use the classical complex derivative: , and the coefficients of the operators (63-65) are defined for j, l = 0, . . ., N − 1, as follows: -For all j ̸ = l, Γ jl and Γ s jl are given by with R jl (z) defined by -For all j, l: -For all j, Λ j is defined by and where R j,t,γ (z) = E[C j,t,γ (0)C j,t,γ (z)] is defined as in (68) upon substitution (t, γ) for l and We give some remarks before focusing our attention on the mode powers.
(1) The convergence result holds in the weak topology.This means that we can only compute quantities of the form E[F (a 0 , . . ., a N −1 , 0 α e,γ a e,γ dγ, 0 α e,γ a ε e,γ dγ, (2) The generator L does not involve ∂ at,γ nor ∂ āt,γ .Therefore (a ε j (z)) N −1 j=0 converges in distribution in C 0 ([0, L], C N ) to the Markov process (a j (z)) N −1 j=0 with generator L. The weak and strong topologies are the same in C N , so we can compute any moment of the form E[F (a 0 , . . ., a N −1 )], which are the limits of E[F (a ε 0 , . . ., a ε N −1 )].(3) L 1 is the contribution of the coupling between guided modes, which gives rise to power exchange between the guided modes (effective diffusion).(4) L 2 is the contribution of the coupling between guided and radiating modes, which gives rise to power leakage from the guided modes to the radiating ones (effective attenuation) and addition of frequencydependent phases on the guided mode amplitudes (effective dispersion).The effective attenuation and dispersion are produced by causal phenomena and they are related to each other through Kramers-Konig relations [28].(5) L 3 is the contribution of the coupling between guided and evanescent modes, which gives rise to additional phase terms on the guided mode amplitudes (effective dispersion).This term is the main effect when the waveguide supports only one propagating mode and the core boundaries are hard or soft so that there is no radiating mode [24].(6) If the generator L is applied to a test function that depends only on the mode powers (P j ) N −1 j=0 , with P j = |a j | 2 , then the result is a function that depends only on (P j ) N −1 j=0 .Thus, the mode powers (P j (z)) N −1 j=0 define a Markov process, with infinitesimal generator defined by (72) below.(7) The radiation mode amplitudes remain constant in L 2 ((0, k 2 )) 2 , equipped with the weak topology, as ε → 0. However, this does not describe the power t∈{e,o} k 2 0 |a ε t,γ | 2 dγ transported by the radiation modes, because the convergence does not hold in the strong topology of L 2 ((0, k 2 )) 2 so we do not have t∈{e,o} (8) When N = 1, then the generator is This shows that a 0 (the amplitude of the unique guided mode) has the same distribution as where W 1 z is a standard Brownian motion.The mode amplitude experiences a random phase modulation and a deterministic damping, which both depend on frequency and two-point statistics of the medium perturbations [24].(9) When N ≥ 2, the limit process (a j (z)) N −1 j=0 can be identified as the solution of a system of stochastic differential equations driven correlated Brownian motions, as formulated in the following corollary.
Then the limit Markov process (a j (z)) N −1 j=0 has the same distribution as the unique solution of starting from a j (z = 0) = a j,s , j = 0, . . ., N − 1, where • stands for the Stratonovich integral.
The proof of the corollary is a straightforward application of Itô's formula.

The effective Markovian dynamics for the mode powers
From the result for the complex mode amplitudes we get the following result.
The coefficients Γ jl describe the effective mode coupling between guided modes due to random scattering.The coefficients Λ j are effective mode-dependent dissipation coefficients and they come from the coupling between guided and radiative modes due to random scattering.
From the form of the generator L P , we can establish that the nth-order moments of the mode powers satisfy closed equations.We will apply this to compute the first moments of P , as well as its second moments later in Section 3.6.
Using (72) we find that the mean mode powers satisfy the closed system of equations starting from Q j (0) = |a j,s | 2 .This equation looks like a discrete radiative transfer equation, with the matrix Γ that plays the role of the scattering cross section.The form of these coupled-mode equations is well-known [17] although the mode-dependent attenuation terms Λ j are usually introduced heuristically.The solution explicitly writes: with the matrix A defined by (δ jl is the Kronecker symbol): We can also remark that the system (74) can be interpreted as the Kolmogorov equation associated to a random walk on the finite space {0, . . ., N − 1}.If we denote by (J z ) z≥0 the Markov process on the state space {0, . . ., N − 1} with infinitesimal generator Γ, then a Feynman-Kac formula gives the following probabilistic representation of the mean mode powers Q j (z): This representation makes it possible to anticipate the result derived below in the high-frequency regime (when the number of modes becomes large), namely that the Q j 's can be approximated by the solution of a diffusion equation, because the normalized random walk (J z /N ) z≥0 can be approximated in distribution by a diffusion process on [0, 1].

Long-range behavior of the mean mode powers
From now on we assume that the symmetric matrix Γ defined by Γ jl given by (66) for j ̸ = l, Γ jj = − l ′ ̸ =j Γ jl ′ , is irreducible.By Perron-Frobenius theorem, the first eigenvalue of A is simple and nonpositive (we denote it by −λ) and the components of the corresponding unit eigenvector V have all the same sign (so we can assume that they are positive).By (75) we get the following result.
Proposition 3.2.The mean mode powers (73) satisfy where (−λ, V ) is the first eigenvalue/eigenvector of A and c V = N −1 l=0 V l |a l,s | 2 .Special cases with explicit expressions are given in [25].The main results are the following ones: (1) If there is no dissipation (i.e.Λ j = 0 for all j), then A = Γ and the first eigenvalue/eigenvector (−λ (0) , V (0) ) of the matrix Γ is λ , which gives the standard equipartition result [15,22,26]: The total input power N −1 l=0 |a l,s | 2 becomes equipartitioned amongst all guided modes.(2) When there is dissipation (i.e.there exists j such that Λ j > 0), then there is effective damping λ > 0 and the first eigenvector V is not the equipartitioned vector V (0) anymore.More exactly: (a) When coupling (described by the matrix Γ) is stronger than dissipation (described by the matrix Φ with Φ jl = Λ j δ jl ), then the effective damping coefficient λ of the mean mode powers is approximately the arithmetic average of the effective mode-dependent damping coefficients Λ j , and the distribution of the mean mode powers is approximately equipartitioned, with slightly reduced allocations for the modes with the strongest damping coefficients.(b) When coupling is weaker than dissipation, then the effective damping of the mean mode powers is approximately the minimum of the effective mode-dependent damping coefficients Λ j .The distribution of the mean mode powers is, moreover, essentially concentrated on the modes corresponding to the minimum of the mode-dependent damping coefficients.

Fluctuation analysis
By (72) we find that the second-order moments of the mode powers satisfy the closed equations This system has the same form as the one found in the literature dedicated to coupled mode theory [15,17].The initial conditions are R jl (0) = |a j,s | 2 |a l,s | 2 .Let us introduce S = (S jl ) 0≤j≤l≤N −1 defined by The S jl 's satisfy the system with the convention that whenever S jl occurs with j > l, it is replaced by S lj .The operator Θ is the infinitesimal generator of a random Markov process (J z , L z ) z≥0 that is a random walk on the discrete triangle {(j, l) ∈ N 2 , 0 ≤ j ≤ l ≤ N − 1}.Using a Feynmac-Kac formula we get the following probabilistic representation of S jl : We can anticipate that, in the continuum limit, the Markov process (J z /N, L z /N ) z≥0 behaves as a diffusion process on the triangle {(u, v) ∈ R 2 , 0 ≤ u ≤ v ≤ 1}, and, therefore, S jl satisfies a diffusion equation on this triangle.
Without dissipation we have the following result for the relative fluctuations of the pointwise intensity: that grows exponentially with the propagation distance.This proves that the Gaussian conjecture is wrong for wave propagation in random waveguides in the presence of radiative damping.The explicit formulas given in [25] for some particular situations show that the growth rate of the relative variance of the intensity increases when the effective modal dissipation coefficients become different from each other and decreases when the number of modes increases.The analysis in the high-frequency regime carried out in the next section confirms that the exponential growth rate of the relative intensity fluctuations vanishes when the number of modes goes to infinity.

The high-frequency regime for wave propagation in random waveguides
In Section 2 we have addressed wave propagation in open random media in the regime λ ≪ r o , ℓ c ≪ L (with r o the radius of the initial beam) and we have found that the Gaussian conjecture could be validated.In Section 3 we have addressed wave propagation in random waveguides in the regime λ ≲ d, ℓ c ≪ L (with d the diameter of the waveguide) and we have found that the Gaussian conjecture was wrong as the relative fluctuations of the intensity are growing exponentially with the propagation distance.In this section we address the high-frequency regime of wave propagation in random waveguides, that is to say, the regime λ ≪ d, ℓ c ≪ L, and we show that we can reconcile the two results.
We consider wave propagation in random waveguides as introduced in Section 3. The following arguments show that, when λ ≪ d, then the number of modes is large, and when λ ≪ ℓ c , then coupling between guided modes is essentially via nearest neighbors (between modes with close modal wavenumbers): 1) The number of modes becomes large when (n 2 − 1)k 2 d 2 ≫ 1.In other words, the number of modes is large when the frequency is large because it is proportional to the ratio of the waveguide diameter over the wavelength.
2) For type II perturbations, coupling via nearest neighbors happens when the fluctuations of the boundaries are smooth so that the Fourier transform of R II decays fast and the correlation radius is much larger than the wavelength.Under such circumstances, we have β j − β j+l ≃ π For type I perturbations, coupling via nearest neighbors happens under similar conditions.Indeed, let us assume that R I (x, x ′ , z) can be factorized as R I (x, x ′ , z) = R I,c (x, x ′ )R I,l (z),

Stability of the intensity fluctuations
In the high-frequency regime, when the number of modes becomes large, we have µ = 2λ and (109) holds.Therefore there is no exponential growth of the fluctuations and we have which corresponds to a relative variance (or scintillation index) equal to one.We recover the standard result that the wavefield, in the point of view of the fourth-order moments, behaves as a Gaussian process with relative variance (scintillation index) equal to one [33].

Conclusion
In this paper we have reviewed an asymptotic theory of wave propagation in random media and we have contrasted the situations in open media and waveguides.We have presented standard results about the first two moments of the wave amplitudes: the mean amplitudes decay exponentially with the propagation distance and the mean Wigner transform (for the open medium case) and the mean mode powers (for the waveguide case) satisfy a radiative transfer equation.The fourth-order moment analysis reveals that the fluctuations of the wave amplitudes may have very different behaviors.The wave amplitudes have asymptotically Gaussian statistics when scattering is strong enough in the open medium case in the white-noise paraxial regime, while the variance of the fluctuations may grow exponentially with the propagation distance in the waveguide case.We have carefully studied the exponential growth rates of the relative variances for different models of randomly perturbed waveguides.We have shown that, in the high-frequency regime, when the number of guided modes increases, the exponential growth rates vanish and the scintillation index (the relative variance of the intensity fluctuations) becomes equal to one, as observed in open media in the white-noise paraxial regime.
These results show that incoherent imaging in a random waveguide (such as a Pekeris waveguide in underwater acoustics) is challenging.Indeed incoherent imaging is based on the use of the cross correlations of the recorded signals [18].The estimation of the second-order moments of the wavefield is, however, extremely difficult because of the large variances of the empirical second-order moments and one may need to average over a lot a samples (while the medium may be not stationary as in underwater acoustics).This is in contrast with the situation in open random media where smoothed Wigner transforms are statistically stable [5,6,33].More generally, the results on the fourth-order moments in the waveguide case show that the predictions of the coupled mode equations (which describe the evolutions of the statistical second-order moments of the wavefield, such as Eq.(74)) may be not easy to exploit experimentally.
Finally, we can now challenge the Gaussian conjecture for wave propagation in random media.This conjecture claims that, after a large propagation distance, the wavefield has Gaussian statistics and therefore we know everything as soon as the first two moments are characterized.We have shown in this paper that the conjecture is sometimes right, but it can be wrong; and then the picture given by the first two moments is not complete.

Figure 1 .
Figure 1.Left: An ideal two-dimensional waveguide.Right: A two-dimensional waveguide with cross-section perturbed by random fluctuations of the top and bottom boundaries.The point source is in the plane z = 0.