Extreme Events through Prevailing Backscattering and Their Suppression by a Focusing Nonlinearity

Vincent H. Schultheiss, Martin Wimmer, Stefan Malzer, and Ulf Peschel Institute of Optics, Information and Photonics, Friedrich-Alexander University Erlangen-Nürnberg (FAU), Haberstraße 9a, 91058 Erlangen, Germany Institute of Condensed Matter Theory and Solid State Optics, Abbe Center of Photonics, Friedrich-Schiller-University Jena, Max-Wien-Platz 1, 07743 Jena, Germany Lehrstuhl für Angewandte Physik, Friedrich-Alexander University Erlangen-Nürnberg (FAU), Staudtstraße 7, 91058 Erlangen, Germany


I. INTRODUCTION
The rare, but inevitable occurrences of high-amplitude events, which significantly peak out of the ensemble average, are a widely experienced side effect of wave dynamics in complex systems [1,2].As statistical extremes they are characterized by their ferocity and suddenness and are thus of vital relevance particularly in the context of damage prevention.The most spectacular and devastating incidents happen in the context of hydrodynamics, namely, ocean wave dynamics, where they are labeled rogue waves due to their unpredictable nature [3,4].But even so, reports from such freak waves were disbelieved by the scientific community for a long time.While being plausible as outliers from a statistical point of view, they were deemed too seldom to be witnessed by humans and, hence, not worthy of any methodical investigation.
The singular giant recorded impact that hit the Draupner platform on New Year's Day in 1995 is widely considered as one of the first scientific proofs that rogue waves actually exist and are not simply the plot point of an old seaman's yarn [5].Since then, oil rig and radar-based measurements have confirmed their omnipresence in the open sea [3].It was found that their frequency exceeds expectations by far.As of now, they are identified as a substantial threat to the commercial shipping industry and have stirred up a lot of research as to their origin, formation, and predictability [6].
The astonishment over rogue waves, i.e., waves with heights exceeding twice the so-called significant wave height H 1=3 (the mean of the highest third of all wave heights), mostly arises because of the conception that ordinary ocean waves are very well approximated by a linearization of the underlying set of hydrodynamic equations and are thus to a fair degree sinelike [7,8].That is why a superposition of random amplitudes and phases is expected to follow the central limit theorem yielding the statistics of wave heights to show characteristics of a Rayleigh distribution [9,10].Taking the square of the wave heights, one should find an exponentially decaying distribution in close correspondence with the probability density function p for intensities I of speckle patterns known from illumination with coherent light after its phase front was randomly disturbed, for example, by the reflection off a rough surface [11]: Here, I 0 is the mean intensity over the considered interval.In this vein, the frequent occurrence of rogue waves or, more generally speaking, any long-tailed deviation of the probability density function from a Rayleigh distribution, remains a mystery.There is an ongoing debate whether the origin of this phenomenon is merely linear as, e.g., a dispersion effect [1] or the formation of caustics by spatial focusing on random impurities [12][13][14][15][16], or whether its explanation requires the inclusion of nonlinear terms into the evolution equations [17][18][19][20][21][22].In the past few years, the systematic investigation of extreme wave events has been aided by the close analytic analogy between hydrodynamic surface gravity waves and the envelope of electromagnetic field oscillations [23,24].It has been shown that growing modulation instabilities can create spontaneous solitary peaks in optical systems during supercontinuum generation in glass fibers [25].Since then, many examples for the occurrence of high-amplitude events were found in optics [26][27][28][29].However, it is widely believed that rogue waves are not based on a single physical mechanism, apart from some sort of disturbance in the system to begin with [30].
In this article, we present evidence that linear wave dynamics can generate histograms with pronounced tails, which are underestimated by the central limit theorem, solely on the basis of weak backscattering.To this end we investigate light propagation in a randomly perturbed planar waveguide structure and consider the intensity profile at the rear facet for low and high input beam powers.In contrast to common perception, we find that nonlinearity (NL) suppresses the formation of rogue waves, and that for increasing power, extreme events become more and more unlikely, thus transforming the light field back into a speckle distribution.Our experimental findings are supported by results obtained from discrete numerical simulations based on a random distribution of scattering matrices [31,32].
This article is structured as follows.After providing details of the used waveguide samples and the experimental setup in Sec.II, experimental findings in the linear regime are presented in Sec.III, accompanied by first simulation results that emerge from our numerical model.In Sec.IV, we give an overview of the 2D diffusion approach used to characterize the samples and necessary to understand the numerical model, which is presented in detail in Sec.V.In Sec.VI, we experimentally and numerically investigate the impact of a Kerr nonlinearity on the intensity distribution.Finally, in Sec.VII, we conclude the article with a short summary of our findings.

