---
abstract: 'The emergence of nonequilibrium phenomena in individual complex wave systems has long been of fundamental interests. Its analytic studies remain notoriously difficult. Using the mathematical tool of the [*concentration of measure*]{} (CM), we develop a theory for structures and fluctuations of waves in individual disordered media. We find that, for both diffusive and localized waves, fluctuations associated with the change in incoming waves (“wave-to-wave” fluctuations) exhibit a new kind of universalities, which does not exist in conventional mesoscopic fluctuations associated with the change in disorder realizations (“sample-to-sample” fluctuations), and originate from the coherence between the natural channels of waves – the [*transmission eigenchannels*]{}. Using the results obtained for wave-to-wave fluctuations, we find the criterion for almost all stationary scattering states to exhibit the same spatial structure such as the diffusive steady state. We further show that the expectations of observables at stationary scattering states are independent of incoming waves and given by their averages with respect to eigenchannels. This suggests the possibility of extending the studies of thermalization of closed systems to open systems, which provides new perspectives for the emergence of nonequilibrium statistical phenomena.'
author:
- Ping Fang
- Liyi Zhao
- Chushun Tian
title: 'Concentration-of-measure theory for structures and fluctuations of waves'
---

Recent studies on the foundations of equilibrium statistical mechanics [@Lebowitz06; @Popescu06; @Deutsch91; @Srednicki94; @Rigol08; @Rigol16; @Borhonovi16] have shed new light on the long-standing problem of nonequilibrium phenomena in individual (quantum) wave systems, where neither fictitious ensembles nor reservoirs exist [@von; @Neumann29]. Many of them [@Deutsch91; @Srednicki94] rely on a conjecture of Berry, i.e., in closed systems random scattering can render waves structureless on large spatial scales [@Berry77]. This wave property gives rise to a basic feature of thermal equilibrium phenomena in individual closed systems, i.e., spatial homogeneity. Yet, a major topic of nonequilibrium statistical mechanics is concerned with various spatial structures in open systems [@Reichl; @Dorfman; @Haenggi05; @Lepri03; @Gaspard96]. This sharp contrast motivates exploring in depth spatial structures of waves, i.e., scattering states, in individual open systems, which may also pave a way for extending the studies of the relations between spatial and entanglement structures in an ensemble of open disordered systems – a new aspect of the fundamentals of nonequilibrium statistical mechanics [@Huse18] – to an individual member.

In fact, spatial structures and fluctuations of waves in open disordered media are central topics of mesoscopic physics [@Rossum99; @Sheng95; @Mirlin00; @Efetov97; @Akkermans]. However, most theoretical efforts have been focused on disorder ensembles; the significance of waves in individual disordered media has been emphasized only recently [@Dogariu09; @Genack13]. The common wisdom of using self-averaging or the ergodic hypothesis to connect certain properties of individual disordered media to their disorder averages [@Hanggi98; @Lifshits88] essentially requires the thermodynamic limit, and cannot be applied to study fluctuations in mesoscopic scales. It remains a challenge to construct a theory for wave statistics in individual mesoscopic systems, where rich fluctuation phenomena of the wave origin can be driven, e.g., by changing the incoming wave. Such wave-to-wave fluctuations differ from well understood sample-to-sample fluctuations [@Altshuler91; @Rossum99; @Sheng95; @Beenakker97; @Hanggi98; @Lifshits88; @Altshuler85; @Lee85; @Mirlin00; @Kamenev; @Efetov97; @Akkermans]. Their in-depth studies are of both fundamental and practical importance. Indeed, in individual mesoscopic systems fluctuations and irreversibility have been known to be closely related [@Altshuler96]. In addition, the strength of wave-to-wave fluctuations determines whether a generic wave can represent the behaviors of most waves, i.e., is typical: in recent years the arising of [*typicality*]{} has appeared as a central topic in statistical mechanics [@von; @Neumann29; @Lebowitz06; @Popescu06], but been restricted to closed systems so far. On the other hand, wave statistics in individual open disordered media has found many optical applications [@Gigan10; @Fink12; @Choi12; @Rotter17].

Recently, the CM [@Milman86; @Ledoux01; @Boucheron13] – “one of the great ideas of analysis in our time” [@Talagrand96] – has been adopted to study statistical phenomena in individual closed classical [@Bievre17] and quantum [@Popescu06; @Winter06; @Guarneri15] systems. The CM is rooted in high-dimensional geometry. The idea can be illustrated by the unit sphere, for which the area of the sphere becomes more and more concentrated around the equator as the dimension increases. Eventually in high dimensions the entire area almost concentrates around the equator. This property can then be visualized by real-valued functions over the sphere with nice continuity properties, through their concentration around some constant value. When the sphere is replaced by a general high-dimensional geometric body (e.g., the Euclidean space) and the area measure by others (e.g., the Gaussian measure), similar results follow. This idea opens new perspectives of probability theory [@Talagrand96; @Boucheron13; @Ledoux01]. It allows us not only to study variables with complicated dependence on random variables instead of being their sum, but also to obtain results which are nonasymptotic, i.e., do not require the limit of large number of variables. A detailed introduction of CM is given in section S0 of supplemental materials (SM) [@SM].

In this work we employ CM to explore universal statistical phenomena of waves in individual open disordered media. We launch a classical wave of circular frequency $\Omega$ and carrying unit energy flux into a finite medium with $N$ ($=$$\frac{\Omega}{\pi}$$\times$the width) channels and length $L$ [@note3; @note5] (Fig. \[fig:2\]). Keeping the disorder realization fixed, but allowing the incoming wave to vary gives rise to various wave-to-wave fluctuations. Below, most attentions are paid to the fluctuations of the spatial structure of scattering state corresponding to the incoming current amplitude $c$, i.e., the depth ($x$) profile $I_x(c)$ of energy density integrated over the cross section. We develop a CM theory of wave-to-wave fluctuations. Physically, it provides information on a single stationary scattering state in a single disordered medium and the differences of behaviors between a disorder ensemble and an individual member; technically, instead of traditional impurity diagrams [@Rossum99; @Sheng95; @Akkermans] and field theories [@Mirlin00; @Kamenev; @Efetov97], its key components are various [*concentration inequalities*]{} [@Boucheron13] of observables \[e.g., $I_x(c)$\]. Armed with the developed theory we achieve the following results:

-   We find that compared to conventional sample-to-sample fluctuations, wave-to-wave fluctuations exhibit a number of “anomalies”. In particular, irrespective of regimes of wave propagation (diffusive, localized, etc.), the distribution of $I_x(c)$ is always [*sub-Gaussian*]{}, i.e., has an (upper) tail decaying at least as fast as a Gaussian tail \[Eq. (\[eq:1\])\]. Contrary to this, for sample-to-sample fluctuations of observables such as total transmission, as waves are more and more localized the distribution tail decays slower and slower, and the shape of the tail changes dramatically [@Mirlin00; @Genack97].

-   Furthermore, we find that the wave-to-wave fluctuations of $I_x(c)$ are governed by an $x$-dependent curve $\|I_x\|_{\rm Lip}$, which arises from the phase coherence between distinct eigenchannels – the natural channels for wave propagation in disordered media [@Choi11; @Tian15]. In contrast, the sample-to-sample fluctuations of $I_x(c)$ ($c$ fixed) are governed by the conductance [@Rossum95] known to equal the number of open eigenchannels [@Imry86]. For diffusive waves we find that the curve is universal with respect to disorder realizations $\omega$ at large $N$ (cf. Fig. \[fig:4\]).

-   We find the criterion (\[eq:7\]) for almost all stationary scattering states to exhibit the same spatial structure, i.e., a nonequilibrium steady state, and show that it can be readily satisfied for diffusive waves.

-   We show that the expectations of generic observables at stationary scattering states are independent of incoming waves and given by their averages with respect to eigenchannels \[Eq. (\[eq:13\])\], and find the corresponding criterion.

The results summarized above are attributed to general wave properties, and thus apply to both classical and quantum waves. (This is similar to Anderson localization applying to both classical and quantum waves [@Mirlin00; @Sheng95].)

While we are not aware of any studies of the applications of CM to wave propagation and scattering in disordered media, we begin with a general discussion on how the high-dimensional geometry emerges from the present setting (Fig. \[fig:2\]), and further provides a basis for applying CM. For simplicity we consider a two-dimensional ($2$D) medium. A disordered dielectric configuration $\delta\epsilon (x,y)$ is embedded into the air background. So the wave field $E(x,y)$ satisfies the Helmholtz equation [@note5], $$\begin{aligned}
\label{eq:15}
    \left\{\partial_x^2 +\partial_y^2 +\Omega^2 [1+
    %\chi_{[0,L]}(x)
    \delta\epsilon (x,y)]\right\}E(x,y)=0.\end{aligned}$$ Given $N$ channel bases, an incoming current amplitude is a projection represented by $N$ complex coefficients: $(c_1$,$c_2$,...,$c_N)$$\equiv$$c$. As the incoming wave carries a unit energy flux, we have $\sum_{n=1}^N$$|$$c_n$$|^2$$=$$1$. So all $c$ constitute a high-dimensional geometric body, the unit sphere $S^{2N-1}$. Next, we discretize the medium into a lattice of $M$ points. The $M$ values $\{\omega_{(x,y)}$$\equiv$$-$$\Omega^2$$\delta$$\epsilon$$($$x$,$y$)$\}$ then constitute the coordinate of another high-dimensional geometric body, the Euclidean space $\mathds{R}^M$, a point in which corresponds to a disorder realization $\omega$$\equiv$$\{$$\omega_{(x,y)}$$\}$. As observables depending on $c$ (respectively $\omega$) define real-valued functions over $S^{2N-1}$ (respectively $\mathds{R}^M$), we can apply CM to them and probe the properties of the wave-to-wave (respectively sample-to-sample) fluctuations by corresponding observables.

![A wave is launched into a disordered medium of $N$ channels. In this setting there are two geometric bodies that provide the basis for applying CM, namely, the sphere $S^{2N-1}$ constituted by distinct incoming current amplitudes $c$ and the Euclidean space $\mathds{R}^M$ by distinct disorder realizations $\omega$.[]{data-label="fig:2"}](typicality_propagation.pdf){width="8.7cm"}

[*Construction of the theory.*]{} We divide the construction into five steps so that the readers may keep track of it. In each step, we outline the derivations and present the key results; motivations and (or) physical implications are also discussed. Technical details and expanded discussions are relegated to the self-contained SM.

[*Step 1 – formulation of the problem.*]{} Below we choose the eigenchannels as the basis. Thus we first introduce the eigenchannel briefly, and this is done in three steps as follows [@Tian15]. (i) Consider the transmission matrix $t\equiv \{t_{ab}\}$, which transforms an incoming current amplitude $c$ into a transmitted current amplitude given by $tc$. Both current amplitude are vectors in the space spanned by the ideal waveguide modes $\varphi_{a}(y)$, with $a$ the mode index. The matrix element $t_{ab}=-i\sqrt{\tilde v_a\tilde v_b} \langle
x=\infty a|G
%(\omega)
|x'=-\infty b\rangle
%\int\!\!\!\!\int dydy' \varphi_a(y)\varphi_b^*(y')G(x=\infty,y,x'=-\infty,y')
$ [@note6], where $G$ is the retarded Green’s function and $\tilde v_a$ is the group velocity of mode $a$. By the singular value decomposition, $t=\sum_{n=1}^N u_n \sqrt{\tau_n} v_n^\dagger$, we obtain a transmission eigenvalue spectrum $\{\tau_n\}$ ($\tau_n$ decreases with $n$.) and two mutual orthogonal unit vectors $\{u_n\}
%\equiv\{u_{an}\}
$ and $\{v_n\}
%\equiv\{v_{an}\}
$. (ii) Replacing $x=\infty$ in $t_{ab}$ above by arbitrary $x\in [0,L]$, we make the extension: $t\rightarrow t(x)\equiv \{t_{ab}(x)\}$. This gives the vector field inside the medium, $E_{n}(x)
%\equiv \{E_{an}(x)\}
=t(x)v_n$, excited by $v_n$. (iii) The triple: $(\tau_n,v_n, E_n(x)\equiv\{E_{na}(x)\})$ defines an eigenchannel, which are completely fixed by $\omega$ and $\Omega$ [@note7]. Each channel transmits waves with transmission coefficient $\tau_n$, and has a specific $2$D spatial structure, namely, the energy density profile $|E_{n}(x,y)|^2$ with $E_n(x,y)=\sum_{a=1}^N E_{na}(x)\varphi_a^*(y)$. Integrating out $y$ we reduce $|E_{n}(x,y)|^2$ to a one-dimensional ($1$D) structure, $$\label{eq:16}
\begin{array}{c}
  W_{\tau_n}(x)\equiv \int dy|E_{n}(x,y)|^2=E_{n}^\dagger(x)\cdot E_{n}(x),
\end{array}$$ where $\cdot$ refers to a scalar product.