II. EXPERIMENTAL SETUP
The waveguides used in the experiments consist of layered stacks of Al x Ga x−1 As, with varying content x of aluminum, deposited on top of a semi-insulating GaAs wafer [see Fig. 1(a)].They were fabricated by means of molecular beam epitaxy.The thickness of the bottom layer is 5.0 μm with aluminum content x ¼ 24%, followed by a 1.5-μm-thick guiding layer with x ¼ 18% and secluded by a 1.0-μm-thick cladding layer with x ¼ 24%.Since the refractive index of AlGaAs decreases with increasing aluminum concentration [33], the light is primarily confined and guided inside the middle layer.Hence, the stacking constitutes a slab waveguide with single mode operation at the working wavelength λ 0 ¼ 1550 nm.Leakage into the substrate is negligible.As the photon energy is well below half the band gap, the waveguide is transparent at this wavelength and two-photon absorption plays no role [34].The nonlinear Kerr index n 2 of AlGaAs is ∼1.5 × 10 −17 m 2 W −1 , corresponding to a nonlinear susceptibility of third order, χ ð3Þ ¼ 5.9 × 10 −19 m 2 V −2 [35,36].
Microscope images taken from above the sample show the existence of elliptically shaped ripples, which were induced by thermal stress during layer growth [visible in Fig. 1 1 10 direction of the GaAs bed due to the orientation dependence of the interfacial energy of Ga on AlGaAs at deposition [37].In what follows, this direction also corresponds to the main direction of light propagation chosen in our experiments.The defects do not deteriorate the guiding properties considerably (scattering out of the waveguide provokes an attenuation length L loss ¼ 5.40 mm; see Sec.IV), but lead to slight anisotropic inhomogeneous variations of the effective refractive index n eff of the waveguide, thus causing weak but noticeable scattering inside the guiding layer.By performing a Fourier analysis of the microscope images, the distance between two adjacent scattering sites, that is, the scattering mean free path length, is determined to be l scat;∥ ≈ 20 μm and l scat;⊥ ≈ 3 μm parallel and orthogonal to the main propagation direction of light, respectively.
A TE-polarized Gaussian laser beam with elliptic intensity profile [full width half maximum of the intensity profile is 53 μm in the x and 1.7 μm in the y direction, see bottom left inset in Fig. 1(c)] is focused on the front face to excite the waveguide mode.The effective refractive index of the guided mode is n eff ¼ 3.33, so that λ eff ¼ λ 0 =n eff ¼ 465 nm and the transverse diffraction length z R ¼ πw 2 0 =λ eff ¼ 31 mm, where w 0 ¼ FWHM= ffiffiffiffiffiffiffiffiffiffi ffi 2 ln 2 p is half the width of the Gaussian field envelope.
We utilize two different types of lasers in the experiments: an infrared diode laser source at λ 0 ¼ 1550 nm producing continuous wave (cw) radiation at power P av ≈ 1 mW for linear reference measurements, and a pulsed laser source with high peak powers of up to P peak ≈ 1 MW at a constant pulse length of 500 fs to trigger nonlinear beam evolution.The intensity profile at the rear facet along the guiding layer of the waveguide is imaged onto a calibrated infrared CCD camera chip [see bottom righthand inset of Fig. 1(c)].The retrieved images form the basis of our experimental measurements.After taking account of a constant background and sporadic bad pixels, the intensity profile IðxÞ is normalized with respect to its mean intensity I 0 .Finally, a histogram of the resulting intensity ratios I=I 0 is recorded with bin size 0.5.
The front facet of the waveguide can be monitored by the insertion of a nonpolarizing beam splitter cube, which directs the fraction of the impinging light that is reflected at the front facet as well as the light that is scattered back from within the waveguide to another CCD camera [see top left-hand inset of Fig. 1(c)].Both components are distinguishable in the images in terms of their lateral extent and brightness.
Different propagation lengths are realized via samples of different sizes, all originating from the same AlGaAs wafer, while a number of ten different realizations within each sample with statistically equivalent features is accomplished by laterally shifting the samples between consecutive measurements.Changing the incoupling position, the guided light interacts with a different set of imperfections.
The histograms plotted throughout this paper represent the cumulative frequencies of the normalized intensities taken from the various measurements.The total number of evaluated intensity values per sample length and input power is about 5 × 10 3 for the short and intermediate samples and 8 × 10 3 for the sample of length L ¼ 30 mm.Reference measurements are taken using an unperturbed AlGaAs sample of length L ¼ 4 mm without index inhomogeneities, but with an identical progression of layers with varying aluminum concentration.

III. RESULTS IN THE LINEAR REGIME
The results of the measurements in the linear regime (achieved with the cw diode laser) are depicted in Fig. 2(a).The histogram plot of the ensemble averaged normalized intensity distribution at the rear facet shows how the typical statistics of an unperturbed Gaussian beam intensity profile (green solid line) is quickly transformed during propagation in the perturbed sample (turquoise to magenta solid lines).Multiple scattering randomizes the intensity profile into a speckle distribution for intermediate propagation lengths L ¼ 5 and 10 mm (turquoise and gray solid lines).However, as the propagation length increases, the tail of the distribution stretches out further, eventually surpassing the specklelike behavior stated by Eq. ( 1) (purple and magenta solid lines; for comparison, see black dashed line).It still resembles an exponentially decaying curve for high intensities, but with a decay factor significantly smaller than the one expected from a random process.Thus, extreme events are much more likely to occur than for a speckle distribution.Since multiple scattering at the random defects is supposed to meet the prerequisites of the central limit theorem and so lead to a randomized distribution, this comes as a complete surprise.It is the aim of this article to elaborate on this unexpected outcome and to provide a comprehensive explanation.
Comparable deviations from Rayleigh statistics have been found recently during linear beam propagation in slowly varying refractive index distributions with weak variations on the one hand [15,16] and randomly distributed Luneburg lens (LL) type defects with very strong index variations on the other [38].In both regimes extreme events could be observed as a result of the formation of caustics due to the (de)focusing effect induced by the spatial index variations.These caustics typically form after a few correlation lengths or LLs, but as was observed for the slowly varying index landscapes, they disappear after many such lengths when the randomization of the light field increases, thus eventually manifesting a speckle distribution.A property of these systems is that the evolution can be traced well with Hamiltonian ray optics, which is a clear indication that it is mainly governed by classical refraction.This is not the case in our setup, where the light undergoes multiple scattering (about 1500 times) before it reaches the end facet, where the intensity profile is finally observed with a CCD camera.The randomization along the propagation is so strong that the ray optics approximation ceases to give an adequate explanation of our experimental situation.
However, finite-difference time-domain (FDTD) simulations for strong index variations (around 30%) as stated in Ref. [38] already hint at the possible importance of backreflection.As we show, our system can be well described as an ensemble of weakly scattering pointlike defects that are separated by l scat , where the ratio between backward and lateral scattering plays a key role in explaining the occurrence of extreme events.
Naturally, it must be checked whether the measurements and the evaluation are sound.Since in the experiment we are only evaluating the intensity statistics on a bounded spatial interval along the rear facets, special care must be taken to ensure that the choice of the considered region does not deteriorate the validity of the statistics.A perfect generic specklelike distribution, for example, can only be expected from unlocalized field distributions.However, the initial Gaussian bell shape of the incoupled beam affects the envelope of the intensity profile even after many scattering events.The actual evaluation of the intensity profile observed at the back side of the waveguide structures is executed in more technical detail in Appendix A. There it is proven that considering only a finite interval of the intensity distribution tends to underestimate the frequency of extreme events and thus cannot explain the observed abundance of high-amplitude events discussed in this article.
In respective systems the scale on which a wave loses all information about its original propagation direction-the so-called transport mean free path length l trans -is equal to the scattering length l scat and of the order of the effective wavelength.Here the extensive amount of lateral scattering during propagation leads to different light paths crossing each other and to a buildup of transverse spatial correlation, which undermines the prerequisites of the central limit theorem and manifests itself in a growing number of highamplitude events.However, this is clearly not the case in our experiments, where scattering is much weaker, resulting in a transport mean free path length of l trans ¼ 1.87 mm as is deduced from transmission measurements in Sec.IV [46][47][48].This is much larger than both l scat and the effective wavelength λ eff inside the waveguide.
Although individual scattering events are weak, multiple scattering still dominates the field evolution in our samples, resulting in a classical diffusion process for light [49].In principle, this should be adequately replicated by numerical simulations utilizing the beam propagation method (BPM) [50], a well-established technique based on the paraxial approximation of scalar and the assumption of solely forward light propagation, yielding a Schrödinger-like differential equation for the components of the electromagnetic field.To our surprise the remarkable deviance from a specklelike intensity distribution observed in the experiment cannot be reproduced by these calculations for a large set of randomly distributed refractive index landscapes (see Appendix B).
In this article we suggest that the small, but noticeable backscattered field components not accounted for in the  1)].For longer propagation lengths (L ¼ 15, 20, 30 mm), extreme events are more probable than anticipated from a random distribution.(b) Histogram of normalized intensity distribution at different propagation lengths z from simulations based on the stray matrix method with plane wave excitation and weak, but prevailing backscattering.The exact scattering amplitudes are stated in the pictogram displayed in the bottom left-hand corner (clockwise, starting from the top): backward jrj, leftward jr ⊥ j, forward jtj, and rightward jr ⊥ j (see also Fig. 5).For increasing z the probability distribution quickly surpasses a speckle distribution, before asymptotically approaching a heavy-tailed limiting function between z ¼ 20 and 30 mm.See also the inset for the dependence of the probability density at a fixed intensity value I ¼ 7.0I 0 on z.(c) Numerical results with less dominant backscattering.For intermediate propagation lengths the frequency of extreme events is greatly diminished.However, partial diffuse reflection at the rear facet increases the probability of extreme events in its vicinity (also visible in the inset).
BPM but observed in our experiments [see top left-hand inset in Fig. 1(c)] cause the detected abundance of extreme events.Our conclusion is supported by numerical results obtained from a discrete simulation scheme.The computed results are depicted in Fig. 2(b).There it can be seen thatalthough being very weak-backscattering accomplishes the expected deviation from a speckle distribution.Provided that it prevails over lateral scattering, it generates an abundance of extreme events in very good accordance with the experimental data.If backscattering is less dominant [see Fig. 2(c)], extreme events are greatly diminished and the tails of the intensity distribution approach the exponential trend line of a speckle pattern.In Sec.V, we give an elaborate introduction to the underlying numerical scheme, which was brought up first by Edrei et al. [31] to investigate the strong scattering case.In this article it is expanded to the weakly scattering regime.
But before we delve into the numerics, first we revisit diffusion theory and properly apply it to the twodimensional waveguide structures at hand.After what follows, we will be able to extract the characteristic features of the waveguide structures by means of experimentally accessible parameters in order to endorse our reasoning and to insert them into the numerical model.

IV. 2D DIFFUSION EQUATION
The aim of this section is to introduce all the critical quantities describing the scattering process analytically and to motivate the way that they are measured experimentally.To this end we first derive the two-dimensional diffusion equation from the respective radiative transfer equation, largely reproducing the derivation of the three-dimensional analogue as stated, for example, in detail in Refs.[39,49], but emphasizing the subtle changes due to the reduced dimensionality.Later we provide solutions of this equation to be compared with experimental results.

A. Derivation of the 2D diffusion equation
The time-dependent, emission-free (no source term) radiative transfer equation for the diffuse radiance Iðr; n; tÞ, corresponding to the energy flux at position r into direction n, is given by [39,51] where t scat and l scat denote the mean time and the mean free path length between two scattering events, respectively, ∇ ¼ ð∂=∂x; ∂=∂yÞ is the two-dimensional gradient operator, and pðn; n 0 Þ is the so-called phase function of one scatterer, describing the scattered intensity into direction n 0 from an incident plane wave traveling in the direction n.
Accordingly, n ¼ ðsin ϑ; cos ϑÞ and n 0 ¼ ðsin ϑ 0 ; cos ϑ 0 Þ stand for arbitrary unit vectors, the angle between which shall be denoted by θ.The integral on the right-hand side of Eq. ( 2) is performed over the whole two-dimensional angular range, so that R Ω dn 0 ≡ R 2π 0 dϑ 0 .Note that here and in the following the normalization factor on the right-hand side differs from the 3D case treated in Ref. [39].
Next always holds for rotationally symmetric scatterers.Hence, we are allowed to multiply the first term on the right-hand side of Eq. ( 4) by n 0 • n 0 ¼ 1 and associate one of the dashed unit vectors with the above integral to find the expression to be the average cosine of the scattering angle.By utilizing this result, we can rewrite Eq. ( 4) as In the following we assume z to be the main propagation direction.The respective component of the current density is given by In the so-called diffusion approximation it is assumed that the diffuse intensity is almost independent of the direction n, except for a small excess in the forward, i.e., positive z direction: where γ is an unknown constant.By inserting this ansatz into the definition of J z stated in Eq. ( 7), we find γ ¼ 1=π.
Finally, after some calculation, Eq. ( 6) takes the form ∇Iðr; tÞ: In the static case there is no time dependence and we can neglect the first term on the left-hand side, so that we get where we introduced the diffusion constant D as Note again the different factor in the denominator compared to the 3D case [39].Here we have also defined the velocity of diffuse light transport v ¼ l scat =t scat and the transport mean free path length: The latter can be used to estimate the (average) phase function of one scatterer (see Sec. IV C).In contrast to the scattering mean free path length l scat , which is the distance between consecutive, potentially weak scattering events (mainly emitting in the forward direction), the transport mean free path length l trans can be thought of as the distance, after which all information about the propagation direction of the incident light beam is completely lost.This might happen only after many consecutive weak scattering events.In the case of strong scattering hcos θi ¼ 0, the transport mean free path length equals the scattering mean free path length.In fact, in scattering theory one often assumes strong scattering in order to significantly simplify the math involved, only to invoke weak scattering in the end result in virtue of Eq. ( 10).
When inserted in Eq. ( 3) and again assuming static behavior, Eq. ( 9) leads to the diffusion equation for the intensity IðrÞ, where we introduce the absorption length L loss so that Keep in mind that in general the diffusion constant D might depend on the position.However, here we treat it as a constant, which can be omitted from Eq. ( 11).

B. Solution of the diffusion equation
Our desired task is to extract the absorption length L loss and transport mean free path length l trans from experimental transmission measurements.To do so, we consider the propagation of light through a slab waveguide of length L with fixed mode profile in the y direction (see Fig. 3).The input and output facets of the waveguide are placed at z ¼ 0 and z ¼ L. The effective refractive index of the waveguide mode is n eff and that of the surrounding environment is n env .Spatially coherent radiation entering the slab from the left is transformed into diffuse light during its propagation through a layer roughly as thick as the transport mean free path length l trans .For the sake of simplicity we neglect the details of this transformation and consider a source of diffuse intensity being placed at z ≈ l trans .
In the static case and in the presence of absorption, the diffuse light intensity thus obeys the diffusion equation [see Eq. ( 11)], with Δ ¼ ∇ 2 and where we assume a delta-peak-shaped source of strength S 0 at z ≈ l trans .Since in our simplified model all diffuse light is originating from the source layer, there is no incoming flux through the boundaries at z ¼ 0 and z ¼ L. Hence, the energy flux is determined by energy loss through the boundary and is expected to be proportional to the diffuse intensity I leading to [46] Iðx; zÞ and Here the so-called extrapolation length z 0 is introduced, which is the hypothetical distance after which the diffuse intensity vanishes due to the light being scattered back to the positive (or negative) z direction completely.The form of the boundary conditions and the value of z 0 , which generally depends on the reflectivity of the interfaces (and therefore the ratio of n eff and n env ), is derived in detail in Appendix C. The resulting profile of the diffuse intensity corresponding to the Green's function of the system, that is, the solution of Eq. ( 12), is given by [39,47] where To calculate the diffuse intensity I out leaving the slab waveguide through the rear interface at z ¼ L, similar to the above consideration, one often assumes that no scattering or considerable absorption takes place in a layer next to the rear interface roughly as thick as l trans and that, hence, the emitted intensity corresponds to the diffusive intensity a distance l trans away from the output facet; formally, For a slab length L ≥ 2z 0 , this results in To gauge the source power S 0 and retrieve an expression for the total transmission TðLÞ of the waveguide, we assume a short section of the waveguide with length L ¼ 2l trans to have a transmission of Such an assumption is in agreement with our model, because the source at z ¼ l trans emits the same amount of power to the front and rear facet resulting in reflection and transmission being equal.As long as the transport length l trans is reasonably smaller than the absorption length L loss , losses are negligible and half of the power entering such a short waveguide section must be transmitted.With this, we finally get for the total transmission through a multiply scattering, fairly lossy two-dimensional waveguide, where I out ðLÞ is given by Eq. ( 16).

C. Characterizing the sample
To determine the scattering strength of the defects in the waveguide structures used in the experiments, measurements of the mean total transmission through the waveguides as a function of the various sample lengths L were performed.When doing so one must consider that only a fraction of the actual impinging power is transmitted through the front facet of the sample due to the Fresnel coefficient for normal incidence and is then converted to diffuse intensity.If this is taken into account, then by comparing the experimental results with Eq. ( 17), one can readily extract both the attenuation length L loss ¼ ð5.41 AE 0.48Þ mm and the transport mean free path length l trans ¼ ð1.87 AE 0.13Þ mm by means of a two-parameter fit [see gray squares and solid red line in Fig. 3(c)].Comparing the quality of the fit to another model based on Eq. ( 15), where possible losses along the order of l trans are specifically taken into account [solid green line in Fig. 3(c)], the assumption l trans < L loss made above seems well justified.
Comparing l trans with the distance between scattering events in the longitudinal direction l scat;∥ ≈ 20 μm, we gain knowledge of the scattering strength by means of Eq. ( 10) [32,39], finally finding hcos θi ¼ 0.989.In general, hcos θi ¼ 0 implies strong (Rayleigh) scattering and hcos θi ≲ 1 represents weak (Mie-like) scattering.So, indeed, scattering in the perturbed AlGaAs samples is very weak, while the scattering length l scat is so short that multiple scattering is ensured, rendering the deviation of the intensity statistics from a random light distribution the more peculiar.
Next we are striving to estimate the average directiondependent scattering amplitude of the defects in the waveguide.However, as we present a discrete numerical model in the next section, we do not need to derive the full angular dependence of the phase function pðcos θÞ, and it will be sufficient to consider the scattering amplitudes into the four major directions only.
In order to roughly approximate the scattering amplitudes and-more importantly-the ratios between them, we assume that the distribution of bright and dark streaks visible in the backlit microscope image depicted in Fig. 1(b) qualitatively mirrors the spatial variation of the refractive index distribution in the perturbed waveguide structure.After applying a suitable high-pass filter and defining a brightness threshold, we are left with a monochromatic distribution [see inset of Fig. 4(a)] that can readily be used as a template for the spatial index variation in a FDTD simulation (commercial software Lumerical FDTD Solutions).The extracted section is 20 μm long and 4 μm wide.To avoid specular reflection at the section boundaries, straight edges and corners are rounded off and the resulting template is convoluted with a Gaussian with width w ¼ 0.5 μm ≈ λ eff .We utilize a plane wave source polarized as in the experiment, i.e., parallel to the plane of incidence (TE polarized in terms of the waveguide structure) and monitor the angle-dependent scattered field amplitude in the far field (use of a total-field scatteredfield source and subsequent far field transformation).An exemplary numerical result is illustrated in Fig. 4(b).The logarithmic intensity plot in the outer frame reveals the phase function.In Fig. 4(c), the angle-dependent scattered intensity distribution is plotted for an exemplary index variation Δn ¼ 0.001 between the background and the  defect.Clearly, most of the radiation is scattered in the forward direction (scattering angle θ ¼ 0°) and only a small part is scattered to the sides, with backscattering (θ ¼ 180°) notably dominating over transverse scattering.
As we are only interested in discrete scattering directions, namely the forward, backward, left, and right direction, we shall integrate over some angular blocks of 90°with the respective major directions in their center.Still, the FDTD simulations suggest that even for small index variations significantly more radiation is scattered in the back direction than to the sides (up to 10 times for very small index variations).
In the following, we conservatively estimate jrj ≈ 2jr ⊥ j, which as will turn out is enough to numerically provoke the abundance of extreme events as observed in the experiments.An even stronger imbalance would result in even more pronounced deviations from a speckle distribution.Now, for discrete scattering the computation of the average cosine of the scattering angle defined in Eq. ( 5) simplifies to To estimate the moduli of the scattering amplitudes in the four major directions, we state which is a direct consequence of the conservation of energy, with scattering out of the waveguide being the only relevant loss mechanism.By combining the determined result for the average cosine of the scattering angle, the estimate for the ratio between the backscattering and the lateral scattering amplitude, and Eqs. ( 18) and ( 19), we eventually find the scattering amplitude in the forward direction to be jtj ¼ 0.997, the scattering amplitude in the backward direction jrj ¼ 0.063, and for transverse scattering jr ⊥ j ¼ 0.032.For comparison, the case of strong scattering would correspond to jtj ¼ jrj ¼ jr ⊥ j ¼ 0.5, stating that all incoming light would be scattered isotropically.