To proceed we introduce the precise definition of $I_x(c)$ [@note2]: $I_x(c)\equiv \int dy E^*(x,y)\hat v_xE(x,y)$, where $\hat v_x\equiv (1-\partial_{\Omega y}^2)^{\frac{1}{2}}$ is a scalar (not vector) operator accounting for the absolute value of the group velocity in waveguide modes. Treating $\omega_{(x,y)}$ as a scattering potential, we apply the scattering theory of waves [@Newton82] to Eq. (\[eq:15\]) and find $\hat v_x^{\frac{1}{2}}E(x,y)=\sum_{a=1}^N(t(x)c)_a\varphi^*_a(y)$. Here $c=\sum_{n=1}^N c_n v_n$ is the incoming current amplitude in the eigenchannel representation. Thus $$\begin{aligned}
\label{eq:10}
\begin{array}{c}
  I_x(c)=\sum_{n,n'=1}^N c_n^*c_{n'} E_{n}^\dagger(x)\cdot E_{n'}(x).
\end{array}\end{aligned}$$ This defines a family of real-valued functions over $S^{2N-1}$, and $x$ labels these functions. Note that at $x<L$ the vectors $E_{n}(x)$ are not orthogonal.

Then the problem is: [*For fixed $\omega$, does $I_x(c)$ exhibit universal behaviors when $c$ varies*]{}? A natural idea is to calculate all the cumulants of $I_x(c)$ and find the distribution. But, one then needs to calculate an infinite number of products of $E_{n}^\dagger\cdot E_{n'}$, and sum up their contributions, which is a formidable task especially for small $N$. The CM allows a different route, which we follow below.

[*Step 2 – Lipschitz continuity: a building block of CM.*]{} This is the concept that formalizes the “nice continuity properties” of real-valued functions mentioned in the introductory part. Let a generic space $\mathscr{C}$ be equipped with the Euclidean metric $\|\cdot\|$. For $f$$:$$\mathscr{C}$$\rightarrow$$\mathds{R}$, if $$\begin{aligned}
\begin{array}{c}
  \|f\|_{{\rm Lip}}\equiv{\rm sup}_{z,z'
  %\in \mathscr{C}
  }\frac{|f(z)-f(z')|}{\|z-z'\|}<\infty\\
  \Leftrightarrow \,|f(z)-f(z')|\leq \|f\|_{{\rm Lip}}\,\|z-z'\|,
\end{array}
\label{eq:11}
    %\|f\|_{{\rm Lip}}\equiv{\rm sup}_{c,c'\in \mathscr{C}}\frac{|f(c)-f(c')|}{\|c-c'\|}<\infty \label{eq:11}\\
    %\Leftrightarrow \,|f(c)-f(c')|\leq \|f\|_{{\rm Lip}}\,\|c-c'\|,\quad\nonumber\end{aligned}$$ where ‘${\rm sup}$’ stands for the least upper bound, then $f(z)$ is said to have the Lipschitz continuity or be Lipschitz, and $\|f\|_{{\rm Lip}}$ is called the Lipschitz constant. As we will see below, even though $f$ has a very complicated dependence on $c$, its wave-to-wave fluctuations are controlled by a single parameter, i.e., $\|f\|_{{\rm Lip}}$.

[*Step 3 – Concentration inequality for $I_x(c)$ and results for general waves.*]{} With the preparations above we are now ready to introduce the following result of CM:\
\
L$\acute{\rm e}$vy’s lemma [@Winter06; @Milman86]. [*Let $\mu$ be the uniform probability measure over $S^{2N-1}$, and $f: S^{2N-1}\rightarrow \mathds{R}$ be Lipschitz. Then the probability for the deviation between $f$ and its mean $\int f d\mu$ to exceed $\varepsilon$ is $\leq 2e^{-
\frac{\delta\varepsilon^2 N}{\|f\|_{{\rm Lip}}^2}}$, where $\delta$ is some positive absolute constant.*]{}\
This means that $f$ concentrates around $\int fd\mu$ with a rate increasing rapidly with $\frac{N}{\|f\|_{{\rm Lip}}^2}$. We stress that $\|f\|_{{\rm Lip}}$ depends on $N$ generally. By the lemma the distribution of $f$ is sub-Gaussian. But, unlike the central limit theorem, the lemma does not require the large $N$ limit, i.e., is nonasymptotic; instead, $N$ can be very small (cf. Fig. \[fig:3\]).

To apply L$\acute{\rm e}$vy’s lemma to $I_x(c)$, in SM (S1) we use Eq. (\[eq:10\]) to derive an analytic expression of the Lipschitz constant $\|I_x\|_{{\rm Lip}}$ of $I_x(c)$. The result reads $$\begin{aligned}
\begin{array}{c}
  \|I_x\|_{{\rm Lip}}={\rm sup}_c L_x(c;\omega),\\
  L_x=\pi \left(\sum_{n=1}^N|\sum_{n'=1}^N \!c_{n'}^* E_{n'}^\dagger(x)\cdot E_{n}(x)|^2-I_x^2\right)^{1/2},
\end{array}
\label{eq:4}\end{aligned}$$ where $L_x$ depends on $c,\omega$ in general. According to Eq. (\[eq:4\]) $\|I_x\|_{{\rm Lip}}<\infty$. Combined with L$\acute{\rm e}$vy’s lemma this gives the following concentration inequality, $$\label{eq:1}
\begin{array}{c}
  {\rm Pr}\left(|I_x(c)-W(x;\omega)|>\varepsilon\right)\leq 2e^{-
    %\frac{2}{9\pi^3\ln 2}
    \frac{\delta\varepsilon^2 N}{\|I_x\|_{{\rm Lip}}^2}}.
\end{array}$$ Here $W(x;\omega)\equiv \int I_x(c)d\mu$ and “Pr” stands for probability. After simple algebra we reduce $W(x;\omega)$ to $$\begin{aligned}
\begin{array}{c}
  W(x;\omega)=\frac{1}{N}\sum_{n=1}^N W_{\tau_n}(x).
\end{array}
\label{eq:12}\end{aligned}$$ The results (\[eq:4\])-(\[eq:12\]) hold for general $N, \omega$, regardless of regimes of wave propagation (diffusive, localized, etc.).

![Using Eqs. (\[eq:4\]) and (\[eq:3\]), we calculate $L_x(c;\omega)$ at $N=800$ for $4$ randomly chosen $c$ at fixed $\omega$ and for $4$ randomly chosen $\omega$ at fixed $c$, respectively, and calculate $\|I_x\|_{\rm Lip}$ for $3$ large $N$. All profiles collapse into the a single curve. $L=50$[]{data-label="fig:4"}](Lipschitz_N_disorder.pdf){width="8.7cm"}

From the inequality (\[eq:1\]) we see that provided $$\label{eq:7}
    W(x;\omega)\gg \|I_x\|_{{\rm Lip}}/\sqrt{N},$$ the wave-to-wave fluctuations of $I_x(c)$ are negligible and almost all incoming waves behave in essentially the same way: their energies are stored in distinct channels with an equal weight of $1/N$ and a nonequilibrium steady state universal with respect to $c$ namely $W(x;\omega)$ results, i.e., $I_x$$($$c$$)$$\approx$$W$$($$x$$;$$\omega$$)$. The state does not carry information on the phase coherence between eigenchannels. The phase coherence, as shown by the sub-Gaussian tail in (\[eq:1\]) and Eq. (\[eq:4\]), enters into $\|I_x\|_{\rm Lip}$ and influences strongly wave-to-wave fluctuations (see [*Step 5*]{} for further discussions).

From the inequality (\[eq:1\]) we also see that the distribution tail of $I_x(c)$ decays at least as fast as a Gaussian tail. In contrast, in sample-to-sample fluctuations ($c$ fixed) the distribution tail decays much slower, known (for $x=L$) to be exponential for diffusive waves [@Rossum95] and log-normal for deeply localized waves [@Beenakker97; @Mirlin00]. [*Step 4 – Diffusive steady state and its fluctuations.*]{} Below we use the general results (\[eq:4\])-(\[eq:12\]) to explore in-depth diffusive waves. Numerical studies have shown that in quasi $1$D media [@note8] the disorder average of Eq. (\[eq:12\]) gives a diffusive steady state [@Tian15], but it is difficult to (dis)prove analytically that without the averaging, this remains true for general geometry. Indeed, for large $N$ ($L$ and $\Omega$ fixed) the medium is a slab [@note8] and thus high-dimensional, but in high dimension the explicit form of $W_{\tau_n}(x)$, even its disorder average, is unknown; for small $N$, i.e., a short quasi $1$D medium [@note8], the impacts of the sample-to-sample fluctuations on $W_{\tau_n}(x)$ have not yet been studied.

To study Eq. (\[eq:12\]) we start from large $N$ and establish a concentration inequality of $W(x;\omega)$. To this end we show in SM (S2) that, even for a single disordered slab, distinct eigenchannel structures $W_{\tau_n}(x)$ are described by a single formula [@Tian15], that depends smoothly on the eigenvalue $\tau_n$ and was derived originally for an ensemble of quasi $1$D disordered media. Using this fact, we show in SM (S3) that $W(x;\omega)$, which is a real-valued function over $\mathds{R}^M\ni\omega$, is Lipschitz, i.e., $$\begin{aligned}
\label{eq:2}
  |W(x;\omega)-W(x;\omega')|\leq
  %\|W_x\|_{{\rm Lip}}
  \tilde c(x)N^{-\frac{1}{2}}\, \|\omega-\omega'\|.\end{aligned}$$ Here $\tilde c(x)={\cal O}(1)$; its explicit form is unimportant and given in SM. The property (\[eq:2\]) allows us to use Pisier’s theorem in CM [@Milman86] to show \[SM (S4)\] the following:\
\
[*If $\omega=\{\omega_{(x,y)}\}$ is drawn randomly from an ensemble of disorder realizations with $\omega_{(x,y)}$ being independent Gaussian variables of zero mean and variance $\sigma^2$, then the probability for the deviation between $W(x;\omega)$ and its disorder mean ${\rm E}[W(x;\omega)]$ to exceed $\varepsilon$ satisfies*]{} $$\label{eq:5}
    {\rm Pr}\left(|W(x;\omega)-{\rm E}[W(x;\omega)]|>\varepsilon\right)\leq 2e^{-\frac{2N\varepsilon^2}{(\tilde c(x)\pi\sigma)^2}}.$$ Thus the concentration is strong for large $N$, i.e., $$\begin{aligned}
W(x;\omega)\approx {\rm E}[W(x;\omega)],\, for\,\, almost\,\, all\,\, \omega.
\label{eq:6}\end{aligned}$$ It is important to remark that the factor $N$ in the sub-Gaussian bound of the concentration inequality (\[eq:5\]) comes from the Lipschitz constant of $W(x;\omega)$, i.e., the coefficient on the right-hand side of the inequality (\[eq:2\]).

![Simulations show that in a single slab $\omega$ the profiles $I_x(c)$ for distinct $c$ concentrate around $W(x;\omega)$ (lower panels), and the data: $(\|I_x\|_{\rm Lip},\sqrt{\pi^2Nvar(I_x)})$ (symbols) collapse into a straight line of unit slope for distinct large $N$ while deviate from this line for small $N$ (upper panel). $L=50$[]{data-label="fig:1"}](typicality_steady_state.pdf){width="8.7cm"}

To proceed to small $N$ we observe that the detailed structures of $\{W_{\tau_n}(x)\}$ enter into the inequality (\[eq:5\]) only through the unimportant factor of $\tilde c(x)$. Thus we conjecture that the inequality applies in this case also. While its rigorous proof is beyond this work, we confirmed the conjecture numerically for $N$ being as small as $20$ \[SM (S4)\]. So Eq. (\[eq:6\]) holds for both large and small $N$.

Due to $I_x(c)\approx W(x;\omega)$ the disorder average of $I_x(c)$ ($c$ fixed) ${\rm E}[I_x(c)]\approx {\rm E}[W(x;\omega)]$. Together with Eq. (\[eq:6\]) this gives $W(x;\omega)\approx {\rm E}[I_x(c)]$. As ${\rm E}[I_x(c)]$ is known to be the solution to the diffusion equation [@Sheng95; @Rossum99; @Akkermans], $W(x;\omega)$ is a diffusive steady state (for almost all $\omega$), which decreases linearly in $x$. This result renders the transport mean free path $\ell$ well defined for single $\omega$ – because for a diffusive steady state the total transmission $W(L;\omega)=\ell/L$ [@Rossum99] – and identical to that defined for a disorder ensemble.

[*Step 5 – Wave-to-wave fluctuations in diffusive regime.*]{} According to the criterion (\[eq:7\]), weak wave-to-wave fluctuations ensure the emergence of a steady state universal with respect to $c$. To study these fluctuations we need to better understand $\|I_x\|_{\rm Lip}$. Let us start from large $N$. In SM (S1) we calculate Eq. (\[eq:4\]) and obtain for this case $$\begin{aligned}
\label{eq:3}
\begin{array}{c}
  \|I_x\|_{{\rm Lip}}=\int L_x(c;\omega) d\mu={\cal O}(1).
\end{array}\end{aligned}$$ The profiles of $\|I_x\|_{{\rm Lip}}$ and $L_x$ are shown in Fig. \[fig:4\]. Surprisingly, we find that the profiles of $\|I_x\|_{{\rm Lip}}$ at distinct $N$ collapse into a single curve; moreover, the profile of $L_x$ is universal with respect to $c$ and $\omega$, and the universal curve is identical to that of $\|I_x\|_{{\rm Lip}}$.

This universality of $L_x$, together with Eq. (\[eq:3\]), implies the universality of $\|I_x\|_{{\rm Lip}}$ with respect to $\omega$. As shown in SM (S1), it even leads to an explicit expression of $\|I_x\|_{{\rm Lip}}^2$: $$\begin{aligned}
\label{eq:14}
\begin{array}{c}
  \|I_x\|_{\rm Lip}^2=\frac{\pi^2}{N+1}\big(\sum_{n=1}^N\left(W_{\tau_n}(x)-W(x;\omega)\right)^2\\
  +\sum_{n\neq n'}^N|\sum_{a=1}^NE_{na}^*(x)E_{n'a}(x)|^2\big).
\end{array}\end{aligned}$$ By using Eq. (\[eq:14\]) we find in SM (S5) that $$\label{eq:20}
    var(I_x)=\|I_x\|_{\rm Lip}^2/(\pi^2 N),$$ and thus includes both incoherent and coherent contributions of eigenchannels, corresponding respectively to the first and second term in Eq. (\[eq:14\]). The incoherent (coherent) contribution dominates the back (front) part of the medium. At the output end, i.e., $x=L$, the second term in Eq. (\[eq:14\]) vanishes and the fluctuations include the incoherent contribution only.

For small $N$ the universality above is violated. But numerical calculations show that $\|I_x\|_{{\rm Lip}}={\cal O}(1)$ still holds (see the symbols corresponding to $N=20$ in the upper panel of Fig. \[fig:1\]). Due to this and $W(x;\omega)={\cal O}(1)$, for diffusive waves the criterion (\[eq:7\]) can be readily satisfied.

![Quasi $1$D simulations show that, for both diffusive (a) and localized (b) waves, the distribution of wave-to-wave fluctuations of $I_x(c)$ (black histograms) displays a tail well fit by a Gaussian distribution (pink dashed lines), while that of sample-to-sample fluctuations (blue histograms) is exponential (a) and stretched-exponential with a stretching exponent of $0.4$ (b), respectively (red dashed lines). In (b) the main panel and inset differ in the horizontal axis. The localization length $N\ell$ [@Dorokhov84; @Mello88] is $130$ in (a) and $65$ in (b). ${\bar I}_x$$=$$\int$$I_xp$$($$I_x$$)$$d$$I_x$.[]{data-label="fig:3"}](typicality_fluctuation.pdf){width="8.7cm"}

[*Numerical confirmations.*]{} We put the theory into numerical tests. The methods of numerical experiments are described in SM (S7). First, we simulate wave propagation in a single slab for $10^4$ randomly chosen $c$. Simulations confirm that the profiles $I_x(c)$ concentrate around a linear decrease for both large and small $N$ (Fig. \[fig:1\], lower panels). This also gives $\ell=13$ for single $\omega$ [@note3]. Moreover, from the distribution $p(I_x)$ of wave-to-wave fluctuations we compute $var(I_x)$, and find that the relation (\[eq:20\]) holds for large $N$ but is violated for small $N$ (Fig. \[fig:1\], upper panel), as expected by our theory. Secondly, we simulate propagation of diffusive and localized waves in quasi $1$D media. We perform the statistics of wave-to-wave and sample-to-sample fluctuations of $I_x$: for the former we fix the disorder realization $\omega$ and randomly choose $10^5$ incoming current amplitudes $c$ and for the latter do the opposite. As shown in Fig. \[fig:3\], for both diffusive and localized waves the distribution of wave-to-wave fluctuations displays a Gaussian tail, in agreement with the inequality (\[eq:1\]). In contrast, the distribution of sample-to-sample fluctuations is much broader, which is exponential for diffusive waves and stretched-exponential for localized waves [@note1].

Our theory provides new perspectives for the long-standing problem of the emergence of irreversibility in individual systems. In particular, in SM (S6) we consider a generic hermitian (which is not essential) operator $\hat O$, and study its expectation value at the stationary scattering state $E$$($$x$,$y$$)$ (determined by the incoming wave $c$) is $O(c)$$\equiv$$\langle E|\hat O|E\rangle$$=$$\sum_{n,n'}$$c_n^*$$c_{n'}$$\langle E_n|\hat O|E_{n'}\rangle$. Recall that $E_n$ is the wave field $E_n$$($$x$,$y$$)$ of the $n$th eigenchannel. Repeating the analysis above we find $$\label{eq:13}
\begin{array}{c}
  O(c)\!\approx\!\frac{1}{N}\sum_{n=1}^N \langle E_n|\hat O|E_n\rangle\equiv {\bar O},\, for\, almost\, all\,c,
\end{array}$$ if ${\bar O}$$\gg$$\frac{\|O\|_{{\rm Lip}}}{\sqrt{N}}$. This incoming wave-independence of observables resembles thermalization in closed systems [@Lebowitz06; @Popescu06; @Deutsch91; @Srednicki94; @Rigol08; @Rigol16; @Borhonovi16; @von; @Neumann29]. But, as the systems here are open conceptual differences exist. Notably, bound states and equilibrium thermal ensembles in closed systems are replaced respectively by stationary scattering states and $\frac{1}{N}$$\sum_{n}$$|E_n$$\rangle$$\langle$$E_n|$, which may be called the eigenchannel ensemble. In the future it is interesting to generalize Eq. (\[eq:13\]) to many-body systems, which would enable us to explore the relations between diffusive steady states and entanglement structures, but, unlike Ref. [@Huse18], requiring neither reservoirs nor disorder ensembles.

C.T. is grateful to J.-C. Garreau, F.-L. Lin, and Z. Q. Zhang, especially to A. Z. Genack and I. Guarneri for inspiring discussions and comments on the manuscript. This work is supported by the National Natural Science Foundation of China (No. 11535011 and No. 11747601).

S. Popescu, A. J. Short, and, A. Winters, [ Nat. Phys.]{} **2**, 754 (2006).

S. Goldstein, J. L. Lebowitz, R. Tumulka, and N. Zanghi, [ Phys. Rev. Lett.]{} **96**, 050403 (2006).

J. Deutsch, [ Phys. Rev. A]{} **43**, 2046 (1991).

M. Srednicki, [ Phys. Rev. E]{} **50**, 888 (1994).

M. Rigol, V. Dunjko, and M. Olshanii, [Nature]{} **452**, 854 (2008).

L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, [ Adv. Phys.]{} **65**, 239 (2016).

F. Borgonovi, F. M. Izrailev, L. F. Santos, and V. G. Zelevinsky, Phys. Rep. **626**, 1 (2016).

J. von Neumann, Z. Phys. **57**, 30 (1929).

M. V. Berry, J. Phys.: Math. Gen. A **10**, 2083 (1977).

L. E. Reichl, [*A Modern Course in Statistical Physics*]{} 4th Ed. (Wiley-VCH, Weinheim, Germany, 2016).

J. R. Dorfman, [*An Introduction to Chaos in Nonequilibrium Statistical Mechanics*]{} (Cambridge University Press, Cambridge, UK, 1999).

, P. H$\ddot{\rm a}$nggi, J. Luzka, and P. Talkner eds., New J. Phys. **7**, (2005).

S. Lepri, R. Livi, and A. Politi, [ Phys. Rep.]{} **377**, 1 (2003).

P. Gaspard, [ Phys. Rev. E]{} **53**, 4379 (1996).

M. J. Gullans and D. A. Huse, arXiv: 1804.00010.

M. C. W. van Rossum and Th. M. Nieuwenhuizen, [Rev. Mod. Phys.]{} **71**, 313 (1999).

P. Sheng, [*Introduction to Wave Scattering, Localization, and Mesoscopic Phenomena*]{} 2nd Ed. (Springer, Berlin, Germany, 2006).

A. D. Mirlin, [ Phys. Rep.]{} **326**, 259 (2000).

K. B. Efetov, [*Supersymmetry in Disorder and Chaos*]{} (Cambridge University, Cambridge, UK, 1997).

E. Akkermans and G. Montambeaux, [*Mesoscopic Physics of Electrons and Photons*]{} (Cambridge University Press, UK, 2007).

J. Broky, K. M. Douglass, J. Ellis, and A. Dogariu, Opt. Express **17**, 10466 (2009).

M. Davy, Z. Shi, J. Wang, and A. Z. Genack, Opt. Express **21**, 10367 (2013).

T. Dittrich, P. H$\ddot{\rm a}$nggi, G.-L. Ingold, B. Kramer, B., G. Sch$\ddot{\rm o}$n, and W. Zwerger, [*Quantum Transport and Dissipation*]{} (Wiley-VCH, Weinheim, Germany, 1998).

I. M. Lifshits, S. A. Gredeskul, and L. A. Pastur, [*Introduction to The Theory of Disordered Systems*]{} (Wiley, New York, 1988).

B. L. Altshuler, Pis’ma Zh. E’ksp. Teor. Fiz. **41**, 530 (1985) \[JETP Lett. **41**, 648 (1985)\].

P. A. Lee and A. D. Stone, Phys. Rev. Lett. **55**, 1622 (1985).

B. L. Altshuler, P. A. Lee, and R. A. Webb, Eds. (North-Holland, Amsterdam, 1991).

A. Kamenev, [*Field Theory of Non-Equilibrium Systems*]{} (Cambridge University Press, Cambridge, UK, 2011).

C. W. J. Beenakker, Rev. Mod. Phys. **69**, 731 (1997).

A. V. Andreev, O. Agam, B. D. Simons, and B. L. Altshuler, Phys. Rev. Lett. **76**, 3947 (1996).

S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, Phys. Rev. Lett. **104**, 100601 (2010).

A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, Nat. Photonics **6**, 283 (2012).

Y. Choi, C. Yoon, M. Kim, T. D. Yang, C. Fang-Yen, R. R. Dasari, K. J. Lee, and W. Choi, Phys. Rev. Lett. **109**, 203901 (2012).

S. Rotter and S. Gigan, [ Rev. Mod. Phys.]{} **89**, 015005 (2017).

V. D. Milman and G. Schechtman, [*Asymptotic Theory of Finite-Dimensional Normed Spaces*]{} (Springer, Berlin, Germany, 1986). M. Ledoux, [*The Concentration of Measure Phenomenon*]{} (AMS, Providence, 2001).

S. Boucheron, G. Lugosi, and P. Massart, [*Concentration Inequalities - A Nonasymptotic Theory of Independence*]{} (Oxford University Press, UK, 2013).

M. Talagrand, [Ann. Probab.]{} **24**, 1 (1996).

S. De Bi$\grave{\rm e}$vre and P. E. Parris, [J. Stat. Phys.]{} **168**, 772 (2017).

P. Hayden, D. W. Leung, and A. Winter, [Comm. Math. Phys.]{}, **265**, 95 (2006).

W. Weiss, G. Benenti, G. Casati, I. Guarneri, T. Calarco, M. Paternostro, and S. Montangero, [New J. Phys.]{} **18**, 013021 (2016).

see Supplemental Materials.

Throughout this work we require the sample size be much larger than the transport mean free path $\ell$. The latter is usually defined for a disorder ensemble. As we show in this work, $\ell$ is well defined for a single disordered medium also, and is identical to that defined for an ensemble. In SM (S2) we further show that, when $N$ is large, for a single disordered medium $\ell$ can be also defined by the transparent eigenchannel structure $W_{\tau=1}(x)$.

The group velocity in the air background is unity. In simulations all lengths are in unit of the inverse wave number in the air background.

M. Stoytchev and A. Z. Genack, Phys. Rev. Lett. **79**, 309 (1997).

W. Choi, A. P. Mosk, Q. H. Park, and W. Choi, [ Phys. Rev. B]{} **83**, 134207 (2011).

M. Davy, Z. Shi, J. Park, C. Tian, and A. Z. Genack, [ Nat. Commun.]{} **6**, 6893 (2015).

Th. M. Nieuwenhuizen and M. C. W. van Rossum, Phys. Rev. Lett. **74**, 2674 (1995).

Y. Imry, Europhys. Lett. **1**, 249 (1986).

Here we write $G$ in the representation which is a mixture of the real coordinate $x$ and the waveguide mode. In other words, we perform the Fourier transform of $\langle xy|G|x'y'\rangle$ with respect to $y,y'$.

To make formulae compact we suppress $\omega$ and $\Omega$ in $I_x(c)$ and $E_n$ throughout.

Strictly speaking, the definition of energy density excludes the scalar operator $\hat v_x$. However, we have found that omitting $\hat v_x$ in the definition only leads to the appearance of an irrelevant overall factor.

R. G. Newton, [*Scattering Theory of Waves and Particles*]{} 2nd Ed. (Springer, New York, 1982).

For a slab the width $\gtrsim L$ while for a quasi $1$D medium the width $\ll L$.

O. N. Dorokhov, Zh. Eksp. Teor. Fiz. **85**, 1040 (1983) \[Sov. Phys. JETP **58**, 606 (1983)\].

P. A. Mello, P. Pereyra, and N. Kumar, [Ann. Phys.]{} **181**, 290 (1988).

For the sample-to-sample fluctuations of $I_{x}(c)$ (with the incoming current amplitude $c$ fixed), we are not aware of any reports of either the exponential distribution of diffusive waves (except for the special case of $x=L$ and $c$ chosen to correpond to a plane wave or a Gaussian beam [@Rossum95]) or the stretched-exponential distribution of localized waves.

[**Supplemental Materials**]{}\
by Ping Fang, Liyi Zhao, and Chushun Tian\*

This appendix is written in a self-contained manner. The CM is introduced in great details. The technical details in quantitative (mathematical and numerical) analysis are presented, and extensive discussions on the results are made. The appendix is organized as follows:

-   Sec. S0: It serves as an introduction of CM. We discuss some basics of CM, and introduce the mathematical results used in this work. Particular emphasis are put on their backgrounds and implications.

-   Sec. S1: We present a detailed analysis of the Lipschitz constant $\|I_x\|_{\rm Lip}$, which leads to the general result of Eq. (\[eq:4\]). Then we calculate this general result to obtain Eq. (\[eq:3\]). We further show Eq. (\[eq:14\]).

-   Sec. S2: We study the $2$D and $1$D eigenchannel structures in a slab. Comparisons with eigenchannel structures in quasi $1$D media are made.

-   Sec. S3: We study in details the Lipschitz continuity of $W(x;\omega)$ as a function of the disorder realization $\omega$, which leads to the inequality (\[eq:2\]). Moreover, we derive analytically the expression of $\tilde c(x)$.

-   Sec. S4: We introduce a basic tool of CM, namely, Pisier’s theorem. Armed with this theorem and the results obtained in Sec. S3, we prove the concentration inequality (\[eq:5\]). We further provide numerical evidences showing that it holds for both large and small $N$.

-   Sec. S5: We derive the relation (\[eq:20\]). We perform further numerical studies of this relation.

-   Sec. S6: We discuss more general observables and derive Eq. (\[eq:13\]).

-   Sec. S7: We describe the method of simulations.

-   Sec. S8: We include the derivations of some integrals used in this appendix.

[**S0 Introduction of CM**]{}\
In this section we will illustrate that CM is a phenomenon of high-dimensional geometry, and introduce the background and implications of the mathematical results used in this work. In particular, we will show that the concentration of the uniform measure over the sphere $\mu$ lays a foundation for L$\acute{\rm e}$vy’s lemma and the concentration of Gaussian measure over the Euclidean space $\gamma$ for Pisier’s theorem. We will discuss that these two results are related to each other. By these discussions we hope that it will become clearer to the readers that CM opens “completely new perspectives” [@Boucheron13SM] of probability theory, and a CM-based probability theory differs from traditional probability theories such as the central limit theorem, with which physicists are familiar, in many fundamental aspects.\
\
[**S0.1 Concentration of uniform measure over the sphere**]{}\
In this part we discuss in details a canonical example of CM. Let us take a unit sphere $S^{2N-1}$ in the Euclidean space $\mathds{R}^{2N}$ and normalize its area to unity. Consider a strip $S$ around and symmetric with respect to an equator (Fig. \[fig:S6\]). $S$ has a range of the latitude $\varphi$ between $-\frac{\epsilon}{2}$ and $\frac{\epsilon}{2}$. The north and south poles correspond the latitude of $-\frac{\pi}{2}$ and $\frac{\pi}{2}$, respectively. Using the parametrization (\[eq:S63\])-(\[eq:S65\]) introduced in Sec. S8, we find that the area of the strip is $$\begin{aligned}
\label{eq:S111}
  \mu(S)&=&\frac{\int_{-\frac{\epsilon}{2}}^{\frac{\epsilon}{2}} \sin^{2N-2}\varphi d\varphi}{\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \sin^{2N-2}\varphi d\varphi}\nonumber\\
  &=&\frac{1}{\pi}\frac{(2N-2)!!}{(2N-3)!!}\int_{-\frac{\pi}{2}}^{\frac{\epsilon}{2}} \sin^{2N-2}\varphi d\varphi.\end{aligned}$$

![A strip $S$ around an equator in the sphere (left) and its projection onto a $2$D plane (right). The latitude of the strip $\varphi \in [-\frac{\epsilon}{2}, \frac{\epsilon}{2}]$.[]{data-label="fig:S6"}](typicality_strip.pdf){width="8.7cm"}

From this we find the value of $\epsilon$ corresponding to $\mu(S)=0.9$ so that $S$ covers $90\%$ of the entire area of the sphere. The result is shown in Table \[tab:1\]. From the table we see that as $N$ increases the entire area of the sphere concentrates more and more around the equator. Eventually for large $N$ the entire area almost concentrates on a very narrow strip around the equator.

        $N$          $1$      $2$      $20$    $200$    $2000$
  ---------------- -------- -------- -------- -------- --------
   $\epsilon/\pi$   $0.90$   $0.59$   $0.17$   $0.05$   $0.02$

  : $\epsilon$-values for $\mu(S)=0.9$

\[tab:1\]

As most areas concentrate around the equator, a spherical cap slightly larger than a hemisphere, which has an area of exactly one half, covers almost the entire area of the sphere. This implies that as a spherical cap becomes bigger and bigger, the cap area undergoes a sharp transition once the cap boundary reaches an equator. This geometric phenomenon can be extended to a general set with an area $\geq \frac{1}{2}$, and has a more general and rigorous statement. Specifically, let $A$ be a subset of $S^{2N-1}$. Its $\varepsilon$-[*enlargement*]{}, denoted as $A_\varepsilon$, is defined as a set including all points which are at a geodesic distance $\leq\varepsilon$ from some point in $A$. (In words, any point in $A_\varepsilon$ is close to some point in $A$.) Then we have the following [@Milman86SM]:\
\
Theorem 0.1. [*Let $A$ be a subset in $S^{2N-1}$, whose area $\mu(A)\geq\frac{1}{2}$, and $N$ be an arbitrary natural number. Then*]{} $$\mu(A_\varepsilon)\geq 1-e^{-(N-1)\varepsilon^2}.
    \label{eq:S113}$$ [*Remarks.*]{} (i) According to this theorem, the entire area of the sphere concentrates around a strip around the equator with a width $\sim \frac{1}{\sqrt{N}}$ for large $N$. (ii) The theorem is the formalization of the concentration of the uniform measure $\mu$ over the sphere. (iii) To prove the theorem one needs to use the [*isoperimetric inequality*]{} discovered by L$\acute{\rm e}$vy [@Milman86SM]. We will not give its rigorous statement here. Instead, we explain the inequality in words. That is, in a sphere the circle enclosing a spherical cap is the curve of given length enclosing the largest area.

Next, we would like to probe the concentration of the measure $\mu$. Consider a real-valued function $f$$:$$S^{2N-1}$$\rightarrow$$\mathds{R}$ which has the Lipschitz continuity, with the Lipschitz constant $$\begin{aligned}
\|f\|_{{\rm Lip}}\equiv{\rm sup}_{c,c'}\frac{|f(c)-f(c')|}{\|c-c'\|}<\infty.
\label{eq:S114}\end{aligned}$$ Recall that $\|\cdot\|$ is the Euclidean metric in $\mathds{R}^{2N}$. Furthermore, we define the [*median*]{} or [*L$\acute{e}$vy’s mean*]{} of $f$, denoted as $M[f]$, as the number such that both $\mu(f(c)\leq M[f])\geq \frac{1}{2}$ and $\mu(f(c)\geq M[f])\geq \frac{1}{2}$ are satisfied. Then by using Theorem 0.1 one can readily prove [@Milman86SM] the following:\
\
Corollary 0.2. [*Let $f:S^{2N-1}\rightarrow\mathds{R}$ be Lipschitz. Then for any $\varepsilon>0$, $${\rm Pr}(|f(c)-M[f]|\geq \varepsilon)\leq 2e^{-\frac{(N-1)\varepsilon^2}{\|f\|_{\rm Lip}^2}}.
    \label{eq:S115}$$* ]{}

[*Remarks.*]{} (i) This corollary is similar to L$\acute{\rm e}$vy’s lemma introduced in the paper. It implies the concentration of function $f$ around the median of $M[f]$, and thus is another version of L$\acute{\rm e}$vy’s lemma. We see again that Lipschitz functions serve as nice observables to probe the concentration of measure $\mu$ through the concentration of functions around the constant value $M[f]$. (ii) An important concept follows from the corollary. That is, the sub-Gaussian tail of the distribution of $f$ inherits from Theorem 0.1 or more precisely the isoperimetric inequality, and thus the sub-Gaussian nature of the distribution is of purely geometric origin and not related to the central limit theorem. (iii) $f(c)$ can have a very complicated dependence on $c$. (iv) However, $M[f]$ differs from the mean: $\int fd\mu$ in general. Due to this the corollary is not useful to us in practice.

A natural question is: how can we obtain L$\acute{\rm e}$vy’s lemma used in the paper, which gives the concentration of $f$ around the mean $\int fd\mu$? To answer this question we make the following observation. We parametrize $N$ complex coefficients $c_n$ ($n=1,2,\cdots, N$), which satisfy $\sum_{n=1}^N |c_n|^2=1$ and thus constitute the coordinate of $S^{2N-1}$, as: $$\begin{aligned}
\label{eq:S58}
    c_n=\frac{a'_n+ib'_n}{\sqrt{\sum_{n=1}^N((a'_n)^2+(b'_n)^2)}},\quad a'_n, b'_n\in \mathds{R}.\end{aligned}$$ Then, it is well known [@Milman86SM] that if the $2N$ real variables $a'_n, b'_n$ are independent standard normal random variables, then $c\equiv(c_1,c_2,\cdots,c_N)$ is distributed uniformly over $S^{2N-1}$. This motivates us to investigate the concentration of Gaussian measures $\gamma$ over the Euclidean space of general dimensions. This is the purpose of the next part. In doing so we will sketch how to obtain L$\acute{\rm e}$vy’s lemma used in the paper.\
\
[**S0.2 Concentration of Gaussian measure over the Euclidean space**]{}\
The discussions above are built upon very few notions, namely, the geometric body $S^{2N-1}$ equipped with the geodesic metric and the uniform measure $\mu$. It is possible to extend the discussions to general geometric body equipped with appropriate metric and measure. We refer the general discussions to standard mathematical literatures [@Milman86SM; @Ledoux01SM]. Here we study in details the Euclidean space $\mathds{R}^k$ of general dimension $k$, equipped with the Euclidean metric $\|\cdot\|$ and the standard Gaussian probability measure: $$\label{eq:S116}
    d\gamma\equiv(2\pi)^{-\frac{k}{2}}e^{-\frac{\|\varpi\|^2}{2}}d\varpi,\quad \varpi\in\mathds{R}^k,$$ which we used in this work.

For this measure, there exists a concentration phenomenon similar to what occurs to the uniform measure over the sphere. That is, a set slightly larger than the half space, which has a Gaussian measure of exactly one half, covers most part of the entire measure of the space, despite that the latter is noncompact. With the introduction of $\varepsilon$-enlargement $A_\varepsilon$ (in the same way as before, except that the geodesic distance is replaced by the Euclidean one), this phenomenon can be formalized as the following [@Ledoux01SM]:\
\
Theorem 0.3. [*Let $A$ be a subset in $\mathds{R}^k$, whose Gaussian measure $\gamma(A)\geq\frac{1}{2}$, and $k$ be an arbitrary natural number. Then*]{} $$\gamma(A_\varepsilon)\geq 1-e^{-\varepsilon^2/2}.
    \label{eq:S117}$$ [*Remark.*]{} Compared with Theorem 0.1 we find that this CM phenomenon has a striking feature: it is [*dimension free*]{}, i.e., the right-hand side of the inequality (\[eq:S117\]) has no $k$-dependence.

As before, the concentration of the measure $\gamma$ can be probed by nice observables. To be specific, we take any real-valued function $f$$:$$\mathds{R}^k$$\rightarrow$$\mathds{R}$ which has the Lipschitz continuity, with the Lipschitz constant $$\begin{aligned}
\|f\|_{{\rm Lip}}\equiv{\rm sup}_{\varpi,\varpi'}\frac{|f(\varpi)-f(\varpi')|}{\|\varpi-\varpi'\|}<\infty.
\label{eq:S118}\end{aligned}$$ Then we have the following [@Milman86SM; @Ledoux01SM]:\
\
(Pisier) Theorem 0.4. [*Let $f: \mathds{R}^k\rightarrow \mathds{R}$ be Lipschitz and $\mathds{R}^k$ be equipped with the standard Gaussian measure $\gamma$. Then for any $\varepsilon>0$, $$\label{eq:S119}
    {\rm Pr}\left(\left|f(\varpi)-{\rm E}[f]\right|>\varepsilon\right)\leq 2e^{-\frac{2\varepsilon^2}{(\|f\|_{\rm Lip}\pi)^2}},$$ where the expectation value ${\rm E}[f]=\int fd\gamma$.*]{}\
[*Remarks.*]{} (i) Unlike Corollary 0.2, this theorem refers to the concentration around the mean, rather than the median. (ii) Inheriting from Theorem 0.3 the concentration of $f$ is dimension free also, i.e., the right-hand side of the concentration inequality (\[eq:S119\]) has no $k$-dependence. This also shows that the sub-Gaussian nature of the distribution of $f(\varpi)$ is purely geometric and not related to the central limit theorem. (iii) $f(\varpi)$ can have a very complicated dependence on $\varpi$.

Now consider $k$$=$$2$$N$ and $\varpi$$\equiv$$($$a'_1$,$b'_1$,$\cdots$,$a'_N$,$b'_N)$, where $a'_n$,$b'_n$ are independent Gaussian variables parametrizing $c_n$ in the way shown by Eq. (\[eq:S58\]). Then, by using Theorem 0.4 one can prove L$\acute{\rm e}$vy’s lemma used in the paper [@Milman86SM]:\
\
(L$\acute{\rm e}$vy) Lemma 0.5. [*Let $f: S^{2N-1}\rightarrow \mathds{R}$ be Lipschitz. Then for any $\varepsilon>0$, $$\label{eq:S120}
    {\rm Pr}\left(\left|f(c)-\int f d\mu\right|>\varepsilon\right)\leq 2e^{-
\frac{\delta\varepsilon^2 N}{\|f\|_{{\rm Lip}}^2}},$$ where $\delta$ is some positive absolute constant.*]{}\
The studies of wave-to-wave fluctuations are based on L$\acute{\rm e}$vy’s lemma. As shown above, this lemma is a nontrivial application of Pisier’s theorem. Later on, we will see that Pisier’s theorem is very useful also in the studies of sample-to-sample fluctuations of $W(x;\omega)$ and $L_x(c;\omega)$.\
\
[**S1 Calculations and properties of $\boldsymbol{\|I_x\|_{\rm Lip}}$**]{}\
This section is divided into five subsections. In Sec. S1.1 we show that if $I_x: S^{2N-1}\rightarrow \mathds{R}$ is continuously differentiable, i.e., $$\label{eq:S76}
    {\rm sup}_c\|\nabla_c I_x\|<\infty,$$ where $\nabla_c$ is the orthographic projection of $\nabla$ in $\mathds{R}^{2N}$ onto $S^{2N-1}$, then $I_x(c)$ is Lipschitz. In other words, the continuous differentiability is stronger than the Lipschitz continuity. In doing so, we connect the Lipschitz constant with respect to distinct metrics, namely, the geodesic metric and the Euclidean metric (in $\mathds{R}^{2N}$). In Sec. S1.2 by deriving Eq. (\[eq:4\]) we confirm that the continuous differentiability, namely, the condition (\[eq:S76\]) is indeed satisfied. In Sec. S1.3 we perform numerical analysis and establish the universality of $\|I_x\|_{\rm Lip}$ (with respect to the Euclidean metric) at large $N$. In Sec. S1.4 we explain this universality by CM. In Sec. S1.5 we use this universality to derive Eq. (\[eq:14\]).\
\
[**S1.1 Results of continuous differentiability**]{}\
We start from showing that the continuous differentiability implies the Lipschitz continuity:\
[Lemma 1.1.]{} [*If $I_x(c)$ is continuously differentiable, i.e., the condition (\[eq:S76\]) holds, then*]{} $$|I_x(c)-I_x(c')|\leq \frac{\pi}{2}{\rm sup}_c\|\nabla_c I_x\|\|c-c'\|.
\label{eq:S112}$$

[*Proof.*]{} Suppose that $\gamma (c,c')$ is the shorter geodesic connecting $c,c'\in S^{2N-1}$. Then $$\label{eq:17}
    |I_x(c)-I_x(c')|=\left|\int_{\gamma (c,c')}ds \frac{\partial I_x}{\partial s}\right|\leq \int_{\gamma (c,c')}ds \left|\frac{\partial I_x}{\partial s}\right|,$$ where $ds$ is the differential element of $\gamma (c,c')$. Because of $|\partial I_x/\partial s|\leq \|\nabla_c I_x\|$, we have $$\begin{aligned}
\label{eq:19}
    \int_{\gamma (c,c')}ds \left|\frac{\partial I_x}{\partial s}\right|&\leq& \int_{\gamma (c,c')}ds \|\nabla_c I_x\|\nonumber\\
    &\leq& {\rm sup}_c\|\nabla_c I_x\| \int_{\gamma (c,c')}ds.\end{aligned}$$ Note that $\int_{\gamma (c,c')}ds$ is the geodesic distance of $\gamma$. It satisfies $$\label{eq:27}
\int_{\gamma (c,c')}ds\leq \frac{\pi}{2}\|c-c'\|.$$ Combining the inequalities (\[eq:17\]), (\[eq:19\]) and (\[eq:27\]), we prove the lemma. $\Box$

[*Remarks.*]{} (i) That the continuous differentiability implies the Lipschitz continuity is actually very general. In particular, one may replace $S^{2N-1}$ by arbitrary space $\mathscr{C}$ (equipped with certain metric). (ii) The lemma shows that $$\label{eq:S77}
    \|I_x\|_{{\rm Lip}}\leq \frac{\pi}{2}{\rm sup}_c\|\nabla_c I_x\|.$$ On the other hand, because of $$\begin{aligned}
    \lim_{c'\rightarrow c}\frac{|I_x(c)-I_x(c')|}{\|c-c'\|}\leq \|\nabla_c I_x\|\leq {\rm sup}_{c}\|\nabla_c I_x\|,
    \label{eq:S78}\end{aligned}$$ we have $$\begin{aligned}
    \|I_x\|_{{\rm Lip}}\geq {\rm sup}_{c}\|\nabla_c I_x\|.
    \label{eq:S101}\end{aligned}$$ The inequalities (\[eq:S77\]) and (\[eq:S101\]) imply that $$\label{eq:S80}
    \|I_x\|_{{\rm Lip}}=a{\rm sup}_c\|\nabla_c I_x\|,\quad 1\leq a\leq \pi/2,$$ where $a$ is some numerical constant. Because $a$ does not grow with $N$ and is close to unity, it does not affect any results obtained from L$\acute{\rm e}$vy’s lemma, except slightly modifying the absolute constant $\delta$ in the sub-Gaussian tail. Thus we do not pay attention to its exact value any more and set $a=\pi/2$, i.e., $$\label{eq:S51}
    \|I_x\|_{{\rm Lip}}=\frac{\pi}{2}{\rm sup}_c\|\nabla_c I_x\|\equiv {\rm sup}_c L_x(c;\omega).$$ This provides an analytic formula of $\|I_x\|_{{\rm Lip}}$, which we will use in subsequent mathematical and numerical analysis.\
\
[**S1.2 Derivations of Eq. (\[eq:4\]): confirming continuous differentiability**]{}\
In this part we use Eq. (\[eq:S51\]) to derive Eq. (\[eq:4\]). Let $c_n=a_n+ib_n$, where $a_n,b_n\in \mathds{R}$. By using the definition (\[eq:10\]) of $I_x(c)$, we obtain $$\begin{aligned}
\label{eq:S54}
    \|\nabla I_x\|^2&=&\sum_{n=1}^N \left(\left(\frac{\partial I_x}{\partial a_n}\right)^2+\left(\frac{\partial I_x}{\partial b_n}\right)^2\right)\nonumber\\
    &=&4\sum_{n=1}^N \frac{\partial I_x}{\partial c_n}\frac{\partial I_x}{\partial c_n^*},\end{aligned}$$ where $$\begin{aligned}
\label{eq:30}
    \frac{\partial I_x}{\partial c_n}&=&\sum_{n'=1}^N \!c_{n'}^* E_{n'}^\dagger(x)\cdot E_{n}(x),\nonumber\\
    \frac{\partial I_x}{\partial c_n^*}&=&\sum_{n'=1}^N \!c_{n'} E_{n}^\dagger(x)\cdot E_{n'}(x).\end{aligned}$$ Denote $\boldsymbol{n}\equiv (a_1,b_1,\cdots,a_N,b_N)$. It is easy to see that $\boldsymbol{n}$ is the radial direction of $S^{2N-1}$. Moreover, because of $\sum_{n=1}^N(a_n^2+b_n^2)=1$ it is a unit vector. As a result, $$\begin{aligned}
\label{eq:S55}
    \boldsymbol{n}\cdot\nabla I_x&=&\sum_{n=1}^N \left(a_n\frac{\partial I_x}{\partial a_n}+b_n\frac{\partial I_x}{\partial b_n}\right)\nonumber\\
    &=&\sum_{n=1}^N \left(c_n\frac{\partial I_x}{\partial c_n}+c_n^*\frac{\partial I_x}{\partial c_n^*}\right).\end{aligned}$$ By the definition (\[eq:10\]), we find $$\begin{aligned}
\label{eq:S56}
    \sum_{n=1}^Nc_n\frac{\partial I_x}{\partial c_n}=\sum_{n=1}^Nc_n^*\frac{\partial I_x}{\partial c_n^*}=I_x.\end{aligned}$$ Substituting it into Eq. (\[eq:S55\]) gives $$\begin{aligned}
\label{eq:S57}
    \boldsymbol{n}\cdot\nabla I_x=2I_x.\end{aligned}$$ Combining Eqs. (\[eq:S54\]) and (\[eq:S57\]) we find $$\begin{aligned}
\label{eq:29}
    L_x(c;\omega)&=&\frac{\pi}{2}\|\nabla_c I_x\|\nonumber\\
    &=&\frac{\pi}{2}\sqrt{\|\nabla I_x\|^2-|\boldsymbol{n}\cdot\nabla I_x|^2}\nonumber\\
    &=&\pi\sqrt{\sum_{n=1}^N\frac{\partial I_x}{\partial c_n}\frac{\partial I_x}{\partial c_n^*}-I_x^2}.\end{aligned}$$ Equations (\[eq:S51\]), (\[eq:30\]) and (\[eq:29\]) give Eq. (\[eq:4\]). According to Eq. (\[eq:4\]), the condition (\[eq:S76\]) is satisfied, i.e., $I_x$ is continuously differentiable.\
\
[**S1.3 Results of Eq. (\[eq:4\])**]{}\
It is difficult to calculate Eq. (\[eq:4\]) explicitly and analytically, especially in view of that $\omega$ is fixed. However, to calculate Eq. (\[eq:4\]) numerically is relatively easy. The point is that, owing to the analytic expression of $L_x(c;\omega)$, we do not need to perform numerical differentiation which is technically demanding. To be specific, given $\omega$ we can simulate Eq. (\[eq:15\]) by using the method described in Sec. S7 and obtain the set of $\{E_n(x)\}$. Upon substituting the result into Eq. (\[eq:4\]) we obtain $L_x(c;\omega)$ and $\|I_x\|_{\rm Lip}$. In this part we put this numerical scheme into practice and calculate Eq. (\[eq:4\]) explicitly. Surprisingly, we will show that for large $N$, $\|I_x\|_{\rm Lip}$ displays universalities, and even stronger results hold for $L_x(c;\omega)$.

Given $N$ we draw $c$ randomly from a uniform distribution over $S^{2N-1}$ by using the algorithm described in Sec. S7. Then we simulate the wave propagation in a fixed disorder realization $\omega$ (the details are described in Sec. S7.), and calculate $\|\nabla_c I_x\|$ or equivalently $L_x(c;\omega)$ by using the method described above. This gives a depth profile of $\|\nabla_c I_x\|$. We repeat the numerical experiment for $10^4$ randomly chosen $c$, and thus have $
10^4$ profiles in total. We have checked that the number of $
10^4$ is large enough to ensure that numerical results converge. Finally, the procedures are repeated for several distinct $N$.

To obtain the depth profile of $\|I_x\|_{\rm Lip}$, we perform the numerical analysis of the profiles of $\|\nabla_c I_x\|$ in the following way. (i) We average the profiles corresponding to the same $N$. The average profile, $\int \|\nabla_c I_x\|d\mu
$, is shown in the main panel of Fig. \[fig:S1\] for three different $N$. We see that these three profiles collapse into a single curve which is continuous in $x$, and the value of $\int \|\nabla_c I_x\|d\mu$ is order of unity at every $x$. (ii) We calculate the maximal deviation $\delta_I(x)$ from the average value $\int \|\nabla_c I_x\|d\mu$ for distinct $x$. As shown in the inset of Fig. \[fig:S1\], $\delta _I(x)$ scales with $N$ as $$\label{eq:S49}
    \delta _I(x)\sim \frac{1}{\sqrt{N}},\quad \forall\, x\in [0,L].$$ This implies that for large $N$, the two profiles, $\|\nabla_c I_x\|$ and $\int \|\nabla_c I_x\|d\mu$, converge, giving $$\label{eq:S50}
    \|\nabla_c I_x\|=\int \|\nabla_c I_x\|d\mu,$$ i.e., $\|\nabla_c I_x\|$ is independent of $c$.

By using Eq. (\[eq:4\]) namely Eq. (\[eq:S51\]) and Eq. (\[eq:S50\]) we find $$\label{eq:S1}
    \|I_x\|_{{\rm Lip}}=\frac{\pi}{2} \int \|\nabla_c I_x\|d\mu.$$ Since as shown in (i) the right-hand side has no $N$-dependence and is order of unity, Eq. (\[eq:3\]) follows. The curves in Fig. \[fig:S1\], when multiplied by $\pi/2$, give the three $\|I_x\|_{\rm Lip}$-curves in Fig. \[fig:4\].

The finding of Eq. (\[eq:S50\]) motivates us to perform a careful study of the universalities of $\|\nabla_c I_x\|=\frac{2}{\pi}L_x(c;\omega)$. We perform simulations for large $N(=800)$. Surprisingly, we find that $L_x(c;\omega)$ is independent not only of $c$, but also of $\omega$, as shown in Fig. \[fig:4\]. In other words, $L_x(c;\omega)$ gives a universal curve parametrized by $x$ for very large $N$.\
\

![Main panel: Simulations show that the profiles of $\int \|\nabla_c I_x\|d\mu$ corresponding to $N=160,240,320$ collapse into a single curve, and $\int \|\nabla_c I_x\|d\mu={\cal O}(1)$ at every $x$. Inset: Simulations show that given $x$ the maximal deviation $\delta_I(x)$ from $\int \|\nabla_c I_x\|d\mu$ decreases with $N$ as $N^{-1/2}$.[]{data-label="fig:S1"}](Lipschitz_constant.pdf){width="8.7cm"}

[**S1.4 Universality of $\boldsymbol{L_x(c;\omega)}$ from CM**]{}\
In this part we employ CM to explain the universality of $L_x(c;\omega)$ with respect to $c$ and $\omega$. We stress that this is not a rigorous mathematical proof, which is far beyond the present work.

First of all, we fix $\omega$. From Eq. (\[eq:4\]) or Eq. (\[eq:29\]) we see that $L_x(c;\omega)$ is continuously differentiable, i.e., $$\label{eq:S106}
    {\rm sup}_c\|\nabla_c L_x\|<\infty,$$ and thus is Lipschitz. Let the Lipschitz constant of $L_x(c;\omega)$ ($\omega$ fixed) be $\|L_x(\omega)\|_{\rm Lip}$. By making use of Lemma 0.5 we obtain $$\begin{aligned}
\label{eq:S107}
  {\rm Pr}\left(\left|L_x(c;\omega)-\int L_x(c;\omega)d\mu\right|>\varepsilon\right)\leq 2e^{-
        \frac{\delta\varepsilon^2 N}{\|L_x(\omega)\|_{{\rm Lip}}^2}}.\nonumber\\\end{aligned}$$ From this we see that provided $$\label{eq:S108}
    \int L_x(c;\omega)d\mu\gg \|L_x(\omega)\|_{{\rm Lip}}/\sqrt{N},$$ strong concentration around the mean: $\int L_xd\mu$ results. This gives the universality with respect to $c$.

Next, we fix $c$. From Eq. (\[eq:4\]) or Eq. (\[eq:29\]) we may expect that $L_x(c;\omega)$ is continuously differentiable with respect to $\omega\equiv\{\omega_{(x,y)}\}$, when the medium space is discretized into $M$ lattice points. Thus we have $$\begin{aligned}
\label{eq:S109}
  |L_x(c;\omega)-L_x(c;\omega')|\leq
  \|L_x(c)\|_{{\rm Lip}}\, \|\omega-\omega'\|,\end{aligned}$$ where $\|L_x(c)\|_{{\rm Lip}}$ is the corresponding Lipschitz constant. The ensemble of disorder realizations $\omega$, by definition, is described by a Gaussian probability measure on $\mathds{R}^M$ of zero mean and variance $\sigma^2$. Let $\varpi$ in (Pisier’s) Theorem 0.4. be $\omega$, $f(\varpi)$ be $L_x(c;\omega)$, and $k=M$. After rescaling we obtain from this theorem $$\label{eq:S110}
    {\rm Pr}\left(|L_x(c;\omega)-{\rm E}[L_x(c;\omega)]|>\varepsilon\right)\leq 2e^{-\frac{2\varepsilon^2}{(\pi\|L_x(c)\|_{{\rm Lip}}\sigma)^2}}.$$ It is important to note that for this concentration inequality the exponent of the sub-Gaussian tail bound has no explicit $N$-dependence. Now, if we assume that $\|L_x(c)\|_{{\rm Lip}}$ decays with $N$, then for sufficiently large $N$ we have strong concentration of $L_x(c;\omega)$ around the disorder average: ${\rm E}[L_x(c;\omega)]$. This explains the universality of $L_x(c;\omega)$ with respect to $\omega$.

Summarizing, the non-rigorous discussions above indicate that the universalities of $L_x(c;\omega)$ may have deep connections to CM.\
\
[**S1.5 Derivations of Eq. (\[eq:14\])**]{}\
In this part we study the consequences of the universalities of $L_x(c;\omega)$. By Jensen’s inequality we have: $$\begin{aligned}
\label{eq:S105}
  \|I_x\|_{{\rm Lip}}^2=\left(\int L_x(c;\omega)d\mu\right)^2\leq \int L_x^2(c;\omega)d\mu\end{aligned}$$ in general. With the universality of $L_x(c;\omega)$ with respect to $c$ taken into account, a stronger result follows, i.e., $$\label{eq:S88}
    \|I_x\|_{{\rm Lip}}^2=\int L_x^2(c;\omega)d\mu.$$ Now we show that this gives Eq. (\[eq:14\]).

Substituting Eq. (\[eq:4\]) into it we obtain $$\begin{aligned}
\label{eq:S89}
\frac{\|I_x\|_{{\rm Lip}}^2}{\pi^2} = \sum_{n,n',m'}^N\!\! \left(\int c_{m'}^*c_{n'} d\mu\right) E_{m'}^\dagger\cdot E_{n}E_{n}^\dagger\cdot E_{n'}\quad\nonumber\\
  -\!\sum_{n,n',m,m'}^N\!\! \left(\int c_{n}^*c_{n'}c_{m}^*c_{m'} d\mu\right) E_{n}^\dagger\cdot E_{n'}E_{m}^\dagger\cdot E_{m'}.\,\,\end{aligned}$$ According to the integral formula (\[eq:S79\]) derived in Sec. S8, for the integral $\int c_{m'}^*c_{n'} d\mu$ not to vanish it is necessary $m'=n'$. For the integral $\int c_{n}^*c_{n'}c_{m}^*c_{m'} d\mu$ not to vanish, it is necessary that one of the following conditions: (i) $n=n', m=m' (n\neq m)$; (ii) $n=m', m=n' (n\neq m)$; (iii) $n=n'=m=m'$ must be satisfied. As a result, $$\begin{aligned}
\label{eq:S90}
    &&\int c_{n}^*c_{n'}c_{m}^*c_{m'} d\mu\nonumber\\
    &=&\int |c_{n}|^2|c_{m}|^2 d\mu(\delta_{nn'}\delta_{mm'}+\delta_{nm'}\delta_{mn'})(1-\delta_{nm})\nonumber\\
    &+&\int |c_{n}|^4 d\mu \delta_{nn'}\delta_{mm'}\delta_{nm}.\end{aligned}$$ Using the integral formulae (\[eq:S60\]) and (\[eq:S61\]) derived in Sec. S8, we obtain $$\begin{aligned}
\label{eq:S91}
    &&\int c_{n}^*c_{n'}c_{m}^*c_{m'} d\mu\nonumber\\
    &=&\frac{1}{N(N+1)}(\delta_{nn'}\delta_{mm'}+\delta_{nm'}\delta_{mn'})(1-\delta_{nm})\nonumber\\
    &+&\frac{2}{N(N+1)}\delta_{nn'}\delta_{mm'}\delta_{nm}.\end{aligned}$$ Using this result we reduce Eq. (\[eq:S89\]) to $$\begin{aligned}
\label{eq:S92}
\frac{\|I_x\|_{{\rm Lip}}^2}{\pi^2} = \frac{1}{N}\sum_{n,m=1}^N |E_{n}^\dagger\cdot E_{m}|^2\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\nonumber\\
  -\frac{1}{N(N+1)}\sum_{n\neq m}^N \left(E_{n}^\dagger\cdot E_{n}E_{m}^\dagger\cdot E_{m}+|E_{n}^\dagger\cdot E_{m}|^2\right)\nonumber\\
  -\frac{2}{N(N+1)}\sum_{n=1}^N |E_{n}^\dagger\cdot E_{n}|^2.\quad\quad\quad\quad\quad\quad\quad\quad\quad\,\,\,\,\end{aligned}$$ From this we obtain $$\begin{aligned}
\label{eq:S121}
  \|I_x\|_{\rm Lip}^2&=&\frac{\pi^2}{N+1}\bigg(\sum_{n=1}^N\left(W_{\tau_n}(x)-W(x;\omega)\right)^2\nonumber\\
  &&+\sum_{n\neq n'}^N|E_{n}^\dagger(x)\cdot E_{n'}(x)|^2\bigg).\end{aligned}$$ after simple algebra. Recall that $E_n(x)\equiv \{E_{na}(x)\}$ and $a$ labels the ideal waveguide modes. We find that Eq. (\[eq:S121\]) is equivalent to Eq. (\[eq:14\]).\
\

![Top: Simulations show that the $2$D structure of transparent eigenchannel ($\tau=1$) in a slab (left, the width-to-length ratio $=48$) is very different from that in quasi $1$D (upper right, the ratio $=0.1$): the former exhibits localization in the transverse ($y$) direction while the latter not. Surprisingly, despite this difference, upon integrating $y$ we find that the $1$D structure $W_\tau(x)$ of a single slab – without ensemble averaging – is well described by Eq. (\[eq:S3\]) for the ensemble of quasi $1$D disorder media (lower right). Bottom: the same as top, but $\tau=0.5$.[]{data-label="fig:S3"}](2Dprofile.pdf){width="8.7cm"}

[**S2 Eigenchannel structures in a slab**]{}\
As the impacts of high dimension on eigenchannel structures $\{W_{\tau_n}(x)\}$ remain largely unexplored, in this section we study the explicit forms of $W_{\tau_n}(x)$ in a single disordered slab and their dependence on $N,\omega$. Specifically, we will use the numerical method described in Sec. S7 to simulate the profiles, $\{|E_{n}(x,y)|^2\}$ and $\{W_{\tau_n}(x)\}$, respectively, and compare the results for slabs with those for quasi $1$D media.

The simulation results are shown in Fig. \[fig:S3\]. We see that an eigenchannel in a slab and the one in a quasi $1$D medium, even they correspond to the same $\tau$, have totally different $2$D structures: in a slab localization structures appear in the transverse ($y$) direction while in quasi $1$D not. This notwithstanding, surprisingly, when integrate out $y$ we find that the ensuing $1$D structure $W_{\tau_n}(x)$ – despite it is for a single slab – is well described by a universal analytic formula $W_\tau(x)$ originally derived for an ensemble of quasi $1$D disordered media [@Tian15SM]. For the convenience below we include the formula here. Let the transmission eigenvalue be parametrized as $$\label{eq:S2}
    \tau=\frac{1}{\cosh^2\phi},\quad \phi\geq 0.$$ Then the analytic expression of $W_\tau(x)$ is as follows: $$\begin{aligned}
  W_\tau(x) = S_\tau(x)W_{\tau=1}(x),
  \label{eq:S103}\end{aligned}$$ where $$\begin{aligned}
  W_{\tau=1}(x) = 1+\frac{\pi x(L-x)}{2L\ell}
  \label{eq:S104}\end{aligned}$$ (recall that $\ell$ is the transport mean free path.) and $$\begin{aligned}
  S_\tau(x) = 2\frac{\cosh^2(h(x')(1-x')\phi)}{\cosh^2(h(x')\phi)}-\tau
  \label{eq:S3}\end{aligned}$$ with $x'=x/L$. Here $h(x')$ is a function that increases monotonically from $h(1)=1$ as $x'$ decreases from $1$. (For diffusive waves) both $W_{\tau=1}(x)$ and $h(x')$ are independent of $N$.

Because in a single disordered slab the transparent eigenchannel structure is in agreement with Eq. (\[eq:S104\]). When fitting the former by the latter, we obtain $\ell$ numerically for a single disorder realization.\
\
[**S3 Lipschitz continuity of $\boldsymbol{W(x;\omega)}$**]{}\
In this section we will study the Lipschitz continuity of $W(x;\omega)$ (as a function of $\omega$). Specifically, we will show the inequality (\[eq:2\]) and derive an analytic expression of $\tilde c(x)$. Below to keep the proof simple we will use the result shown numerically in Sec. S2, i.e., that the analytic expression of $W_\tau(x)$ derived originally for quasi $1$D disordered media [@Tian15SM] applies also for a single disordered slab. However, as we will discuss in the end of this section, the key result actually does not depend on the details of the expression, but on its general properties regarding the continuity.

The proof includes following three steps.\
\
[**S3.1 Step I**]{}\
We introduce a family of real-valued functions (labeled by $x$) defined as follows: $$\label{eq:S5}
    g_{z}(x)\equiv W_\tau(x)=W_{z^2}(x),\quad z\equiv\sqrt{\tau}.$$ In this step we will study the Lipschitz continuity of this family of functions.

As discussed in Sec. S1.1, if a function is continuously differentiable, then it is Lipschitz. Thus we consider the derivative of $g_{z}$, which is $$\label{eq:S6}
    \frac{d}{dz}g_{z}(x)=\frac{d}{dz}W_{z^2}(x).$$ Substituting Eqs. (\[eq:S103\])-(\[eq:S3\]) into Eq. (\[eq:S6\]), we obtain $$\begin{aligned}
\label{eq:S7}
  \frac{d}{dz}g_{z}(x)=2F(\phi)W_{\tau=1}(x),\end{aligned}$$ where $$\begin{aligned}
\label{eq:S8}
  F(\phi)=\frac{\cosh^2\phi}{\sinh\phi}
  \frac{\cosh((1-x')h(x')\phi)}{\cosh^3(h(x')\phi)}h(x')\quad\quad\nonumber\\
  \times (x'\sinh((2-x')h(x')\phi)+(2-x')\sinh(x'h(x')\phi))\nonumber\\
  -\frac{1}{\cosh\phi}.\quad\quad\quad\quad\quad\end{aligned}$$ For fixed $x'\in [0,1]$ the function $F(\phi)$ is continuous and finite at $\phi=0$. On the other hand, it is well known [@Dorokhov84SM; @Mello88SM] that $\phi$ has an upper bound $\phi_m$ which is independent of $N$. Thus the maximal value of $|F(\phi)|$ is independent of $N$. Taking this into account, we find from Eq. (\[eq:S7\]) that the derivative of $g_{z}(x)$, with $$\label{eq:S24}
    z\in \mathscr{C}\equiv[1/\cosh\phi_m,1],$$ has a least upper bound which is independent of $N$. This bound gives the Lipschitz constant, $$\label{eq:S82}
    \|g_x\|_{\rm Lip}={\rm sup}_{z\in \mathscr{C}}\left|\frac{d}{dz}g_z(x)\right|,$$ of the function $g_{z}(x)$. As a result, $$\label{eq:S9}
    |g_{z}(x)-g_{z'}(x)|\leq \|g_x\|_{\rm Lip}|z-z'|,\, \forall\, z,z'\in \mathscr{C}.$$ We thus achieve the result of the first step.\
\
[**S3.2 Step II**]{}\
To proceed we introduce the Hilbert-Schmidt (HS) norm of a matrix $A=\{A_{ij}\in \mathds{C}\}$, defined as $$\label{eq:S10}
    \|A\|_{\rm HS}\equiv \sqrt{\sum_{ij}|A_{ij}|^2}=\sqrt{{\rm Tr}(A^\dagger A)}.$$ In this step we will bound $|W(x;\omega)-W(x;\omega')|$ by $\|t(\omega)-t(\omega')\|_{\rm HS}$.

By using the Cauchy-Schwartz inequality, we have $$\begin{aligned}
\label{eq:S13}
    &&|W(x;\omega)-W(x;\omega')|\nonumber\\
    &=&\frac{1}{N}\left|\sum_{n=1}^N\left(W_{\tau_n(\omega)}(x)-W_{\tau_n(\omega')}(x)\right)\right|\nonumber\\
    &\leq&\frac{1}{\sqrt{N}}\left(\sum_{n=1}^N\left(W_{\tau_n(\omega)}(x)-W_{\tau_n(\omega')}(x)\right)^2\right)^{1/2}.\quad\end{aligned}$$ With the HS norm we can rewrite the inequality (\[eq:S13\]) as $$\begin{aligned}
\label{eq:S14}
    &&|W(x;\omega)-W(x;\omega')|\nonumber\\
    &\leq&\frac{1}{\sqrt{N}}\left\|W_{tt^\dagger(\omega)}(x)-W_{tt^\dagger(\omega')}(x)\right\|_{\rm HS}.\end{aligned}$$ Throughout this section we put explicitly the dependence of the $N\times N$ matrix $t$ and the eigenvalue spectrum $\{\tau_n\}$ on the disorder realization $\omega$.

To proceed we show the following\
\
[Lemma 3.1.]{} [*Let $\mathscr{C}$ be a subset of $\mathds{R}$ and $f:\mathscr{C}\rightarrow \mathds{R}$ be Lipschitz, i.e., $$\label{eq:S11}
    |f(z)-f(z')|\leq \|f\|_{{\rm Lip}}\,|z-z'|,\quad \forall\, z,z'\in \mathscr{C}.$$ Then for any pair $A,A'$ of $N\times N$ hermitian matrices whose eigenvalues are in $\mathscr{C}$, we have*]{} $$\label{eq:S12}
    \|f(A)-f(A')\|_{\rm HS}\leq \|f\|_{{\rm Lip}}\,\|A-A'\|_{\rm HS}.$$

[*Proof.*]{} We follow the scheme described in Ref. [@Haaergrup03SM]. Because $A$ is hermitian we can write it as $$\label{eq:S15}
    A=\sum_{n=1}^N \lambda_n {\cal P}_n,\quad {\cal P}_n\equiv |n\rangle\langle n|,$$ where $\{{\cal P}_n\}$ are $1$D mutually orthogonal projectors, and $|n\rangle$ is the eigenvector corresponding to the eigenvalue $\lambda_n\in \mathscr{C}$. For $A'$ we have the eigenvalue equation, $$\label{eq:S16}
    A'|\psi_n\rangle=\mu_n |\psi_n\rangle.$$ In the bases of $\{|n\rangle\}$, the eigenvector $|\psi_n\rangle$ corresponding to $\mu_n\in \mathscr{C}$ can be written as $$\label{eq:S17}
    |\psi_n\rangle=\sum_{n'=1}^N C_{nn'}|n'\rangle, \quad \sum_{n'=1}^N |C_{nn'}|^2=1.$$ From this we have the following expression for $A'$, $$\label{eq:S18}
    A'=\sum_{m=1}^N \mu_m {\cal P}'_m,\, {\cal P}'_m\equiv \sum_{n,n'=1}^N C_{mn}C_{mn'}^* |n\rangle\langle n'|.$$

By using Eqs. (\[eq:S15\]) and (\[eq:S18\]), we obtain $$\begin{aligned}
\label{eq:S19}
  && \|f(A)-f(A')\|_{\rm HS}^2\nonumber\\
  &=& {\rm Tr}\bigg(\sum_{m=1}^N \left(f^2(\lambda_m){\cal P}_m+f^2(\mu_m){\cal P}_m'\right)\nonumber\\
  &&-2\sum_{m,m'=1}^N f(\lambda_m)f(\mu_{m'}){\cal P}_{m}{\cal P}_{m'}'\bigg).\end{aligned}$$ With the insertion of the identity: $$\label{eq:S81}
    \sum_{m=1}^N{\cal P}_{m}=\sum_{m=1}^N{\cal P}_{m}'=1,$$ Eq. (\[eq:S19\]) is reduced to $$\begin{aligned}
\label{eq:S20}
  && \|f(A)-f(A')\|_{\rm HS}^2\nonumber\\
  &=& \sum_{m,m'=1}^N \left(f(\lambda_m)-f(\mu_m)\right)^2{\rm Tr}\left({\cal P}_{m}{\cal P}_{m'}'\right).\end{aligned}$$ On the other hand, $$\begin{aligned}
\label{eq:S21}
  {\rm Tr}\left({\cal P}_{m}{\cal P}_{m'}'\right)
  &=& {\rm Tr} \sum_{n,n'=1}^N \left(C_{m'n}C_{m'n'}^* |m\rangle\langle m |n\rangle\langle n'|\right)\nonumber\\
  &=& |C_{m'm}|^2\nonumber\\
  &\geq& 0.\end{aligned}$$ Thus all the summands in Eq. (\[eq:S20\]) are nonnegative. With this result and the condition: $\lambda_m,\mu_m\in \mathscr{C}$, we can apply the inequality (\[eq:S11\]) to Eq. (\[eq:S20\]) and obtain $$\begin{aligned}
\label{eq:S22}
  &&\|f(A)-f(A')\|_{\rm HS}^2\nonumber\\
  &\leq& \|f\|_{\rm Lip}^2
  \sum_{m,m'=1}^N \left(\lambda_m-\mu_m\right)^2{\rm Tr}\left({\cal P}_{m}{\cal P}_{m'}'\right).\end{aligned}$$ For the sum we can apply the identity (\[eq:S20\]) with $f$ being the identity map. As a result, $$\begin{aligned}
\label{eq:S23}
  \sum_{m,m'=1}^N \left(\lambda_m-\mu_m\right)^2{\rm Tr}\left({\cal P}_{m}{\cal P}_{m'}'\right)=\|A-A'\|_{\rm HS}^2.\end{aligned}$$ Combining it with the inequality (\[eq:S22\]) we complete the proof of lemma. $\Box$

[*Remark.*]{} The lemma differs from a well-known result [@Haaergrup03SM] in the conditions. That is, we require $\mathscr{C}$ to be a subset of $\mathds{R}$, instead of $\mathds{R}$, and correspondingly, the eigenvalues of $A,A'$ to be in this subset. This modification is important for subsequent analysis.

Let $\mathscr{C}$ in Lemma 3.1 be $[1/\cosh\phi_m,1]$ and $f(z)$ be $g_z(x)$. Consider the hermitian matrices, (the derivations below follow Ref. [@Haaergrup03SM].) $$\begin{aligned}
\label{eq:S25}
  \tilde t(\omega)=\left(
                     \begin{array}{cc}
                       0 & t^\dagger(\omega) \\
                       t(\omega) & 0 \\
                     \end{array}
                   \right),\,
  \tilde t(\omega')=\left(
                     \begin{array}{cc}
                       0 & t^\dagger(\omega') \\
                       t(\omega') & 0 \\
                     \end{array}
                   \right).\,\,\end{aligned}$$ It is easy to see that their eigenvalues are $\{\sqrt{\tau_n(\omega)}\}$ and $\{\sqrt{\tau_n(\omega')}\}$, respectively, and thus are in the set of $\mathscr{C}$. Using the inequality (\[eq:S9\]) and Lemma 3.1 we find that $$\label{eq:S26}
    \|g_{\tilde t(\omega)}(x)-g_{\tilde t(\omega')}(x)\|_{\rm HS}
    \leq \|g_x\|_{\rm Lip}\,\|\tilde t(\omega)-\tilde t(\omega')\|_{\rm HS}.$$ By definition (\[eq:S5\]) we have $$\begin{aligned}
\label{eq:S27}
  &&g_{\tilde t(\omega)}-g_{\tilde t(\omega')} \nonumber\\
  &=& W_{\tilde t^2(\omega)}-W_{\tilde t^2(\omega')} \nonumber\\
  &=& \left(
        \begin{array}{cc}
          W_{t^\dagger t(\omega)}-W_{t^\dagger t(\omega')} & 0 \\
          0 & W_{tt^\dagger(\omega)}-W_{tt^\dagger(\omega')} \\
        \end{array}
      \right),\,\,\,\end{aligned}$$ where to make the formula compact we have suppressed the parameter $x$. Substituting it into the inequality (\[eq:S26\]) we obtain $$\begin{aligned}
\label{eq:S28}
  \|W_{t^\dagger t(\omega)}-W_{t^\dagger t(\omega')}\|_{\rm HS}^2+\|W_{tt^\dagger(\omega)}-W_{tt^\dagger(\omega')}\|_{\rm HS}^2\quad\nonumber\\
  \leq \|g_x\|_{\rm Lip}^2\left(\|t(\omega)-t(\omega')\|_{\rm HS}^2+\|t^\dagger(\omega)-t^\dagger(\omega')\|_{\rm HS}^2\right),\,\,\end{aligned}$$ from which $$\begin{aligned}
\label{eq:S29}
  &&\|W_{tt^\dagger(\omega)}(x)-W_{tt^\dagger(\omega')}(x)\|_{\rm HS}\nonumber\\
  &\leq& \|g_x\|_{\rm Lip}\|t(\omega)-t(\omega')\|_{\rm HS}\end{aligned}$$ follows.

Combining the inequalities (\[eq:S14\]) and (\[eq:S29\]), we obtain $$\begin{aligned}
\label{eq:S30}
|W(x;\omega)-W(x;\omega')|
\leq\frac{\|g_x\|_{\rm Lip}}{\sqrt{N}}\left\|t(\omega)-t(\omega')\right\|_{\rm HS}.\end{aligned}$$ We thus achieve the result of the second step.\
\
[**S3.3 Step III**]{}\
In the last step we will further bound $\left\|t(\omega)-t(\omega')\right\|_{\rm HS}$ on the right-hand side of the inequality (\[eq:S30\]) by $\|\omega-\omega'\|$. To this end we note that the matrix element of $t(\omega)$ is given by $$\begin{aligned}
\label{eq:S31}
  t_{ab}(\omega)
  &=& -i\sqrt{\tilde v_a\tilde v_b} \int\!\!\!\!\int dydy' \varphi_a(y)\varphi_b^*(y')\nonumber\\
  &\times& G(x=\infty,y,x'=-\infty,y')\nonumber\\
  &\equiv&-i\sqrt{\tilde v_a\tilde v_b}\langle \infty a|G(\omega)|-\infty b\rangle,\end{aligned}$$ where in the last line we have put the $\omega$-dependence of Green’s function $G$ explicitly. With the substitution of this expression, $\left\|t(\omega)-t(\omega')\right\|_{\rm HS}$ is written as $$\begin{aligned}
\label{eq:S32}
  && \left\|t(\omega)-t(\omega')\right\|_{\rm HS}^2\nonumber\\
  &=& \sum_{a,b=1}^N \tilde v_a\tilde v_b |\langle \infty a|G(\omega)-G(\omega')|-\infty b\rangle|^2.\end{aligned}$$ With the help of the identity, $$\begin{aligned}
\label{eq:S33}
G(\omega)-G(\omega')
&=&G(\omega)\left(G^{-1}(\omega')-G^{-1}(\omega)\right)G(\omega')\quad\nonumber\\
&=&G(\omega)\left(\omega-\omega'\right)G(\omega'),\end{aligned}$$ we write Eq. (\[eq:S32\]) as $$\begin{aligned}
\label{eq:S34}
  && \left\|t(\omega)-t(\omega')\right\|_{\rm HS}^2\nonumber\\
  &=& \sum_{a,b=1}^N \tilde v_a\tilde v_b |\langle \infty a|G(\omega)\left(\omega-\omega'\right)G(\omega')|-\infty b\rangle|^2.\quad\quad\end{aligned}$$ After simple algebra we further reduce it to $$\label{eq:S35}
    \left\|t(\omega)-t(\omega')\right\|_{\rm HS}^2=(\omega-\omega')^{\rm T}{\cal M}(\omega,\omega')(\omega-\omega').$$ Here ${\cal M}(\omega,\omega')$ is an $M\times M$ (recall that $M$ is the number of lattice points when the disordered medium is discretized.) hermitian matrix. The corresponding matrix element ${\cal M}_{rr'}(\omega,\omega')$, with $r\equiv (x,y)$ and $r'\equiv (x',y')$ being the discrete lattice points, reads $$\begin{aligned}
\label{eq:S36}
    &&{\cal M}_{rr'}(\omega,\omega')\nonumber\\
    &=&\left(G^\dagger(\omega){\cal V}_LG(\omega)\right)_{r'r}\left(G(\omega'){\cal V}_0G^\dagger(\omega')\right)_{rr'},\end{aligned}$$ where $$\begin{aligned}
\label{eq:S102}
  {\cal V}_{0}&=&\sum_{a=1}^N \tilde v_a|-\infty a\rangle\langle -\infty a|,\nonumber\\
  {\cal V}_{L}&=&\sum_{a=1}^N \tilde v_a|\infty a\rangle\langle \infty a|.\end{aligned}$$ Let $\lambda_1(\omega,\omega')$ be the largest eigenvalue of ${\cal M}(\omega,\omega')$ (which obviously is nonnegative). Then by using Eq. (\[eq:S35\]) the following inequality, $$\label{eq:S37}
    \left\|t(\omega)-t(\omega')\right\|_{\rm HS}\leq \sqrt{\lambda_1(\omega,\omega')}\|\omega-\omega'\|,$$ can be readily shown.

Because $M_{rr'}$ decays rapidly as the distance between $r$ and $r'$ increases, we have $$\begin{aligned}
{\rm Tr}{\cal M}^n={\cal O}(M),\quad \forall n\in \mathds{N}.
\label{eq:S42}\end{aligned}$$ Therefore, all the absolute values of the eigenvalues cannot grow with $M$. Thanks to $M\propto N$ we have $$\label{eq:S43}
    \lambda_1(\omega,\omega')={\cal O}(1),$$ from which $$\label{eq:S44}
    \tilde \lambda_1\equiv {\rm sup}_{\omega,\omega'}\lambda_1(\omega,\omega')={\cal O}(1)$$ follows.

Taking the inequality (\[eq:S37\]) and the definition of $\tilde \lambda_1$ \[cf. Eq. (\[eq:S44\])\] into account we obtain $$\label{eq:S39}
    \left\|t(\omega)-t(\omega')\right\|_{\rm HS}\leq \tilde\lambda_1^{1/2}\|\omega-\omega'\|.$$ Combining it with the inequality (\[eq:S30\]) we obtain $$\begin{aligned}
\label{eq:S40}
|W(x;\omega)-W(x;\omega')|
\leq \|W_x\|_{\rm Lip}\|\omega-\omega'\|\end{aligned}$$ and the analytic expression, $$\label{eq:S41}
    \|W_x\|_{\rm Lip}=\frac{\tilde c(x)}{\sqrt{N}},\quad \tilde c(x)\equiv \tilde\lambda_1^{1/2} \|g_x\|_{\rm Lip}.$$ Equations (\[eq:S44\]) and (\[eq:S41\]) show that $$\label{eq:S83}
    \|W_x\|_{\rm Lip}={\cal O}(1/\sqrt{N})$$ and its depth profile is governed by the eigenchannel structure $W_\tau(x)$. Thus the inequality (\[eq:2\]) is justified. In the next section we further provide the numerical confirmation of the scaling (\[eq:S83\]).\
\
[**S3.4 Discussions**]{}\
We have seen that the exact expression of $W_\tau(x)$ is not essential to the proof. What we really used in the proof are the following two properties. (i) The structures $\{W_{\tau_n}(x)\}$ have a universal expression $W_\tau(x)$, which is a function of $\tau$; (ii) $\frac{d}{dz}W_{z^2}(x)$ is bounded. Provided (i) and (ii) are satisfied (according to previous experiences [@Tian15SM] and the results in Sec. S2, this seems to be the case.), even though the exact expression of $W_\tau(x)$ is unknown, we can still repeat the proof and obtain the key results, i.e., the inequality (\[eq:2\]) and the scaling (\[eq:S83\]). Of course we cannot determine the explicit form of $\|g_x\|_{\rm Lip}$ and $\tilde c(x)$ in this case.\
\

![We simulate $W(x;\omega)$ for $10^4$ randomly chosen disorder realizations $\omega$. For both large (a) and small (b) $N$ the distributions of $W(x;\omega)$ at distinct $x$ (histograms) are well fit by (one-sided) Gaussian distributions (dashed lines). The variance scales with $N$ as $N^{-1}$. $L=50$[]{data-label="fig:S2"}](distribution_steady_state.pdf){width="8.7cm"}

![The same as Fig. \[fig:S2\]. The simulated profiles of the steady state $W(x;\omega)$ concentrate strongly around their disorder mean ${\rm E}[W(x;\omega)]$.[]{data-label="fig:S5"}](typicality_W.pdf){width="8.7cm"}

![The depth profiles of $\pi^2Nvar(I_x)$ obtained from simulating $10^4$ randomly chosen incoming waves and $\|I_x\|_{\rm Lip}^2$ obtained from the analytic formula (\[eq:14\]) collapse into a single curve. The blue (red) curve is the (in)coherent part of $\|I_x\|_{\rm Lip}^2$.[]{data-label="fig:S4"}](typicality_variance_Lip.pdf){width="8.7cm"}

[[**S4 Proof and numerical confirmation of the concentration inequality (\[eq:5\])**]{}]{}\
[[**S4.1 Proof of the inequality**]{}]{}\
[*Proof.*]{} We use (Pisier’s) Theorem 0.4. Let $\varpi$ in that theorem be $\omega$, $f(\varpi)$ be $W(\omega)$, and $k=M$. Due to the inequality (\[eq:2\]) and recalling that the disorder ensemble is described by a Gaussian probability measure on $\mathds{R}^M$ of zero mean and variance $\sigma^2$, we find that the conditions of Pisier’s theorem are satisfied. After rescaling we obtain from the theorem the inequality (\[eq:5\]). $\Box$\
\
[[**S4.2 Numerical confirmation**]{}]{}\
We simulate the eigenchannel structures $\{W_{\tau_n}(x)\}$ by using the method described in Sec. S7 for a large number of randomly chosen $\omega$. For each $\omega$ we substitute the simulation results of $\{W_{\tau_n}(x)\}$ into the expression: $$\begin{aligned}
W(x;\omega)=\frac{1}{N}\sum_{n=1}^N W_{\tau_n}(x)
\label{eq:S84}\end{aligned}$$ and obtain $W(x;\omega)$. As shown in Fig. \[fig:S2\], for both large (a) and small (b) $N$ the distributions of $W(x;\omega)$ at distinct $x$ are all well fit by (one-sided) Gaussian distributions. Moreover, the variance of these distributions are found to change with $N$ as (c): $$\label{eq:S85}
    \ln var(W(x;\omega))=-\ln N+\tilde C(x),$$ where $\tilde C(x)$ is some function independent of $N$. Thus simulations confirm the scaling (\[eq:S83\]) and the concentration inequality (\[eq:5\]). In addition, they confirm that the scaling and the inequality hold for both large and small $N$, as long as the medium is diffusive. As shown in Fig. \[fig:S5\], simulations confirm that for both large and small $N$, the profiles of the steady state $W(x;\omega)$ concentrate strongly around their disorder mean ${\rm E}[W(x;\omega)]$ which is a linear decrease, namely, the diffusive steady state.\
\
[**S5 Derivations and numerical studies of Eq. (\[eq:20\])**]{}\
In the first part of this section we derive Eq. (\[eq:20\]). For this purpose we decompose $I_x(c)$ into two parts, $$\begin{aligned}
    I_x(c)&=& I_{inc,x}(c)+I_{coh,x}(c),\label{eq:S70}\\
    I_{inc,x}(c)&\equiv&\sum_{n=1}^N |c_n|^2 W_{\tau_n}(x),\label{eq:S71}\\
    I_{coh,x}(c)&\equiv&\sum_{n\neq n'}^N c_n^*c_{n'} E_{n}^\dagger(x)\cdot E_{n'}(x).
    \label{eq:S72}\end{aligned}$$ The first part, $I_{inc,x}(c)$, accounts for the incoherent contributions of eigenchannels, while the second part, $I_{coh,x}(c)$, for the contributions arising from the coherence between distinct eigenchannels. By using the integral (\[eq:S79\]) in Sec. S8, we find $$\begin{aligned}
\label{eq:S73}
  var(I_x) = \int \left(I_x(c)-W(x;\omega)\right)^2 d\mu\quad\quad\quad\nonumber\\
  = \int \left(\left(I_{inc,x}(c)-W(x;\omega)\right)+I_{coh,x}(c)\right)^2 d\mu\quad\quad\nonumber\\
  = \int \left(I_{inc,x}(c)-W(x;\omega)\right)^2 d\mu+\int I_{coh,x}^2(c) d\mu.\,\end{aligned}$$ This shows that the variance includes the incoherent and coherent contributions, corresponding to the first and second terms, respectively. Substituting Eqs. (\[eq:S71\]) and (\[eq:S72\]) into these two terms, and using the integrals (\[eq:S79\])-(\[eq:S61\]), we obtain $$\begin{aligned}
\label{eq:S74}
  &&\int \left(I_{inc,x}(c)-W(x;\omega)\right)^2 d\mu\nonumber\\
  &=& \frac{1}{N(N+1)}\sum_{n=1}^N\left(W_{\tau_n}(x)-W(x;\omega)\right)^2\end{aligned}$$ and $$\begin{aligned}
\label{eq:S75}
  \int I_{coh,x}^2(c)d\mu=\frac{1}{N(N+1)}\sum_{n\neq n'}^N|E_{n}^\dagger(x)\cdot E_{n'}(x)|^2.\end{aligned}$$ Adding Eqs. (\[eq:S74\]) and (\[eq:S75\]) together we find $$\begin{aligned}
\label{eq:S93}
  var(I_x) &=& \frac{1}{N(N+1)}\bigg(\sum_{n=1}^N\left(W_{\tau_n}(x)-W(x;\omega)\right)^2\nonumber\\
  &&\quad\quad\quad\,\,\,\,\,\,\,\,+\sum_{n\neq n'}^N|E_{n}^\dagger(x)\cdot E_{n'}(x)|^2\bigg).\end{aligned}$$ Comparing it with Eq. (\[eq:14\]) we justify Eq. (\[eq:20\]).

We have provided in Fig. \[fig:1\] simulation results for distinct $N$ and $x$, which confirm the relation (\[eq:20\]). In the second part of this section we perform further numerical studies of $var(I_x)$ and $\|I_x\|_{\rm Lip}^2$. To be specific, we fix $N$ to be $400$, and use the method described in Sec. S7 to simulate wave propagation in a single disordered slab ($L=50$) for $10^4$ randomly chosen incoming waves. The simulated depth profile of $\pi^2Nvar(I_x)$ is shown in Fig. \[fig:S4\]. Furthermore, for the same disordered medium we use the analytic formula (\[eq:14\]) to calculate $\|I_x\|_{\rm Lip}^2$, and its incoherent and coherent parts, respectively. Corresponding profiles are shown in Fig. \[fig:S4\] also. We see that the profiles of $\pi^2Nvar(I_x)$ and $\|I_x\|_{\rm Lip}^2$ are identical, and the incoherent (coherent) part of $\|I_x\|_{\rm Lip}^2$ dominates in the back (front) part of the medium.\
\
[**S6 Derivations of Eq. (\[eq:13\])**]{}\
In most of this work we have focused on the energy density $I_x(c)$ integrated over the transverse coordinate $y$. In this section we consider a generic observable, namely, an operator $\hat O$. It is important that this observable needs not to be an integral over $y$. For simplicity we consider $\hat O$ which is hermitian. Note that this condition is not essential. For a non-hermitian operator its expectation value at a stationary scattering state has an imaginary part in general. In this case we only need to discuss the real and imaginary parts separately, and repeat the analysis below.

Treating $\omega$ as a scattering potential, we apply the scattering theory of waves [@Newton82SM] to Eq. (\[eq:15\]) and find $$\label{eq:S122}
    \hat v_x^{\frac{1}{2}}E(x,y)=\sum_{a=1}^N(t(x)c)_a\varphi^*_a(y).$$ Recall that $$\label{eq:S123}
    c=\sum_{n=1}^N c_n v_n$$ is the incoming current amplitude in the eigenchannel representation, and $\hat v_x\equiv (1-\partial_{\Omega y}^2)^{\frac{1}{2}}$ is a scalar (not vector) operator accounting for the group velocity in the waveguide modes. Using Eqs. (\[eq:S122\]) and (\[eq:S123\]), we find that at the stationary scattering state $E(x,y)$ corresponding to $c$, the expectation value of $\hat O$ is $$\begin{aligned}
\label{eq:S124}
  O(c)\equiv \langle E|\hat O|E\rangle=\sum_{n,n'=1}^N c_n^*c_{n'} \langle E_n|\hat O|E_{n'}\rangle.\end{aligned}$$ Recall that $E_n(x,y)$ is the $2$D wave field of the $n$th eigenchannel. Strictly speaking, like the definition of $I_x(c)$, in the last equality above there are two scalar operators $\hat v_x^{-\frac{1}{2}}$ sandwiching $\hat O$. However, because these two operators account only for an irrelevant overall factor, we omit them hereafter. Naturally, the problem now is: [*For fixed $\omega$, does $O(c)$ exhibit universal behaviors when $c$ varies*]{}? Below we use L$\acute{\rm e}$vy’s lemma to study this problem.

Equation (\[eq:S124\]) defines a real-valued function over $S^{2N-1}$. If the coefficients $\langle E_n|\hat O|E_{n'}\rangle$ do not diverge, (this condition can be achieved readily in physical systems.) then $O(c)$ is continuously differentiable and thus Lipschitz, according to the discussions in Sec. S1.1. With the help of (L$\acute{\rm e}$vy) Lemma 0.5, we obtain the following concentration inequality: $$\label{eq:S125}
    {\rm Pr}\left(\left|O(c)-{\bar O}\right|>\varepsilon\right)\leq 2e^{-
\frac{\delta\varepsilon^2 N}{\|O\|_{{\rm Lip}}^2}},$$ where $\|O\|_{{\rm Lip}}$ is the Lipschitz constant of $O$. This means that $O(c)$ concentrates around its mean: ${\bar O}\equiv\int Od\mu$. The latter can be readily calculated by using the definition (\[eq:S124\]). As a result, $$\label{eq:S126}
    {\bar O}=\frac{1}{N}\sum_{n=1}^N \langle E_n|\hat O|E_n\rangle.$$ The wave-to-wave fluctuations of $O$ is governed by $\|O\|_{{\rm Lip}}$. According to the inequality (\[eq:S126\]), provided $$\label{eq:S127}
    {\bar O}\gg\frac{\|O\|_{{\rm Lip}}}{\sqrt{N}},$$ the concentration effect is strong and the wave-to-wave fluctuations are negligible. Thus we have Eq. (\[eq:13\]). Note that the condition (\[eq:S127\]) is the generalization of the criterion (\[eq:7\]).\
\
[**S7 Method of numerical simulations**]{}\
We simulate the wave propagation by using Eq. (\[eq:15\]), where $\delta \epsilon (x,y)$ is drawn from a uniform distribution over an interval $[-\delta \epsilon_0,\delta \epsilon_0]$. Here $\delta \epsilon_0\in (0,1)$ governs the disorder strength and is set to $0.97$ in simulations. The values of $\delta \epsilon (x,y)$ at distinct spatial points are chosen in the same way and independently. For simulations, Eq. (\[eq:15\]) is discretized on a square grid, with the grid spacing being the inverse wave number in the ideal waveguide, and $L$ is defined as the number of discrete points in $x$-direction. We solve this equation by using the recursive Green’s function method [@MacKinnon85SM]. We first calculate $\langle x=L,y|G|x'=0,y'\rangle$, from which we obtain the transmission matrix $t$. By performing the singular value decomposition of $t$, we obtain $\{v_n\}$ and $\{\tau_n\}$. Next, we calculate $\langle x,y|G|x'=0,y'\rangle$, from which we obtain the matrix $t(x)$. Then the wave field in the interior of the medium is given by $E(x,y)=\sum_{a=1}^N (t(x)c)_a\varphi^*_a(y)$. For a single disorder realization, small oscillations in $x$ are often superposed on non-oscillatory backgrounds. These oscillations occur in the wavelength scale and are unimportant. They are removed by performing the local (in $x$) average over a window with a width of wavelength.

To draw a point $c\equiv(c_1,c_2,\cdots,c_N)$ from a uniform distribution over $ S^{2N-1}$, we use the standard method [@Milman86SM; @Marsaglia72SM] (cf. Sec. S0.1). First, we generate $2N$ independent standard normal random variables, $(a'_n,b'_n)$ with $n=1,2,\cdots, N$. Secondly, we normalize each $a'_n$ ($b'_n$) by $\sqrt{\sum_{n=1}^N((a'_n)^2+(b'_n)^2)}$. Define the normalized $a'_n$ ($b'_n$) as the real (imaginary) part of $c_n$ \[cf. Eq. (\[eq:S58\])\]. We generate a desired random point $c$.\
\
[**S8 Derivations of three integrals**]{}\
In the discussions above we have used the following three integrals: $$\begin{aligned}
  \int c_nc_{n'}^*d\mu &=& \frac{1}{N}\delta_{nn'},\label{eq:S79}\\
  \int |c_n|^4d\mu &=& \frac{2}{N(N+1)},\label{eq:S60}\\
  \int |c_n|^2|c_{n'}|^2d\mu &=& \frac{1}{N(N+1)},\quad n\neq n'.
  \label{eq:S61}\end{aligned}$$ Recall that $\mu$ is the uniform probability measure over $S^{2N-1}$, and $(c_1,c_2,\cdots,c_N)$ are the coordinates of $S^{2N-1}$. In this section we derive these integrals.

The first integral is easy to show. For $n\neq n'$ the left-hand side of Eq. (\[eq:S79\]) vanishes after integrating out the relative phase between $c_n$ and $c_{n'}$. Thus $n,n'$ must be paired, i.e., $n=n'$ for it not to vanish. By the symmetry Eq. (\[eq:S79\]) then follows. The remainder of this section is to derive the last two integrals.

We parametrize $c_n$ by \[cf. Eq. (\[eq:S58\])\]: $$\begin{aligned}
    &&c_n=a_n+ib_n,\quad a_n,b_n\in \mathds{R}\nonumber\\
    &&\sum_{n=1}^N \left(a_n^2+b_n^2\right)=1,
    \label{eq:S62}\end{aligned}$$ where $(a_1,b_2,\cdots,a_N,b_N)$ constitute the real coordinates of $S^{2N-1}$. Substituting it into the left-hand sides of Eqs. (\[eq:S60\]) and (\[eq:S61\]) and using the symmetry, we obtain $$\begin{aligned}
\int |c_n|^4 d\mu &=& 2\left(\int a_1^4d\mu+\int a_1^2b_1^2 d\mu\right),\label{eq:S98}\\
\int |c_n|^2|c_{n'}|^2d\mu &=&4\int a_1^2b_1^2d\mu.
  \label{eq:S99}\end{aligned}$$ So the calculations of the two integrals are reduced to the calculations of the integrals $\int a_1^4d\mu$ and $\int a_1^2b_1^2 d\mu$.

To calculate the last two integrals we define $$\label{eq:S100}
    r_{2n-1}\equiv a_n,\quad r_{2n}\equiv b_n,\quad n=1,2,\cdots,N$$ and parametrize $(r_1,r_2,\cdots,r_{2N-1},r_{2N})$ as follows [@Prudnikov81]: $$\begin{aligned}
r_1&=&\cos\varphi_1,\nonumber\\
r_2&=&\sin\varphi_1\cos\varphi_2,\nonumber\\
r_3&=&\sin\varphi_1\sin\varphi_2\cos\varphi_3,\nonumber\\
&&\cdots\cdots\cdots\cdots\cdots\nonumber\\
r_{2N-1}&=&\sin\varphi_1\sin\varphi_2\cdots\sin\varphi_{2N-2}\cos\varphi_{2N-1},\nonumber\\
r_{2N}&=&\sin\varphi_1\sin\varphi_2\cdots\sin\varphi_{2N-2}\sin\varphi_{2N-1},
\label{eq:S63}\end{aligned}$$ where $$\label{eq:S64}
    0\leq \varphi_1,\varphi_2,\cdots,\varphi_{2N-2}\leq \pi,\quad 0\leq \varphi_{2N-1}<2\pi.$$ Correspondingly, the Jacobian $$\label{eq:S65}
    J=\sin^{2N-2}\varphi_1\sin^{2N-3}\varphi_2\cdots \sin^{2}\varphi_{2N-3}\sin\varphi_{2N-2}.$$ With the substitution of this parametrization we obtain $$\begin{aligned}
  \int a_1^4d\mu &=& \int r_1^4d\mu = \frac{\int_0^{\frac{\pi}{2}}\cos^4\varphi_1\sin^{2N-2}\varphi_1 d\varphi_1}{\int_0^{\frac{\pi}{2}}\sin^{2N-2}\varphi_1 d\varphi_1}\nonumber\\
  &=& \frac{B\left(\frac{2N-1}{2},\frac{5}{2}\right)}{B\left(\frac{2N-1}{2},\frac{1}{2}\right)},
  \label{eq:S66}\end{aligned}$$ where $B(x,y)$ is the beta function. Using the identity: $$\label{eq:S67}
    B(x,y)=\frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)},$$ with $\Gamma(x)$ being the gamma function, we reduce Eq. (\[eq:S66\]) to $$\begin{aligned}
  \int a_1^4d\mu=\frac{3}{4N(N+1)}.
  \label{eq:S94}\end{aligned}$$ With the substitution of this parametrization we also obtain $$\begin{aligned}
  &&\int a_1^2b_1^2d\mu=\int r_1^2r_2^2d\mu \nonumber\\
  &=& \frac{\int_0^{\pi}\cos^2\varphi_1\sin^{2N}\varphi_1 d\varphi_1}{\int_0^{\pi}\sin^{2N-2}\varphi_1 d\varphi_1}
  \frac{\int_0^{\pi}\cos^2\varphi_2\sin^{2N-3}\varphi_2 d\varphi_2}{\int_0^{\pi}\sin^{2N-3}\varphi_2 d\varphi_2}\nonumber\\
  &=& \frac{B\left(\frac{2N+1}{2},\frac{3}{2}\right)}{B\left(\frac{2N-1}{2},\frac{1}{2}\right)}
  \frac{B\left(\frac{2N-2}{2},\frac{3}{2}\right)}{B\left(\frac{2N-2}{2},\frac{1}{2}\right)}.
  \label{eq:68}\end{aligned}$$ Using the identity (\[eq:S67\]) we reduce it to $$\begin{aligned}
\int a_1^2b_1^2d\mu &=&\frac{1}{4N(N+1)}.
\label{eq:S95}\end{aligned}$$ Substituting Eqs. (\[eq:S94\]) and (\[eq:S95\]) into Eqs. (\[eq:S98\]) and (\[eq:S99\]) we justify Eqs. (\[eq:S60\]) and (\[eq:S61\]).

S. Boucheron, G. Lugosi, and P. Massart, [*Concentration Inequalities - A Nonasymptotic Theory of Independence*]{} (Oxford University Press, UK, 2013).

V. D. Milman and G. Schechtman, [*Asymptotic Theory of Finite-Dimensional Normed Spaces*]{} (Springer, Berlin, Germany, 1986). M. Ledoux, [*The Concentration of Measure Phenomenon*]{} (AMS, Providence, 2001).

M. Davy, Z. Shi, J. Park, C. Tian, and A. Z. Genack, [ Nat. Commun.]{} **6**, 6893 (2015).

O. N. Dorokhov, [Solid State Commun.]{} **51**, 381 (1984).

P. A. Mello, P. Pereyra, and N. Kumar, [Ann. Phys.]{} **181**, 290 (1988).

U. Haagerup and S. Thorbjornsen, Expo. Math. **21**, 293 (2003).

R. G. Newton, [*Scattering Theory of Waves and Particles*]{} 2nd Ed. (Springer, New York, 1982).

A. MacKinnon, [Z. Phys. B]{} **59**, 385 (1985).

G. Marsaglia, [Ann. Math. Stat.]{} **43**, 645 (1972).

A. P. Prudnikov, Yu. A. Brychkov, and O. I. Marichev, [*Integraly i ryady, tom 1. E’lementarnye funkcii*]{} (Nauka, Moscow, 1981).