V. STRAY MATRIX APPROACH
After having properly described the samples from the experiments by their characteristic features, next we strive to replicate the experimental situation with a numerical model.To study the impact of backscattering we make use of the stray matrix approach proposed by Edrei et al. [31] (see also [32]) and represent our randomly perturbed waveguide samples in a very simplified manner by a discrete two-dimensional square lattice of pointlike scatterers (see Fig. 5).This approach has a tremendous advantage over rigorous finite element calculations regarding computation time.For a domain size as large as considered in this analysis there is no alternative, which still incorporates multiple scattering that is inherent to our experimental situation.Also note the necessity of an appropriate ensemble size for valid statistical conclusions.
In the following we close some holes in the justification of the applicability of the simulation scheme not touched upon in the original references and transfer the scheme to the limiting case of weak scattering.
Each scatterer is connected only to its adjacent neighbors on the left, right, front, and back, so that the scattering sites constitute a discretized network with four bonds leaving each scatterer.The effect of one single scattering event and, hence, the field evolution after one time step, is described by a linear set of equations a ¼ Ŝb, where b is an input vector at a specific scattering site (consisting of the field components injected from the back, left, right, and front) and a is the respective output vector (consisting of the field components leaving the scattering site to the back, left, right, and front).Where no ambiguities arise, the notion of the position dependence [indicated by the indices n, m as in Fig. 5(c)] will be dropped.The 4 × 4 matrix Ŝ has the components S pq ¼ m pq e iðφ p þφ q Þ ; with p, q ¼ 1, 2, 3, 4 [see Fig. 5(b)].The factors φ p and φ q represent the phase accumulated along the bonds p and q, respectively (note that each scattering matrix Ŝ contains only half the total bond length between adjacent scattering centers).The matrix m is unitary and symmetric (ensuring conservation of energy and time-reversal symmetry) and contains the complex valued scattering amplitudes, 5 ; where r ¼ jrje iφ r is the reflection coefficient, t ¼ jtje iφ t is the transmission coefficient, and r ⊥ ¼ jr ⊥ je iφ ⊥ is the transverse scattering amplitude.Here we have imposed left-right symmetry, which should be present at least in average.Note that in this form the scatterers are assumed to be rotationally invariant so that, for example, t is the scattering amplitude of the nondeflected field component no matter from which direction the initial field was approaching.
Since only relative phases are of importance, for the sake of simplicity we set φ ⊥ ¼ 0.Then, due to m being unitary, we end up with the following constraints: As we show in the next section, if all scatterers which comprise the network, are represented by identical matrices and have uniform distances between them, this system reproduces wave diffraction in an unperturbed planar waveguide well (see Fig. 6 and the discussion of the dispersion relation in the following section).On the other hand, in order to introduce randomness into our model and simulate a chaotic distribution of scatterers, we assume the accumulated phase to fluctuate around a mean value with a normal distribution derived from our experiments (see Fig. 7 for the effect of phase variations).
In the simulations the size of the computational domain is N L ¼ 1500 scatterers in length and N W ¼ 501 scatterers in width.The edges of the computational domain have to be treated separately and require the implementation of boundary conditions.All numerical results presented in this article are achieved with plane wave excitation and periodic boundary conditions on the sides.The front and rear facets are partially reflective according to the reflection coefficient for diffuse intensity of the interface between the waveguiding medium and the environment.The reflection coefficient for diffuse intensity differs from the Fresnel coefficient for specular reflection (see Appendix C).For the evaluation of the intensity profiles we have always utilized the forward propagating field component on each site.Numerical results as pictured throughout the article are built from the average outcome of 256 different bond phase realizations, where each simulation was allowed to run until a steady state was reached (about 2 hours per realization on a desktop computer).
Each iteration of the simulation consists of multiplying the input vector b at each scattering site with its respective scattering matrix Ŝ and then copying the components of the output vector to the corresponding input vector slots of all four adjacent scattering sites.This procedure is repeated until variations of the intensity profile between single iterations become negligibly small (this is the steady state).In case of nonlinear propagation, as discussed in Sec.VI, no steady state is reached in this strict sense.Still, after a certain amount of iterations the average overall spatial intensity distribution does not change anymore and only fluctuates locally.This is considered sufficient to analyze its probability density function of the intensities.

A. Spatial dispersion relation
To test the applicability of our model and to determine appropriate input parameters, we study the case of the unperturbed slab waveguide first.Monochromatic wave propagation in uniform two-dimensional media is completely determined by the so-called spatial dispersion (or diffraction) relation, which imposes a dependency between the longitudinal component k ∥ and the transverse component k ⊥ of the wave vector.
While in the isotropic case ⊥ (which is the definition of a circle) always holds, this is not necessarily true in periodic media [52].Given some k ⊥ , the effective propagation direction is perpendicular to the graph of the spatial dispersion relation at that point.Its second derivative hints at the divergence, i.e., defocusing behavior of plane waves with slightly varying values of k ⊥ .In the following we show that-depending on the specific choice of the bond phase parameters φ-our discrete network model also displays a circular diffraction relation and is thus indeed capable of simulating homogeneous space.
In order to do so, we first stipulate a regular (translational and rotational invariant) lattice with φ 1 ¼ φ 2 ¼ φ 3 ¼ φ 4 ≡ φ bond for all scatterers.We may put φ bond ¼ kl scat =2, with the wave number k ¼ 2πn=λ 0 (with a certain refractive index n and the wavelength λ 0 ) and a physical distance l scat between scatterers.
The time-step-dependent global field distribution along the mesh of scatterers is characterized by the output and input vectors a t nm and b t nm at sites n, m, respectively.Since output vectors at a time t generate input vectors of adjacent scatterers at time t þ δt [see Fig. 5(c)], we can formulate the global field distribution after one time step δt generally as a function of the field distribution before scattering as Note that in the perturbed case the scattering matrix Ŝ ¼ Ŝnm would be site dependent.
While boundary conditions must be implemented in the limited spatial domain of the simulation, in this analysis we assume an infinitely extended domain in order to derive the dispersion relation.Assuming no time evolution and inserting the plane wave ansatz b nm ¼ b 0 e iαn e iβm (with α ¼ k ⊥ l scat and β ¼ k ∥ l scat ), we arrive at This homogeneous set of equations has only nontrivial solutions for detð ĜÞ ¼ 0. In general, detð ĜÞ is complex, but we may consider the real or the imaginary part separately, since both give the same results.We can then calculate the resulting dispersion relation β ¼ βðαÞ.

Strong scattering
For the case of strong scattering (jrj ¼ jtj ¼ jr ⊥ j ¼ 0.5), the structure of Ĝ simplifies drastically and we are left with The shape of the dispersion relation depends on φ bond .If cos 2 φ bond ≈ 1, this equation has real valued solutions only for α ≪ 1 leading to β ≪ 1, so that we are able to approximate leading to which is the equation of a circle with radius 2 ffiffi ffi 2 p φ bond .Hence, for cos φ bond ≈ 1, the dispersion relation of a homogeneous system with the effective wave number k is reproduced.Remember that for strong scattering, l scat ¼ l trans .

Weak scattering
In the weakly scattering case the expression for the dispersion relation is in general quite complicated, but can still be solved analytically.As in the case of strong scattering, a critical phase φ bond;crit can be found at which the dispersion relation dwindles to a point, but for all jrj ≥ jr ⊥ j is asymptotically circular as φ bond approaches φ bond;crit , thus ensuring physically correct diffraction.For jrj < jr ⊥ j, the dispersion relation is distorted into a diamondlike shape, leading to physically incorrect diffraction.But as has been confirmed by means of FDTD simulations in Sec.IV, assuming jrj > jr ⊥ j is very appropriate in the experimental context and thus poses no restriction on the applicability of the discrete numerical scheme.In Fig. 6(a), we have plotted the critical phase for different values of the modulus of the scattering amplitude in the forward direction jtj under the assumption jrj ¼ 2jr ⊥ j.
Just like multiple weak scattering has the same effect as strong scattering after the light has passed N ¼ l trans =l scat ¼ 1=ð1 − hcos θiÞ ¼ 1=ð1 − jtj 2 þ jrj 2 Þ scattering sites, imposing weak scattering in the unperturbed simulation can be thought of as effectively increasing the resolution of the simulation by a factor of N in comparison with the strongly scattering case.This is reflected in the approximately identical values of the effective wave numbers k eff taken from the radius of the asymptotically circular dispersion relation in units π=l trans displayed in Fig. 6(b).The comparison of two numerical results implementing isotropic and weak scattering can be found in Figs.6(c) and 6(d).
At first sight, any criterion for φ bond leading to a circular dispersion relation implements a notable restriction.However, we show now that this restriction is quite weak, when the numerical model is compared with the experimental situation.

B. Distribution of the mean free path lengths
The elliptically shaped defects in the AlGaAs waveguide structures used in our experiments impose an anisotropic distribution of scattering mean free path lengths.This can be accommodated in the simulation by considering two different phase factors, φ 1 ¼ φ 4 ¼ φ long ≈ 43π and φ 2 ¼ φ 3 ¼ φ short ≈ 6.5π, corresponding to l scat;∥ and l scat;⊥ with respect to the effective wavelength λ eff (see Sec. II for a description of the setup).As a consequence, the dispersion relation loses its fourfold rotational symmetry, but remains axially symmetric.
However, it should be noted that the dispersion relation solely consists of trigonometric functions and, hence, depends only on mod ðφ p ; 2πÞ (with p ¼ 1, 2, 3, 4).Consequently, in this model any anisotropy arises solely due to differing phases between the longitudinal and the lateral direction in the order of 2π.Since the distance between adjacent scattering sites in the utilized waveguides corresponds to a phase accumulation of several times 2π and the uncertainty in their measurements exceeds π=2, it is in fact reasonable to roughly neglect any deviance in mod ðφ p ; 2πÞ.Hence, within this uncertainty the dispersion relation can be regarded as isotropic.
Practically, depending on the chosen values of the scattering amplitudes jrj, jtj and jr ⊥ j, the bond phase φ bond is gauged so as to result in the spatial dispersion relation of the unperturbed discrete system to be circular.This ensures the system to simulate homogeneous space as apparent from the normal diffraction of a Gaussian beam as depicted in Fig. 7(a).
Still, the two different length scales enter the simulation through the variation of the bond lengths in the perturbed case, which are considered to be distributed normally with a standard deviation σ taken to be 10% of the respective mean values φlong=short , namely, σðφ bond;long Þ ¼ 4.3π and σðφ bond;short Þ ¼ 0.7π.Although it is our aim here to alter the mean free path length between scattering sites thus randomizing the discrete network, it should be noted that this  approach is formally identical to a two-dimensional refractive index variation inside the waveguide structure.It has been found that a certain amount of randomness is crucial and that too small variations will not lead to the deviation of the intensity distribution from a speckle pattern, as found in our experiments.A series of simulations for a varying standard deviation σ of the Gaussian phase distribution is depicted in Fig. 7. Depending in detail on the actual scattering amplitudes, a standard deviation σ of about 0.1% gives rise to the formation of large speckles only.This case is comparable to beam propagation in a slowly varying, but otherwise random refractive index landscape with a macroscopic correlation length.In contrast, a σ amounting to 10% leaves adjacent scatterers mostly uncorrelated, but allows for sporadic resonant substructures, in which the intensity can become excessively high.Further simulations suggest that a σ of about 1% of the mean value can be considered random enough to provoke the emergence of outstanding highamplitude events.
We end this section by noting that a cross-check has been performed on the numerical results to extract the transport mean free path length and the attenuation length in return from the simulations that emerge from our now established model.The same two-parameter fit based on Eqs. ( 16) and ( 17) that was applied earlier to the experimental results yields l trans ¼ ð1.84 AE 0.12Þ mm and L loss ¼ ð5.32 AE 0.44Þ mm.Since these values are in very good accordance with the experiments, we are confident that the numerical model mimics the key features of the setup adequately.

C. Linear scattering
Certainly, there are many ways to implement the stray matrix scheme numerically.To perceive the field distribution after a single time step, one might iterate over all the scattering sites performing a nm ¼ Ŝnm b nm with 1 < n < N W and 1 < m < N L individually, then copy each component of the output vector a nm to the appropriate counterpart of the input vectors b nAE1mAE1 of the next time step of the adjacent sites and progress to the next scatterer.However, this procedure turns out to be relatively time consuming.
In contrast, it is particularly enlightening to achieve all the necessary computation for one time step with a single matrix multiplication.To do so, all local input vectors b nm are assembled systematically in one large input column vector b total of dimension 4 • N L • N W .In the following, it should be kept in mind that this column vector can be reordered into a two-dimensional grid at any time, hence illustrating the complete field distribution along the bonds over the whole considered domain.
In this section, we focus on linear propagation and neglect scattering losses.Both the individual scattering processes and the propagation along the bonds can be described by where M is a ð4 mainly consists of groups of constant values along the main diagonal.Hence, although it is very large, most of the matrix's entries are in fact zero so that it can be conveniently represented by a sparse matrix, which still fits into the memory of most computers even for large computational domains.Not only does this approach improve computation speed by an order of magnitude, the consideration of one sole evolution matrix is an interesting mathematical concept by itself as it opens up a way to understand our system in terms of linear algebra.We may strive now to understand the matrix's action on a given field distribution encoded in b total in terms of its eigenvectors.The eigenvectors of M represent stable field profiles across the computational domain which do not change after one iteration apart from being multiplied uniformly with a complex number, which is given by the corresponding eigenvalues λ.Hence, these particular field profiles are only attenuated and/or phase shifted between consecutive time steps.They correspond to configurations where numerous local cavities are excited and coupled in just the right way so as to form a macroscopic stationary state.The modulus of λ determines how fast its field decays for vanishing excitation.Those eigenvectors with eigenvalues of modulus close to 1 behave like cavity modes with a quality factor Q ¼ 2πn eff d=λ 0 ð1 − ajλj 2 Þ, where λ 0 =n eff is the effective wavelength in the waveguide, d ¼ l scat;∥ is the lattice period in the simulations, and a is the albedo.In our case the Q factors of the cavity modes reach values up to about 10 5 (see discussion below).The resonance frequency of such modes is encoded in the phase of λ.For example, a cavity mode with very high Q factor resonating at the excitation frequency corresponds to a real eigenvalue slightly smaller than þ1.
There are in fact 4 • N L • N W distinctive eigenvectors.For plane wave excitation at the front facet, all of them are excited approximately uniformly.In what follows we investigate the statistical properties of the eigenvectors for varying backscattering strength of the local scatterers.If backscattering and lateral scattering are equal but weak (jrj ¼ jr ⊥ j ≪ jtj), the field statistics is Gaussian, yielding a specklelike intensity distribution [see blue lines in Figs.8(a) and 8(b)].The intensity statistics of the eigenvectors in this case corresponds very well to the iterative simulation results of the total intensity distribution as displayed in Fig. 2(c).
As can be seen in Fig. 8(c), the intensity profile of the respective cavity modes is not just chaotic, but shows an increased correlation along the propagation direction resulting in pronounced scars caused by the dominant forward scattering.
If backscattering prevails lateral scattering several times, as it is the case in our samples, the field statistics gradually transforms from a Gaussian to a sech shape [see red and yellow solid lines in Fig. 8(a)].Thus, heavy tails form and extreme intensity events are promoted [see Fig. 8(b)].Again, as visible in Figs.8(d) and 8(e), a certain coherence length can be witnessed in the intensity profile of the respective eigenvectors represented by highly pronounced scars.This extreme focusing in longitudinal cavities formed via backscattering is promoted by the relatively small lateral scattering, which hinders light to escape in the sideward direction.
If backscattering completely dominates lateral scattering (10 times or more), the system resolves into an ensemble of weakly coupled one-dimensional transmission lines.In this case the field components are Laplace distributed [see green line in Fig. 8(a)].Looking at the intensity profile at a fixed propagation length [see Fig. 8(f)], we find the characteristic log-normal intensity distribution, which is well known for strictly one-dimensional randomly disturbed channels [53][54][55].
This analysis shows that the abundance of extreme events in our system can be traced back to the level of the eigenvectors of the evolution matrix M and thus the excitation of quasi-one-dimensional cavity states.It illustrates that our experiment takes place in between isotropic two-dimensional and highly directional one-dimensional scattering.
But how is it then in the light of the central limit theorem that a superposition of many such field distributions as suggested by the uniform excitation of all eigenvectors can still yield extreme events as evident from the intensity histogram plots in Figs.2(a) and 2(b) and is not simply specklelike, as the central limit theorem would suggest?Besides a small amount of leakage to the surroundings during scattering, which is treated here as affecting all scatterers equally and thus is accounted for by a scalar factor a that is not of further interest, losses mainly occur at the front and end facet only.As a consequence, all eigenvalues of M have moduli very close to but slightly smaller than 1 [see an exemplary plot of eigenvalues for a computational domain with N L ¼ 100 and N W ¼ 31 in Fig. 9(a)].
Under consecutive excitation with a plane wave, in the steady state (b tþδt total ¼ b t total ) the amplitude of each eigenvector is proportional to 1=ð1 − λÞ.Hence, from all excited eigenvectors only those survive (or are severely enhanced) for which the real part of λ is very close to þ1.As it turns out, only a handful of eigenvectors actually contribute to the field distribution then [see Fig. 9(b), where despite of 12 400 eigenvectors in total, the steady state is constituted of only about a dozen eigenvectors with even less standing out significantly].Under this condition the central limit theorem does not hold and the intensity distribution of individual eigenvectors dominates the intensity distribution of the steady state.
If, however, losses due to scattering or absorption become too big and the moduli of the eigenvalues drop significantly below 1, more and more eigenvectors contribute to the steady state equally.Then indeed it follows from the central limit theorem that the superposition of quasirandom field components adds up to a total field whose distribution is Gaussian [see Fig. 9(c)].This is confirmed by the iterative approach, where small additional scattering losses lead to a "delay" of the formation of the heavy-tailed probability distribution (see Appendix D for numerical results).

VI. RESULTS IN THE NONLINEAR REGIME
While we already see extreme field enhancement in the linear case, according to common expectation this effect should be intensified by incorporating a focusing nonlinearity of the Kerr type, often thought of as playing a key role in generating extreme events [1][2][3]8,[18][19][20][21][22][23][24][25]27,30,38,56].Here we assume the simplest case of a Kerr NL and suppose the refractive index to be changed proportional to the intensity of the propagating light like n ¼ n 0 þ n 2 I.This is well justified for our waveguide having a fixed mode profile and polarization structure.In the numerical model this can be thought of as implementing an additional phase accumulation between scattering centers proportional to the local intensity along the bonds.Since the Kerr NL of AlGaAs has a positive value, one expects it to counteract diffraction, eventually leading to an intensitydependent self-focusing of the local field amplitude.
The impact of nonlinear beam evolution in a random medium has been at the core of only a few experimental or theoretical considerations.In these related systems, the interplay of disorder and nonlinearity was shown to affect weak [57][58][59] and strong [60,61] localization as well as pulse broadening [56].

A. Experiments
Overall the results of the experiments executed with the pulsed laser source for low pulse powers align very well and cross over continuously with the respective measurements conducted with the cw laser [reproduced as dark blue dashed curves in Figs.10(a)-10(c)].From this, we conclude that the shorter coherence length of the pulses plays no significant role.Increasing the laser power leaves the pulse length constant.Indeed, we see a power driven rise of extreme events for the shortest perturbed sample with L ¼ 5 mm [see Fig. 10(a)].Here, modulation instabilities are expected to grow, finally leading to solitons that arise from the weakly disturbed field distribution.This results in an increase of the frequency of high-amplitude events for growing input powers.But for longer samples with a nonspecklelike distribution in the linear regime, the trend is inverted: For increasing input power extreme events become less likely and the tail of the distribution restores a specklelike behavior in complete accordance with Eq. ( 1) [see Figs.10(b) and 10(c)].For the longest sample with L ¼ 30 mm, the abundance of extreme events that can be observed in the linear regime is completely suppressed.An overview of the power-dependent development of the probability density function for varying sample lengths at one fixed intensity I ¼ 7.0I 0 is plotted in Fig. 10(d).
Although it still possesses focusing properties, the nonlinearity merely tends to randomize the intensity pattern by destroying the coherent superposition of light scattered back and forth between the defects.Consequently, resonant field enhancement is prevented, a trend which becomes more and more pronounced for increasing sample length.This is illustrated most clearly in an exemplary plot of an experimentally determined intensity profile for increasing input powers and L ¼ 30 mm in Fig. 10(e).As the nonlinearity grows, the intensity distribution becomes more chaotic and the power is allocated over an increasing number of spikes, so that the most outstanding ones are diminished.

B. Simulations
Despite the simplicity of our discrete numerical model, it reproduces self-focusing in the homogeneous waveguide with uniform bond phases very well (see Fig. 11).The NL successfully compensates for the divergence of the injected Gaussian beam, up to the formation of a solitary solution.After that successful test of our model we are ready to investigate the nonlinear propagation in the perturbed waveguide.Here we consider individual scattering events as linear [57] and assume that the nonlinearity only causes an intensity-dependent phase accumulation during propagation between scatterers.To this end additional phase shifts proportional to the squared field modulus at the input (right before scattering) and output (right after scattering) ports are added to the linear phase distribution.The experimentally determined power-dependent intensity profile as well as the statistics match well with the results of the simulations depicted in Fig. 12.Since the linear phase is already randomized, one might expect the influence of additional phase terms to be negligible in average.But as can be seen by comparing Fig. 12(a profile at low and high NL].It coarsens the field shape significantly, leading to a kind of thermalization [see also Fig. 12(d)].Apparently, the formation of randomly spread resonating substructures, which create coherent amplification by chance in the linear regime, seems to get lost for increasing levels of nonlinearity.

VII. CONCLUSION
In conclusion, we observe an abundance of extreme wave events at low power in a slab waveguide with embedded random but only weakly scattering obstacles.Obvious deviations from a specklelike intensity distribution are caused by backscattering, which although small in amplitude dominates over scattering in the transverse direction and so occasionally gives rise to the emergence of resonant cavities.This mechanism is in stark contrast to the buildup of spatial correlation through the crossing of light paths in strongly scattering media, but its impact on the intensity probability distribution is quite comparable.In addition, against the common notion, here a focusing nonlinearity of the Kerr type suppresses those extreme events by preventing resonant field enhancement in the random cavities, thus reinforcing a specklelike intensity pattern.
The findings we present in this paper are applicable to all linear wave dynamical systems subject to weak multiple scattering in stationary random environments, under the prerequisites of reduced lateral scattering amplitudes (effectively forming a 1.5-dimensional system) and only minor losses.This includes optics and microwave dynamics, acoustics, seismic waves, electron transport, Bose-Einstein condensates, and interplanetary scintillation.Understanding the generation of nonspeckle intensity distributions is thought to be of high value in examining biological tissue and improving the yield of solar cells and broad area laser diodes.Anisotropic phase functions as considered in this work appear naturally for Mie-like scattering or can be realized intentionally with photonic antennas.Our results give evidence that the role of nonlinearity concerning the origin and formation of extreme wave events is still poorly understood and even counterintuitive in the regime of weak multiple scattering.Since until now weakly scattering systems were considered to only exhibit exponentially damped highamplitude events, our findings underline the need to rethink the risk potential of hazardous extreme waves in such environments.

ACKNOWLEDGMENTS
We thank Roberto Morandotti for providing the unperturbed AlGaAs reference sample.We gratefully acknowledge support by Deutsche Forschungsgemeinschaft (DFG) by providing a pulsed laser (GZ: INST 90/822-1 FUGG).This work was also supported by the Cluster of Excellence Engineering of Advanced Materials (EAM).

APPENDIX A: INFLUENCE OF A LOCALIZED SOURCE ON THE INTENSITY STATISTICS
In the experimental setup an elliptically shaped Gaussian laser beam is coupled into a highly perturbed waveguide structure and acts as a coherent light source.As the light propagates further and further through the medium, coherent radiation is gradually converted into diffuse intensity, forming an extended specklelike pattern with an ever fading localized Gaussian background.This can in fact be thought of as the sum of different field profiles: the Gaussian and the specklelike field distribution.Here we show that a slight Gaussian background tends to suppress extreme events, if the profile is normalized with respect to the mean intensity within a reasonable interval, and that the number of extreme events found in our experiments might thus be slightly underestimated.
As schematically described in Fig. 1(c), in the experiments the intensity profile is monitored at the rear facet of the waveguide samples.By laterally shifting the individual samples between single measurements, ten different realizations of the refractive index landscape with equal statistical features and sample length are obtained.To standardize the evaluation, the various measured intensity profiles IðxÞ from a certain sample length are centered with respect to their first moment and then superimposed to determine the average width of the extended speckle pattern.The resulting mean intensity profile ĪðxÞ is fitted with a bell-shaped Gaussian to extract (half) its width w 0 .For the subsequent evaluation of the individual intensity profiles, a region of interest with the transverse coordinate x ∈ ½−2.5w 0 ; 2.5w 0 symmetrically located around the first moment is chosen.As we demonstrate, this region of interest renders the evaluation sound in terms of the statistics of extreme events and ensures the general validity of our statements.
First, imagine a Gaussian bell-shaped field distribution with amplitude A and width ffiffi ffi 2 p w 0 (so that the intensity profile has width w 0 ) along the x direction: Since we are going to treat E 1 like a random number, we suppress the notion of the dependence of E 1 on the position x in the following.The (normalized) probability density distribution of this supposed random number E 1 in some interval x ∈ ½−mw 0 ; mw 0 is given by Here, θðξÞ denotes the Heaviside step function, which limits the range of possible values for E 1 .Now consider a specklelike field distribution E 2 ¼ E real þ iE imag , where the real and imaginary parts are uncorrelated and normally distributed random variables with standard deviation σ according to We are interested in the probability density distribution of the squared modulus of the sum of the amplitudes E 1 and E 2 , which (apart from some factors) corresponds to the intensity distribution of the coherent superposition; that is, The probability density function pðI=I 0 Þ of the intensity I of the total field as defined in Eq. (A1), with I 0 being the mean value over the considered interval, can then be calculated following these general rules for the manipulation of random numbers (see, for example, Ref. [62]).
( (3) In addition, we would like to consider the normalization of a distribution of random numbers with respect to their mean value.The mean value of a set of random numbers X is given by so that the probability density function of the normalized random number Y ¼ X=x 0 is calculated as p Y ðζÞ ¼ x 0 p X ðx 0 ζÞ: While there is no closed analytical solution for pðI=I 0 Þ, numerical results are presented in Fig. 13(a) for m ¼ 2.5 (which is used in the evaluation of the experimental data depicted throughout this article) and varying ratios A=σ.
The statistics of the fluctuating field in which we are interested (in this case a speckle distribution) is only revealed for coherent residuals of the input beam being A pronounced Gaussian background tends to suppress extreme events, if the profile is normalized with respect to its mean intensity I 0 in an interval x ∈ ½−2.5w 0 ; 2.5w 0 , where w 0 is (half) the width of the squared modulus of the Gaussian field profile.For an increasing ratio of the amplitude A of the Gaussian background and the standard deviation σ of the normally distributed random variables, the tail of the probability density drops down significantly already at small intensities.noticeably smaller than the field fluctuations.To some degree the suppression of high-amplitude events for ratios A=σ > 0 is a result of the intensity profile not dropping down to zero in the central region of the considered interval [see Fig. 13(b)].This affects the mean intensity I 0 , whose reciprocal value corresponds to the decay factor of the exponentially decaying probability density function of the intensity I of an ideal speckle pattern as stated by Eq. ( 1) represented by the dashed black line in Fig. 13(a).
Since the Gaussian background suppresses extreme events by trend, the frequency of extreme events visible in our experimental results might be underestimated.

APPENDIX B: NUMERICAL RESULTS OBTAINED WITH THE BEAM PROPAGATION METHOD
The experimental results presented in this article could not be reproduced with numerical schemes only considering forward propagation, like the beam propagation method, which is based on the discretized paraxial approximation of the Helmholtz equation.In Fig. 14 we have plotted some BPM results for a normally distributed refractive index landscape with correlation lengths corresponding to our experimental situation (l corr;∥ ¼ 20 μm in the direction parallel to the main propagation of light and l corr;⊥ ¼ 3 μm in the transverse direction) for varying strength of the refractive index modulations Δn.As can be seen, there is no significant abundance of extreme events relative to a speckle distribution (indicated by the dashed black line) even for strong variations of the refractive index.
For a weak modulation the central limit theorem does not yet hold as the output field still carries the imprint of the initial Gaussian beam.On the length scale of a few correlation lengths (for z < 1 mm) and for large index variations, deviations from speckle statistics are to be expected as a result of a slightly distorted phase front, effectively leading to focusing spatial regions and the consecutive formation of caustics.This phenomenon was reported in detail in Refs.[15,16].However, this is an effect that can be attributed solely to the spatial coherence of the impinging beam.It reaches its maximum a few times l scat behind the front facet, but then quickly fades away as soon as wave dynamical effects like random interference play an increasing role and dominate the intensity distribution to ultimately form a speckle pattern after all.The modulus of the current flowing in the negative z direction is given by integrating over the negative halfspace and multiplying the result by −1: According to our model there is no source at z ≤ 0 and all current flowing in the positive z direction is caused by a reflection of the current flowing in the negative direction; formally, J þz ð0Þ ≡ RJ −z ð0Þ: Inserting Eqs.(C1) and (C2), we find Ið0Þ − z 0 I 0 ð0Þ ¼ 0; with the extrapolation length Here, we introduced the effective reflectivity R for diffuse light, which must not be confused with the angle-(and polarization-)dependent Fresnel coefficient R Fresnel ðϑÞ for specular reflection.For R ¼ 0, we find z 0 ¼ ðπ=4Þl trans ≈ 0.785l trans , which is in very good agreement with the theoretical value 0.7104l trans found by calculating the Milne equation in three dimensions [39,40].
The effective reflectivity R takes into account the distribution of propagation directions of the diffusive light impinging on the interfaces.To derive z 0 from Fresnel coefficients, we calculate the reflected current flowing in the positive z direction as From J refl þz ð0Þ ≡ J þz ð0Þ, we finally find Equating the right-hand sides of Eqs.(C3) and (C4), we can readily calculate the effective reflectivity R as Since in our experiments we consider only the propagation of TE-polarized light (in terms of the waveguide structures), we have to take the Fresnel component R Fresnel;∥ for light polarized in the plane of incidence.Performing the integrals for n eff ¼ 3.33 (effective refractive index of the ground mode of the AlGaAs waveguide structure at λ 0 ¼ 1550 nm) and n env ¼ 1.00 (refractive index of air), we finally find R ¼ 0.75 and z 0 ≈ 5.5l trans .

APPENDIX D: ROBUSTNESS OF THE EFFECT TO LOSSES
The abundance of extreme events in the simulations and the shape of its accompanying probability density function are to some extent robust with regards to losses (absorption or, as in our case, scattering losses to the environment).In a lossy environment the buildup of high-amplitude events, visible as a clear deviation from the speckle distribution, takes a longer propagation length to reach its full extent.See Fig. 15 for a compilation of numerical results for different intrinsic losses.Figure 15(b) emulates the experimental results depicted in Fig. 2(a) even slightly better than the simulation displayed in Fig. 2(b), for which the attenuation length L loss as obtained from a fit to the total transmittance of the samples was used (see Sec. IV).This is in accordance with the discussion in Sec.V C, which is based on the fact that losses reduce the eigenvalues of the system in general.As a result, the pole of the function 1=ð1 − λÞ, which describes the amplitude of an eigenvector with corresponding eigenvalue λ in the steady state, has a less filtering effect.Since a multitude of eigenvectors equally contribute to the steady state, the central limit theorem kicks in, rendering the intensity distribution more specklelike.The fact that deviations from a speckle distribution can still be witnessed in Fig. 15 for long propagation lengths is most likely due to the iterative simulation having not quite reached the steady state yet.
FIG. 1. Experimental setup.(a) Waveguide structures consisting of a stacking of Al x Ga x−1 As layers with varying refractive index (depending on the content x of aluminum).(b) Because of thermal stress during fabrication, ripples form in a preferred direction.The red arrow indicates the main propagation direction of the guided light.Scale bar is 20 μm.(c) Light is coupled into the waveguide by focusing an elliptically shaped laser beam onto the front facet of the samples (see bottom left-hand inset).A beam splitter is used to monitor the front facet with a CCD camera (see top left-hand inset).Note the differing magnification factor in the horizontal direction due to the cylindrical lens.Besides the backreflected impinging beam, the backscattered light from inside the sample is clearly visible.The intensity distribution at the rear facet is monitored with another CCD camera (see bottom right-hand inset).

FIG. 2 .
FIG.2.Comparison of experimental and numerical results in the linear regime.(a) Histogram of the normalized intensity distribution as obtained experimentally from the rear facets of the samples.The shortest sample (L ¼ 4 mm) is unperturbed, all other samples are perturbed.For L ¼ 5 and 10 mm, the trend follows the exponentially decaying curve of a speckle distribution [dashed black line; compare with Eq. (1)].For longer propagation lengths (L ¼ 15, 20, 30 mm), extreme events are more probable than anticipated from a random distribution.(b) Histogram of normalized intensity distribution at different propagation lengths z from simulations based on the stray matrix method with plane wave excitation and weak, but prevailing backscattering.The exact scattering amplitudes are stated in the pictogram displayed in the bottom left-hand corner (clockwise, starting from the top): backward jrj, leftward jr ⊥ j, forward jtj, and rightward jr ⊥ j (see also Fig.5).For increasing z the probability distribution quickly surpasses a speckle distribution, before asymptotically approaching a heavy-tailed limiting function between z ¼ 20 and 30 mm.See also the inset for the dependence of the probability density at a fixed intensity value I ¼ 7.0I 0 on z.(c) Numerical results with less dominant backscattering.For intermediate propagation lengths the frequency of extreme events is greatly diminished.However, partial diffuse reflection at the rear facet increases the probability of extreme events in its vicinity (also visible in the inset).

FIG. 4 .
FIG. 4. Estimation of scattering amplitudes.(a) Microscope image of one of the samples illuminated from the back by an infrared diode.Inset: Filtered and extracted structure used in the FDTD simulation, where the light gray section is considered as a part of the background and the dark gray section constitutes the actual defect with varying refractive index.(b) Numerical result of the FDTD simulation for plane wave excitation (contour of defect plotted for reference) and index variation Δn ¼ 0.001.The absolute value of the total field is displayed in the inner rectangle, while the scattered field is shown in the outer frame (see respective color scale).Most of the radiation is scattered to the front.Only a small part is backscattered and the radiation scattered to the sides is negligible.(c) Plot of the angle-dependent scattering amplitudes for the index structure depicted in (a) for illumination with a plane wave polarized parallel to the plane of incidence (TE polarized in terms of the waveguide structure).
FIG.3.Sketch of the diffuse intensity.(a) A waveguide structure consisting of a guiding layer packed between two cladding layers.The mode profile in the y direction is localized in the guiding layer (light blue solid line).The excited mode defines an effective refractive index n eff .The waveguide is of length L in the z direction and practically infinite in the x direction.(b) Schematic normalized diffuse intensity for a source located at l trans ¼ 1.The length of the waveguide L is taken to be 10 and z 0 ¼ 1. Depicted are the cases for vanishing absorption (dark blue solid line) and L loss ¼ 3 (red solid line).(c) Fit (red solid line) to the experimentally determined mean total transmission TðLÞ ¼ P out =P in (gray squares) for cw operation as measured at the rear facets of the samples, based on Eq.(17).We obtain L loss ¼ 5.40 mm and l trans ¼ 1.87 mm.Note the logarithmic scale, which amplifies errors at low transmission values.The fit function is dashed for L < 2l trans , since Eq.(16) is reasonable only for L ≥ 2l trans .The light red corridor illustrates the functional prediction bounds for 1σ.A fit based on a model with nonnegligible losses along l trans (green solid line) does not result in a better fit.

FIG. 5 .
FIG.5.Numerical model.(a) Approximation of the sample as a network of scatterers.The initial field distribution is continuously fed to the first row as the field components propagate further and further through the network with every simulation step (see red arrows) until a steady state is reached.(b) One building block, i.e., a single scattering event, is implemented by a stray matrix m composed of the complex valued scattering amplitudes in the forward (t), backward (r), and transverse (r ⊥ ) directions.Illustrated is the effect on a discrete field component u in entering from the top.The scattered components are multiplied with a phase factor φ p (with p ¼ 1, 2, 3, 4) to account for the propagation between scattering sites.To simulate wave propagation in a randomly perturbed waveguide, these phases are picked from a Gaussian distribution with mean value kl scat .The effect of scattering and intermediate propagation is combined in the matrix Ŝ.The pictogram on the bottom right illustrates the short notion utilized in the other figures.(c) The transverse coordinate is denoted with n (with 1 ≤ n ≤ N W ) and the longitudinal coordinate is denoted with m (with 1 ≤ m ≤ N L ).

18 FIG. 7 .
FIG.7.Exemplary simulated intensity distributions inside the sample for an excitation with a Gaussian beam.The intensity is normalized row-wise and color scaling is the same for all images [note that in (b) the intensity has been divided by 2 and in (c) by 10 to comply with the color bar].(a) For vanishing standard deviation σ of the phases φ bond accumulated between scattering events, the system simulates homogeneous space so that the Gaussian beam diffracts normally.(b) For a slight variation of the phases the beam is disturbed during propagation, leading to the formation of large speckles.(c) For a stronger variation of the phases the intensity distribution qualitatively changes to a scarlike pattern, with single pixel values standing out up to a factor of 100 relative to the background (see inset).These randomly distributed resonant substructures are responsible for the deviation of the probability distribution from a speckle distribution.

50 FIG. 6 .
FIG. 6.Comparison between the simulation of isotropic and weak scattering.(a) Plot of the critical phase φ bond;crit , for which the dispersion relation is pointlike, depending on the scattering amplitude in the forward direction jtj under the assumption that jrj ¼ 2jr ⊥ j.The dispersion relation is asymptotically circular for φ bond ≲ φ bond;crit .(b) Exemplary dispersion relations for isotropic (jtj ¼ jrj ¼ jr ⊥ j ¼ 0.5, red line) and weak (jtj ¼ 0.9, jrj ¼ 2jr ⊥ j, green line) scattering.Depending on the exact value of φ bond , the two graphs are almost congruent in the Brillouin zone in terms of the transport mean free path length l trans ¼ Nl scat .(c),(d) Results of numerical simulations for isotropic (c) and weak (d) scattering.Weak scattering can be thought of as increasing the spatial resolution.The initial width of the Gaussian field distribution and the size of the considered region are the same in both cases in effective transverse and longitudinal distance measures x ¼ n=N and z ¼ m=N, respectively, which are the multiples of traversed transport mean free path lengths.

FIG. 8 .
FIG. 8. Field and intensity statistics of the cavity modes for the big domain (N L ¼ 1500, N W ¼ 501).(a) Probability density functions, i.e., relative frequencies of the real valued field amplitudes over the whole spatial domain for an ensemble of 256 eigenvectors and varying ratios between the backscattering amplitude jrj and the lateral scattering amplitude jr ⊥ j.For jrj ¼ jr ⊥ j, the statistics of the field amplitudes is Gaussian (blue solid line), for jrj ¼ 2jr ⊥ j it is sech shaped (red solid line), while for jrj ¼ 10jr ⊥ j a Laplace distribution evolves (green solid line).(b) Corresponding histogram of the normalized intensity distribution for the 256 eigenvectors evaluated at a propagation length z ¼ 10 mm [see dashed white lines in (c)-(f)].For prevailing backscattering the probability distribution increasingly differs from a speckle distribution (black dashed line).Exemplary plots of the intensity distribution of individual eigenvectors, for (c) jrj ¼ jr ⊥ j, (d) jrj ¼ 2jr ⊥ j, (e) jrj ¼ 3jr ⊥ j, and (f) jrj ¼ 10jr ⊥ j.Note that in (d)-(f) the intensities have been divided by 2, 3, and 20, respectively, to comply with the color bar.Scarlike patterns aligned along the propagation direction of light can be seen in all four cases, but the contrast to the background changes dramatically.

FIG. 9 .
FIG.9.Eigenvalues and evolution coefficients in the steady state for a small domain (N L ¼ 100, N W ¼ 31).(a) The eigenvalues λ (blue dots) are situated close to a circle of radius 1 in the complex plane, since the evolution matrix M is almost unitary.Additional losses reduce the radius of this circle (light blue dots in the inset).Red circles of increasing intensity illustrate increments of the function 1=ð1 − λÞ by an order of 10 around its pole.(b) Moduli of the amplitudes of the eigenvectors in the steady state versus the real part of the eigenvalues.For negligible losses, a few coefficients exceed all others (blue crosses); hence, few respective eigenvectors dominate the steady state.For increasing losses this effect is damped (light blue circles).(c) Comparison of the field statistics of the superposition in the steady state for negligible losses (blue squares, solid line is fit with sech function) and additional losses (light blue squares, solid line is fit with Gaussian).
FIG.10.Experimental results in the nonlinear regime.(a)-(c) Full histograms of the normalized intensity distribution as obtained at the rear facets of the samples for increasing pulse powers P pulse and L ¼ 5 mm (a), L ¼ 10 mm (b), and L ¼ 30 mm (c).In each case the respective linear measurement taken from Fig.2(a) is depicted for reference (dark blue dashed curve).(a) While for low pulse powers the intensity distribution is specklelike, increasing the nonlinearity induces more extreme events due to self-focusing.(b) After an initial rise of extreme events, further increasing the pulse power leads to a regression.(c) The abundance of extreme events at low pulse powers is progressively suppressed as the pulse power increases, finally restoring a specklelike distribution.(d) Comparison of the probability density of a specific intensity I ¼ 7.0I 0 for various sample lengths L. (e) Exemplary intensity profile as visible at the rear facet for L ¼ 30 mm. High pulse powers cause a thermalization of the intensity distribution by diminishing outstanding intensity peaks (see inset for a magnification).

5 FIG. 11 .FIG. 12 .
FIG. 11.Nonlinear light propagation in the unperturbed mesh.Panels (a)-(d) show the numerical outcomes with increasing focusing nonlinearity of the Kerr type from left to right.The diffraction of the initial Gaussian beam profile is progressively compensated for, up to a point where the profile comes close to that of a soliton being localized in the central region.

FIG. 13 .
FIG.13.Influence of a Gaussian background on the intensity statistics.(a) Probability density function pðI=I 0 Þ, whereI ¼ jE 1 þ E real þ iE imag j 2 .A pronounced Gaussian background tends to suppress extreme events, if the profile is normalized with respect to its mean intensity I 0 in an interval x ∈ ½−2.5w 0 ; 2.5w 0 , where w 0 is (half) the width of the squared modulus of the Gaussian field profile.For an increasing ratio of the amplitude A of the Gaussian background and the standard deviation σ of the normally distributed random variables, the tail of the probability density drops down significantly already at small intensities.(b)-(d) Exemplary numerical realizations of intensity profiles for A=σ ¼ 5.0 (b), A=σ ¼ 1.0 (c), and A=σ ¼ 0.1 (d).
FIG.13.Influence of a Gaussian background on the intensity statistics.(a) Probability density function pðI=I 0 Þ, whereI ¼ jE 1 þ E real þ iE imag j 2 .A pronounced Gaussian background tends to suppress extreme events, if the profile is normalized with respect to its mean intensity I 0 in an interval x ∈ ½−2.5w 0 ; 2.5w 0 , where w 0 is (half) the width of the squared modulus of the Gaussian field profile.For an increasing ratio of the amplitude A of the Gaussian background and the standard deviation σ of the normally distributed random variables, the tail of the probability density drops down significantly already at small intensities.(b)-(d) Exemplary numerical realizations of intensity profiles for A=σ ¼ 5.0 (b), A=σ ¼ 1.0 (c), and A=σ ¼ 0.1 (d).