Anisotropic Field Theory of Wave Transmission Statistics in Disordered Media

David Gaspard [email protected]    Arthur Goetschy [email protected] Institut Langevin, ESPCI Paris, PSL University, CNRS, 75005 Paris, France
(November 18, 2024)
Abstract

We present a field theory to characterize the mean distribution of transmission eigenvalues for coherent wave propagation through a disordered medium. Unlike the Dorokhov-Mello-Pereyra-Kumar (DMPK) theory, our approach does not rely on the isotropy hypothesis, which assumes uniform field intensity across scattering channels. As a result, it accurately predicts the transmission eigenvalue distribution both in the diffusive and quasiballistic regime. In addition, it proves to be more versatile than the DMPK theory, allowing for the incorporation of physical effects such as absorption or partial channel control with minimal adjustments to the formalism. Although these effects are commonly encountered in experiments on complex media, until now there has been no ab initio theory addressing them. Our predictions are numerically validated against transmission eigenvalue distributions computed from the microscopic wave equation.

Wave propagation; Transmission matrix; Disordered media; Complex media; Transmission eigenvalues; Full counting statistics; Eilenberger equation; Usadel equation; Nonlinear sigma model; Wavefront shaping; Quasiballistic regime; Diffusive regime; Absorption; Partial channel control; Dorokhov-Mello-Pereyra-Kumar theory; Recursive Green’s algorithm;

Coherent wave propagation in disordered media arises across various fields of science and technology, from optical and acoustic imaging to studies of coherent electron transport in mesoscopic conductors [1, 2, 3]. In many situations it is relevant to describe this propagation using the transmission coefficient matrix, 𝗍𝗍\mathsf{t}sansserif_t, which connects the wavefronts at the system input to the wavefronts at the output. The eigenvalues of 𝗍𝗍superscript𝗍𝗍\mathsf{t}^{\dagger}\mathsf{t}sansserif_t start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT sansserif_t, always contained between 00 and 1111, determine the transmittances of the corresponding eigenmodes. It has been known for decades that, in the diffusive regime, the distribution of these eigenvalues follows the bimodal law, ρ(T)=T¯/(2T1T)𝜌𝑇¯𝑇2𝑇1𝑇\rho(T)=\bar{T}/(2T\sqrt{1-T})italic_ρ ( italic_T ) = over¯ start_ARG italic_T end_ARG / ( 2 italic_T square-root start_ARG 1 - italic_T end_ARG ), T¯¯𝑇\bar{T}over¯ start_ARG italic_T end_ARG being the mean transmission [4, 5, 6, 7]. Surprisingly, this distribution allows the existence of quasitransparent channels (near T=1𝑇1T=1italic_T = 1) despite the on-average opacity of the disordered medium (T¯1much-less-than¯𝑇1\bar{T}\ll 1over¯ start_ARG italic_T end_ARG ≪ 1). This rather unusual property shows that the incident wavefront can be shaped to minimize the impact of scattering.

The theoretical study of wave transmission through disordered systems was marked by the parallel development of two competing theories: the nonlinear sigma model [8, 9, 10, 11, 12, 13, 14, 15], and the Dorokhov-Mello-Pereyra-Kumar (DMPK) theory [16, 17, 18, 19, 7]. The nonlinear sigma model, on one hand, refers to a class of field-theoretical models defined by constraining the field to a specific target manifold. These models are typically obtained by a systematic approach for averaging over the disorder initially proposed by Wegner and Schäfer [8, 9] and previous authors [20, 21, 22, 23] in the framework of electron transport. This approach originally exploited the replica method to achieve the disorder averaging but was soon endowed of a supersymmetric space by Efetov [10, 11, 12, 13] to avoid the use of replicas. In the following years, many authors used the supersymmetric method to study the transport properties of waves in disordered media especially in the strongly localized regime [24, 25, 26, 27, 28, 29, 30, 31, 32]. Nazarov was the first to derive the transmission eigenvalue distribution in the nonlocalized diffusive regime using a nonlinear sigma model without supersymmetry [33, 34]. His approach, apparently unrelated to that of Efetov, is based on an analogy with kinetic equations for type-II superconductors [35, 36, 37, 38, *Larkin1973, *Larkin1975, *Larkin1977, 42] where the disorder caused by impurities is known to play an important role. Nazarov’s approach did not gain widespread adoption in the literature on disordered systems, probably due to shortcuts in mathematical reasoning. However, given its universality, this approach holds potential applications beyond electron transport, including fields like optics and acoustics. The primary challenge lies in extending the approach to account for quasiballistic transport, absorption, and incomplete channel control—regimes that Nazarov’s method originally did not address.

The DMPK theory, on the other hand, proceeds very differently. It is based on a perturbative treatment of an infinitesimal slice of disordered medium. This approach provides a Fokker-Planck-type equation, known as the DMPK equation, for the joint distribution of transmission eigenvalues. The DMPK theory became very successful for the simplicity, elegance, and predictive power of its formalism [7]. However, it has notable limitations in scenarios of high experimental interest.

A first limitation of the DMPK approach comes from the isotropy hypothesis [17, 7] which assumes that the field intensity is uniformly distributed across the different channels due to scattering by each infinitesimal slice of the disordered medium. It turns out that this isotropy hypothesis is only valid in the multiple-scattering regime and not in the quasiballistic regime. This distorts the predictions of the DMPK theory in the quasiballistic regime, which is experimentally highly relevant. Incidentally, nonlinear sigma models face the same limitation due to the equivalence between the two theories [43, 44, 45]. Another limitation of the DMPK theory reported by Brouwer [46] comes from the inherent difficulties in including absorption, a crucial phenomenon in optics and acoustics. Last but not least, the DMPK theory is focused on the transmission matrix and cannot be modified to study the statistics of other observables such as the field intensity inside the system for example [47, 48], or to account for typical experimental configurations such as partial channel control which results from the limited numerical aperture of measuring devices or the finite spatial extent of the incident beam [49, 50, 51]. All these limitations of the DMPK theory can be overcome by the anisotropic field theory presented in this Letter and derived in the companion paper [52].

To set the framework of our approach, we consider a disordered region of thickness L𝐿Litalic_L, characterized by the random potential U(𝐫)𝑈𝐫U(\boldsymbol{\mathrm{r}})italic_U ( bold_r ), contained within a perfect waveguide. We assume that the wave is emitted in the transverse mode χn(𝐲)subscript𝜒𝑛𝐲\chi_{n}(\boldsymbol{\mathrm{y}})italic_χ start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ( bold_y ) from a fictitious surface with longitudinal coordinate xasubscript𝑥ax_{\mathrm{a}}italic_x start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT and measured in the mode χm(𝐲)subscript𝜒𝑚𝐲\chi_{m}(\boldsymbol{\mathrm{y}})italic_χ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( bold_y ) at another surface at xbsubscript𝑥bx_{\mathrm{b}}italic_x start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT. It is not necessary that these surfaces lie outside the disordered region [33, 52]. The transmission matrix coefficient tmnsubscript𝑡𝑚𝑛t_{mn}italic_t start_POSTSUBSCRIPT italic_m italic_n end_POSTSUBSCRIPT is then determined by the Fisher and Lee relation [53, 54, 55],

tmn=2ikx,mkx,nχm,xb|G^+|χn,xa,subscript𝑡𝑚𝑛2isubscript𝑘𝑥𝑚subscript𝑘𝑥𝑛quantum-operator-productsubscript𝜒𝑚subscript𝑥bsuperscript^𝐺subscript𝜒𝑛subscript𝑥at_{mn}=2\mathrm{i}\sqrt{k_{x,m}k_{x,n}}\left\langle\chi_{m},x_{\mathrm{b}}% \right|\hat{G}^{+}\left|\chi_{n},x_{\mathrm{a}}\right\rangle\>,italic_t start_POSTSUBSCRIPT italic_m italic_n end_POSTSUBSCRIPT = 2 roman_i square-root start_ARG italic_k start_POSTSUBSCRIPT italic_x , italic_m end_POSTSUBSCRIPT italic_k start_POSTSUBSCRIPT italic_x , italic_n end_POSTSUBSCRIPT end_ARG ⟨ italic_χ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT , italic_x start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT | over^ start_ARG italic_G end_ARG start_POSTSUPERSCRIPT + end_POSTSUPERSCRIPT | italic_χ start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT , italic_x start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT ⟩ , (1)

where kx,nsubscript𝑘𝑥𝑛k_{x,n}italic_k start_POSTSUBSCRIPT italic_x , italic_n end_POSTSUBSCRIPT is the longitudinal wavenumber of mode n𝑛nitalic_n. In Eq. (1) and those that follow, we use Dirac’s notation: hats denote operators acting in the position space. The waveguide supports Npsubscript𝑁pN_{\rm p}italic_N start_POSTSUBSCRIPT roman_p end_POSTSUBSCRIPT propagating modes, but only Nasubscript𝑁aN_{\mathrm{a}}italic_N start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT of them are controlled at the input, and Nbsubscript𝑁bN_{\mathrm{b}}italic_N start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT at the output. Therefore, the dimensions of 𝗍𝗍\mathsf{t}sansserif_t are Nb×Nasubscript𝑁bsubscript𝑁aN_{\mathrm{b}}\times N_{\mathrm{a}}italic_N start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT × italic_N start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT. In Eq. (1), G^+superscript^𝐺\hat{G}^{+}over^ start_ARG italic_G end_ARG start_POSTSUPERSCRIPT + end_POSTSUPERSCRIPT is the Green’s operator defined by

G^±=[^𝐫2+k2±iεU(𝐫^)]1,superscript^𝐺plus-or-minussuperscriptdelimited-[]plus-or-minussubscriptsuperscript^2𝐫superscript𝑘2i𝜀𝑈^𝐫1\hat{G}^{\pm}=\big{[}\hat{\nabla}^{2}_{\boldsymbol{\mathrm{r}}}+k^{2}\pm% \mathrm{i}\varepsilon-U(\hat{\boldsymbol{\mathrm{r}}})\big{]}^{-1}\>,over^ start_ARG italic_G end_ARG start_POSTSUPERSCRIPT ± end_POSTSUPERSCRIPT = [ over^ start_ARG ∇ end_ARG start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT bold_r end_POSTSUBSCRIPT + italic_k start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ± roman_i italic_ε - italic_U ( over^ start_ARG bold_r end_ARG ) ] start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT , (2)

where k=2π/λ𝑘2𝜋𝜆k=2\pi/\lambdaitalic_k = 2 italic_π / italic_λ is the wavenumber, and ε𝜀\varepsilonitalic_ε the Laplace variable associated with time (tεsubscript𝑡𝜀\partial_{t}\leftrightarrow\varepsilon∂ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ↔ italic_ε). The random potential U(𝐫)𝑈𝐫U(\boldsymbol{\mathrm{r}})italic_U ( bold_r ) in Eq. (2) is normally distributed according to 𝒫[U(𝐫)]exp(dd𝐫U(𝐫)22α(𝐫))proportional-to𝒫delimited-[]𝑈𝐫subscriptsuperscript𝑑differential-d𝐫𝑈superscript𝐫22𝛼𝐫\textstyle\mathcal{P}[U(\boldsymbol{\mathrm{r}})]\propto\exp(-\int_{\mathbb{R}% ^{d}}\mathop{}\!\mathrm{d}\boldsymbol{\mathrm{r}}\frac{U(\boldsymbol{\mathrm{r% }})^{2}}{2\alpha(\boldsymbol{\mathrm{r}})})caligraphic_P [ italic_U ( bold_r ) ] ∝ roman_exp ( - ∫ start_POSTSUBSCRIPT blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT end_POSTSUBSCRIPT roman_d bold_r divide start_ARG italic_U ( bold_r ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_α ( bold_r ) end_ARG ). The disorder strength α𝛼\alphaitalic_α determines the value of the scattering mean free path, s=k/πναsubscripts𝑘𝜋𝜈𝛼\ell_{\rm s}=k/\pi\nu\alpharoman_ℓ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT = italic_k / italic_π italic_ν italic_α, where ν𝜈\nuitalic_ν is the local density of states per unit k2superscript𝑘2k^{2}italic_k start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT.

The central quantity is the distribution of transmission eigenvalues defined by

ρ(T)=def1NaTrδ(T𝗍𝗍)=1πT2ImF(1T+0i),𝜌𝑇def1subscript𝑁adelimited-⟨⟩Tr𝛿𝑇superscript𝗍𝗍1𝜋superscript𝑇2Im𝐹1𝑇0i\rho(T)\overset{\mathrm{def}}{=}\frac{1}{N_{\mathrm{a}}}\left\langle% \operatorname{Tr}\delta(T-\mathsf{t}^{\dagger}\mathsf{t})\right\rangle=\frac{1% }{\pi T^{2}}\operatorname{Im}F(\tfrac{1}{T}+0\mathrm{i})\>,italic_ρ ( italic_T ) overroman_def start_ARG = end_ARG divide start_ARG 1 end_ARG start_ARG italic_N start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT end_ARG ⟨ roman_Tr italic_δ ( italic_T - sansserif_t start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT sansserif_t ) ⟩ = divide start_ARG 1 end_ARG start_ARG italic_π italic_T start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG roman_Im italic_F ( divide start_ARG 1 end_ARG start_ARG italic_T end_ARG + 0 roman_i ) , (3)

where delimited-⟨⟩\left\langle\cdot\right\rangle⟨ ⋅ ⟩ denotes the average over the random potential U(𝐫)𝑈𝐫U(\boldsymbol{\mathrm{r}})italic_U ( bold_r ). As indicated by the second equality of Eq. (3), the distribution ρ(T)𝜌𝑇\rho(T)italic_ρ ( italic_T ) can be extracted from the generating function:

F(γ)=1NaddγlnZ(γ),𝐹𝛾1subscript𝑁add𝛾delimited-⟨⟩𝑍𝛾F(\gamma)=\frac{1}{N_{\mathrm{a}}}\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!% \mathrm{d}\gamma}\left\langle\ln Z(\gamma)\right\rangle\>,italic_F ( italic_γ ) = divide start_ARG 1 end_ARG start_ARG italic_N start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT end_ARG divide start_ARG roman_d end_ARG start_ARG roman_d italic_γ end_ARG ⟨ roman_ln italic_Z ( italic_γ ) ⟩ , (4)

where Z(γ)=𝒩det(1γ𝗍𝗍)1𝑍𝛾𝒩superscript1𝛾superscript𝗍𝗍1Z(\gamma)=\mathcal{N}\det(1-\gamma\mathsf{t}^{\dagger}\mathsf{t})^{-1}italic_Z ( italic_γ ) = caligraphic_N roman_det ( 1 - italic_γ sansserif_t start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT sansserif_t ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT, with 𝒩𝒩\mathcal{N}caligraphic_N a γ𝛾\gammaitalic_γ-independent arbitrary prefactor. It is useful to write this partition function in the form (see Sec. II B of the companion paper [52]):

Z(γ)=𝒩det(1γG^+K^aG^K^b)1=𝒟[ϕ(𝐫)]eid𝐫ϕ(𝐫)[𝖶^0U(𝐫)]ϕ(𝐫),𝑍𝛾𝒩superscript1𝛾superscript^𝐺subscript^𝐾asuperscript^𝐺subscript^𝐾b1𝒟delimited-[]bold-italic-ϕ𝐫superscripteidifferential-d𝐫superscriptbold-italic-ϕ𝐫delimited-[]subscript^𝖶0𝑈𝐫bold-italic-ϕ𝐫\begin{split}Z(\gamma)&=\mathcal{N}\det(1-\gamma\hat{G}^{+}\hat{K}_{\mathrm{a}% }\hat{G}^{-}\hat{K}_{\mathrm{b}})^{-1}\\ &=\int\mathcal{D}[\boldsymbol{\mathrm{\phi}}(\boldsymbol{\mathrm{r}})]\mathop{% }\!\mathrm{e}^{\mathrm{i}\int\mathop{}\!\mathrm{d}\boldsymbol{\mathrm{r}}\>% \boldsymbol{\mathrm{\phi}}^{\dagger}(\boldsymbol{\mathrm{r}})[\hat{\mathsf{W}}% _{0}-U(\boldsymbol{\mathrm{r}})]\boldsymbol{\mathrm{\phi}}(\boldsymbol{\mathrm% {r}})}\>,\end{split}start_ROW start_CELL italic_Z ( italic_γ ) end_CELL start_CELL = caligraphic_N roman_det ( 1 - italic_γ over^ start_ARG italic_G end_ARG start_POSTSUPERSCRIPT + end_POSTSUPERSCRIPT over^ start_ARG italic_K end_ARG start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT over^ start_ARG italic_G end_ARG start_POSTSUPERSCRIPT - end_POSTSUPERSCRIPT over^ start_ARG italic_K end_ARG start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL = ∫ caligraphic_D [ bold_italic_ϕ ( bold_r ) ] roman_e start_POSTSUPERSCRIPT roman_i ∫ roman_d bold_r bold_italic_ϕ start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ( bold_r ) [ over^ start_ARG sansserif_W end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - italic_U ( bold_r ) ] bold_italic_ϕ ( bold_r ) end_POSTSUPERSCRIPT , end_CELL end_ROW (5)

where ϕ(𝐫)=[ϕ1(𝐫),ϕ2(𝐫)]bold-italic-ϕ𝐫superscriptsubscriptitalic-ϕ1𝐫subscriptitalic-ϕ2𝐫\boldsymbol{\mathrm{\phi}}(\boldsymbol{\mathrm{r}})={[\phi_{1}(\boldsymbol{% \mathrm{r}}),\phi_{2}(\boldsymbol{\mathrm{r}})]}^{\intercal}bold_italic_ϕ ( bold_r ) = [ italic_ϕ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( bold_r ) , italic_ϕ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ( bold_r ) ] start_POSTSUPERSCRIPT ⊺ end_POSTSUPERSCRIPT is a fictitious two-component complex field, and 𝖶^0subscript^𝖶0\hat{\mathsf{W}}_{0}over^ start_ARG sansserif_W end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT is the disorder-free matrix Hamiltonian

𝖶^0=(^𝐫2+k2+iεγK^aK^b^𝐫2+k2iε),subscript^𝖶0matrixsubscriptsuperscript^2𝐫superscript𝑘2i𝜀𝛾subscript^𝐾asubscript^𝐾bsubscriptsuperscript^2𝐫superscript𝑘2i𝜀\hat{\mathsf{W}}_{0}=\begin{pmatrix}\hat{\nabla}^{2}_{\boldsymbol{\mathrm{r}}}% +k^{2}+\mathrm{i}\varepsilon&\gamma\hat{K}_{\mathrm{a}}\\ \hat{K}_{\mathrm{b}}&\hat{\nabla}^{2}_{\boldsymbol{\mathrm{r}}}+k^{2}-\mathrm{% i}\varepsilon\end{pmatrix}\>,over^ start_ARG sansserif_W end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = ( start_ARG start_ROW start_CELL over^ start_ARG ∇ end_ARG start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT bold_r end_POSTSUBSCRIPT + italic_k start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + roman_i italic_ε end_CELL start_CELL italic_γ over^ start_ARG italic_K end_ARG start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL over^ start_ARG italic_K end_ARG start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT end_CELL start_CELL over^ start_ARG ∇ end_ARG start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT bold_r end_POSTSUBSCRIPT + italic_k start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - roman_i italic_ε end_CELL end_ROW end_ARG ) , (6)

K^asubscript^𝐾a\hat{K}_{\mathrm{a}}over^ start_ARG italic_K end_ARG start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT being a current density operator defined by

K^a=Θ^ai^xδ(x^xa)δ(x^xa)Θ^ai^x,subscript^𝐾asubscript^Θaisubscript^𝑥𝛿^𝑥subscript𝑥a𝛿^𝑥subscript𝑥asubscript^Θaisubscript^𝑥\hat{K}_{\mathrm{a}}=-\hat{\Theta}_{\mathrm{a}}\mathrm{i}\hat{\nabla}_{x}% \delta(\hat{x}-x_{\mathrm{a}})-\delta(\hat{x}-x_{\mathrm{a}})\hat{\Theta}_{% \mathrm{a}}\mathrm{i}\hat{\nabla}_{x}\>,over^ start_ARG italic_K end_ARG start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT = - over^ start_ARG roman_Θ end_ARG start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT roman_i over^ start_ARG ∇ end_ARG start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT italic_δ ( over^ start_ARG italic_x end_ARG - italic_x start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT ) - italic_δ ( over^ start_ARG italic_x end_ARG - italic_x start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT ) over^ start_ARG roman_Θ end_ARG start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT roman_i over^ start_ARG ∇ end_ARG start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT , (7)

where Θ^a=Θ(ma𝛀^)subscript^ΘaΘsubscript𝑚anormsubscript^𝛀perpendicular-to\hat{\Theta}_{\mathrm{a}}=\Theta(m_{\mathrm{a}}-\|\hat{\boldsymbol{\mathrm{% \Omega}}}_{\perp}\|)over^ start_ARG roman_Θ end_ARG start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT = roman_Θ ( italic_m start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT - ∥ over^ start_ARG bold_Ω end_ARG start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT ∥ ) with 𝛀^=1ik^𝐲subscript^𝛀perpendicular-to1i𝑘subscript^bold-∇𝐲\hat{\boldsymbol{\mathrm{\Omega}}}_{\perp}=\tfrac{1}{\mathrm{i}k}\hat{% \boldsymbol{\mathrm{\nabla}}}_{\boldsymbol{\mathrm{y}}}over^ start_ARG bold_Ω end_ARG start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG roman_i italic_k end_ARG over^ start_ARG bold_∇ end_ARG start_POSTSUBSCRIPT bold_y end_POSTSUBSCRIPT, and similarly for K^bsubscript^𝐾b\hat{K}_{\mathrm{b}}over^ start_ARG italic_K end_ARG start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT. The Heaviside step function ΘΘ\Thetaroman_Θ has the effect of reducing the number of modes under control when ma<1subscript𝑚a1m_{\mathrm{a}}<1italic_m start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT < 1 or mb<1subscript𝑚b1m_{\mathrm{b}}<1italic_m start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT < 1. The parameters masubscript𝑚am_{\mathrm{a}}italic_m start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT and mbsubscript𝑚bm_{\mathrm{b}}italic_m start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT are the numerical apertures of the input and output ports, respectively. In two dimensions, they are given by ma=Na/Npsubscript𝑚asubscript𝑁asubscript𝑁pm_{\mathrm{a}}=N_{\mathrm{a}}/N_{\rm p}italic_m start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT = italic_N start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT / italic_N start_POSTSUBSCRIPT roman_p end_POSTSUBSCRIPT and mb=Nb/Npsubscript𝑚bsubscript𝑁bsubscript𝑁pm_{\mathrm{b}}=N_{\mathrm{b}}/N_{\rm p}italic_m start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT = italic_N start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT / italic_N start_POSTSUBSCRIPT roman_p end_POSTSUBSCRIPT [49].

The average over the disorder in Eq. (4) can be calculated by the standard replica method [14], lnZ=limR0RZRdelimited-⟨⟩𝑍subscript𝑅0subscript𝑅delimited-⟨⟩superscript𝑍𝑅\textstyle\left\langle\ln Z\right\rangle=\lim_{R\rightarrow 0}\partial_{R}% \left\langle Z^{R}\right\rangle⟨ roman_ln italic_Z ⟩ = roman_lim start_POSTSUBSCRIPT italic_R → 0 end_POSTSUBSCRIPT ∂ start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT ⟨ italic_Z start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT ⟩, where R𝑅Ritalic_R is the number of replicas. Introducing the replicated field, denoted as 𝚽(𝐫)=[ϕ1(𝐫),,ϕR(𝐫)]𝚽𝐫superscriptsubscriptbold-italic-ϕ1𝐫subscriptbold-italic-ϕ𝑅𝐫\boldsymbol{\mathrm{\Phi}}(\boldsymbol{\mathrm{r}})={[\boldsymbol{\mathrm{\phi% }}_{1}(\boldsymbol{\mathrm{r}}),\ldots,\boldsymbol{\mathrm{\phi}}_{R}(% \boldsymbol{\mathrm{r}})]}^{\intercal}bold_Φ ( bold_r ) = [ bold_italic_ϕ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( bold_r ) , … , bold_italic_ϕ start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT ( bold_r ) ] start_POSTSUPERSCRIPT ⊺ end_POSTSUPERSCRIPT, yields

ZR=𝒟[U]𝒟[𝚽]ed𝐫(U22α+i𝚽(𝖶^0U)𝚽).delimited-⟨⟩superscript𝑍𝑅𝒟delimited-[]𝑈𝒟delimited-[]𝚽superscriptedifferential-d𝐫superscript𝑈22𝛼isuperscript𝚽subscript^𝖶0𝑈𝚽\left\langle Z^{R}\right\rangle=\int\mathcal{D}[U]\mathcal{D}[\boldsymbol{% \mathrm{\Phi}}]\mathop{}\!\mathrm{e}^{\int\mathop{}\!\mathrm{d}\boldsymbol{% \mathrm{r}}\left(-\frac{U^{2}}{2\alpha}+\mathrm{i}\boldsymbol{\mathrm{\Phi}}^{% \dagger}(\hat{\mathsf{W}}_{0}-U)\boldsymbol{\mathrm{\Phi}}\right)}\>.⟨ italic_Z start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT ⟩ = ∫ caligraphic_D [ italic_U ] caligraphic_D [ bold_Φ ] roman_e start_POSTSUPERSCRIPT ∫ roman_d bold_r ( - divide start_ARG italic_U start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_α end_ARG + roman_i bold_Φ start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ( over^ start_ARG sansserif_W end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - italic_U ) bold_Φ ) end_POSTSUPERSCRIPT . (8)

The integral over U(𝐫)𝑈𝐫U(\boldsymbol{\mathrm{r}})italic_U ( bold_r ) being Gaussian, we get

ZR=𝒟[𝚽]ed𝐫(i𝚽𝖶^0𝚽α2(𝚽𝚽)2).delimited-⟨⟩superscript𝑍𝑅𝒟delimited-[]𝚽superscriptedifferential-d𝐫isuperscript𝚽subscript^𝖶0𝚽𝛼2superscriptsuperscript𝚽𝚽2\left\langle Z^{R}\right\rangle=\int\mathcal{D}[\boldsymbol{\mathrm{\Phi}}]% \mathop{}\!\mathrm{e}^{\int\mathop{}\!\mathrm{d}\boldsymbol{\mathrm{r}}\left(% \mathrm{i}\boldsymbol{\mathrm{\Phi}}^{\dagger}\hat{\mathsf{W}}_{0}\boldsymbol{% \mathrm{\Phi}}-\frac{\alpha}{2}(\boldsymbol{\mathrm{\Phi}}^{\dagger}% \boldsymbol{\mathrm{\Phi}})^{2}\right)}\>.⟨ italic_Z start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT ⟩ = ∫ caligraphic_D [ bold_Φ ] roman_e start_POSTSUPERSCRIPT ∫ roman_d bold_r ( roman_i bold_Φ start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT over^ start_ARG sansserif_W end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT bold_Φ - divide start_ARG italic_α end_ARG start_ARG 2 end_ARG ( bold_Φ start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT bold_Φ ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) end_POSTSUPERSCRIPT . (9)

The Lagrangian in the exponential argument possesses a 𝚽4superscript𝚽4\boldsymbol{\mathrm{\Phi}}^{4}bold_Φ start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT term giving rise to Goldstone modes at the saddle points [9]. Since these modes vary in position space much more slowly than the wavelength, they can be approximated semiclassically as we will see shortly. The 𝚽4superscript𝚽4\boldsymbol{\mathrm{\Phi}}^{4}bold_Φ start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT term can be eliminated by a Hubbard-Stratonovich transformation [13] which consists in the introduction of an auxiliary matrix field of dimensions 2R×2R2𝑅2𝑅2R\times 2R2 italic_R × 2 italic_R denoted 𝖰(𝐫)𝖰𝐫\mathsf{Q}(\boldsymbol{\mathrm{r}})sansserif_Q ( bold_r ):

ZR=𝒟[𝖰]ed𝐫(α2Tr[𝖰(𝐫)2]𝐫|Trln(𝖶^0+iα𝖰(𝐫^))|𝐫).delimited-⟨⟩superscript𝑍𝑅𝒟delimited-[]𝖰superscriptedifferential-d𝐫𝛼2Tr𝖰superscript𝐫2quantum-operator-product𝐫Trsubscript^𝖶0i𝛼𝖰^𝐫𝐫\left\langle Z^{R}\right\rangle=\int\mathcal{D}[\mathsf{Q}]\mathop{}\!\mathrm{% e}^{\int\mathop{}\!\mathrm{d}\boldsymbol{\mathrm{r}}\left(\frac{\alpha}{2}% \operatorname{Tr}[\mathsf{Q}(\boldsymbol{\mathrm{r}})^{2}]-\left\langle% \boldsymbol{\mathrm{r}}\right|\operatorname{Tr}\ln(\hat{\mathsf{W}}_{0}+% \mathrm{i}\alpha\mathsf{Q}(\hat{\boldsymbol{\mathrm{r}}}))\left|\boldsymbol{% \mathrm{r}}\right\rangle\right)}\>.⟨ italic_Z start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT ⟩ = ∫ caligraphic_D [ sansserif_Q ] roman_e start_POSTSUPERSCRIPT ∫ roman_d bold_r ( divide start_ARG italic_α end_ARG start_ARG 2 end_ARG roman_Tr [ sansserif_Q ( bold_r ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ] - ⟨ bold_r | roman_Tr roman_ln ( over^ start_ARG sansserif_W end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + roman_i italic_α sansserif_Q ( over^ start_ARG bold_r end_ARG ) ) | bold_r ⟩ ) end_POSTSUPERSCRIPT . (10)

It is possible to simplify the Lagrangian in Eq. (10) by assuming the diffusive regime (Lsmuch-greater-than𝐿subscriptsL\gg\ell_{\rm s}italic_L ≫ roman_ℓ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT). This would lead to a nonlinear sigma model. However, we do not make this assumption because it is equivalent to the isotropy hypothesis and is approximate. Instead, we consider the saddle-point equation of the Lagrangian in Eq. (10), which has the form of a nonlinear wave equation for the matrix Green’s function Γ(𝐫,𝐫)sans-serif-Γ𝐫superscript𝐫\mathsf{\Gamma}(\boldsymbol{\mathrm{r}},\boldsymbol{\mathrm{r}}^{\prime})sansserif_Γ ( bold_r , bold_r start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ),

[𝐫2+k2+iεΛ3+iα𝖰(𝐫)+γaK^aΛ++γbK^bΛ]Γ(𝐫,𝐫)=i𝟣2δ(𝐫𝐫),delimited-[]subscriptsuperscript2𝐫superscript𝑘2i𝜀subscriptsans-serif-Λ3i𝛼𝖰𝐫subscript𝛾asubscript^𝐾asubscriptsans-serif-Λsubscript𝛾bsubscript^𝐾bsubscriptsans-serif-Λsans-serif-Γ𝐫superscript𝐫isubscript12𝛿𝐫superscript𝐫\begin{split}\Big{[}&\nabla^{2}_{\boldsymbol{\mathrm{r}}}+k^{2}+\mathrm{i}% \varepsilon\mathsf{\Lambda}_{3}+\mathrm{i}\alpha\mathsf{Q}(\boldsymbol{\mathrm% {r}})\\ &+\gamma_{\mathrm{a}}\hat{K}_{\mathrm{a}}\mathsf{\Lambda}_{+}+\gamma_{\mathrm{% b}}\hat{K}_{\mathrm{b}}\mathsf{\Lambda}_{-}\Big{]}\mathsf{\Gamma}(\boldsymbol{% \mathrm{r}},\boldsymbol{\mathrm{r}}^{\prime})=\mathrm{i}\mathsf{1}_{2}\delta(% \boldsymbol{\mathrm{r}}-\boldsymbol{\mathrm{r}}^{\prime})\>,\end{split}start_ROW start_CELL [ end_CELL start_CELL ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT bold_r end_POSTSUBSCRIPT + italic_k start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + roman_i italic_ε sansserif_Λ start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT + roman_i italic_α sansserif_Q ( bold_r ) end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL + italic_γ start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT over^ start_ARG italic_K end_ARG start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT sansserif_Λ start_POSTSUBSCRIPT + end_POSTSUBSCRIPT + italic_γ start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT over^ start_ARG italic_K end_ARG start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT sansserif_Λ start_POSTSUBSCRIPT - end_POSTSUBSCRIPT ] sansserif_Γ ( bold_r , bold_r start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ) = roman_i sansserif_1 start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT italic_δ ( bold_r - bold_r start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ) , end_CELL end_ROW (11)

with the self-consistency condition 𝖰(𝐫)=Γ(𝐫,𝐫)𝖰𝐫sans-serif-Γ𝐫𝐫\mathsf{Q}(\boldsymbol{\mathrm{r}})=\mathsf{\Gamma}(\boldsymbol{\mathrm{r}},% \boldsymbol{\mathrm{r}})sansserif_Q ( bold_r ) = sansserif_Γ ( bold_r , bold_r ) and the Pauli matrices Λ1,2,3,±subscriptsans-serif-Λ123plus-or-minus\mathsf{\Lambda}_{1,2,3,\pm}sansserif_Λ start_POSTSUBSCRIPT 1 , 2 , 3 , ± end_POSTSUBSCRIPT. As explained in the companion paper [52], the saddle-point approximation is well justified in the nonlocalized regime (Lξmuch-less-than𝐿𝜉L\ll\xiitalic_L ≪ italic_ξ, where ξ𝜉\xiitalic_ξ is the localization length). Additionally, coupling between replicas can be neglected in this regime. This allows us to approximate ZRZ¯Rsimilar-to-or-equalsdelimited-⟨⟩superscript𝑍𝑅superscript¯𝑍𝑅\left\langle Z^{R}\right\rangle\simeq\bar{Z}^{R}⟨ italic_Z start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT ⟩ ≃ over¯ start_ARG italic_Z end_ARG start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT where Z¯¯𝑍\bar{Z}over¯ start_ARG italic_Z end_ARG is the partition function (10) corresponding to a single replica (R=1𝑅1R=1italic_R = 1). Since the problem is reduced to one replica, the field 𝖰(𝐫)𝖰𝐫\mathsf{Q}(\boldsymbol{\mathrm{r}})sansserif_Q ( bold_r ) becomes a 2×2222\times 22 × 2 matrix. After this simplification, it is relevant to approximate Eq. (11) semiclassically using the Wigner transform. In this regard, we define the matrix radiance 𝗀(𝛀,𝐫)𝗀𝛀𝐫\mathsf{g}(\boldsymbol{\mathrm{\Omega}},\boldsymbol{\mathrm{r}})sansserif_g ( bold_Ω , bold_r ) as

𝗀(𝛀,𝐫)=2kπkdpdd𝐬Γ(𝐫+𝐬2,𝐫𝐬2)eip𝛀𝐬.𝗀𝛀𝐫2𝑘𝜋subscript𝑘differential-d𝑝subscriptsuperscript𝑑differential-d𝐬sans-serif-Γ𝐫𝐬2𝐫𝐬2superscriptei𝑝𝛀𝐬\mathsf{g}(\boldsymbol{\mathrm{\Omega}},\boldsymbol{\mathrm{r}})=\frac{2k}{\pi% }\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$% \scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.291% 66pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\!\int_{k}\mathop{}% \!\mathrm{d}{p}\int_{\mathbb{R}^{d}}\mathop{}\!\mathrm{d}{\boldsymbol{\mathrm{% s}}}\>\mathsf{\Gamma}(\boldsymbol{\mathrm{r}}+\tfrac{\boldsymbol{\mathrm{s}}}{% 2},\boldsymbol{\mathrm{r}}-\tfrac{\boldsymbol{\mathrm{s}}}{2})\mathop{}\!% \mathrm{e}^{-\mathrm{i}p\boldsymbol{\mathrm{\Omega}}\cdot\boldsymbol{\mathrm{s% }}}\>.sansserif_g ( bold_Ω , bold_r ) = divide start_ARG 2 italic_k end_ARG start_ARG italic_π end_ARG - ∫ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT roman_d italic_p ∫ start_POSTSUBSCRIPT blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT end_POSTSUBSCRIPT roman_d bold_s sansserif_Γ ( bold_r + divide start_ARG bold_s end_ARG start_ARG 2 end_ARG , bold_r - divide start_ARG bold_s end_ARG start_ARG 2 end_ARG ) roman_e start_POSTSUPERSCRIPT - roman_i italic_p bold_Ω ⋅ bold_s end_POSTSUPERSCRIPT . (12)

The dashed integral on p𝑝pitalic_p in Eq. (12) only retains the Cauchy principal value on the wavenumber shell 𝐩=knorm𝐩𝑘\left\|\boldsymbol{\mathrm{p}}\right\|=k∥ bold_p ∥ = italic_k. As shown in Sec. III B of the companion paper [52], the result of the Wigner transform of Eq. (11) is a nonlinear matrix transport equation for 𝗀(𝛀,𝐫)𝗀𝛀𝐫\mathsf{g}(\boldsymbol{\mathrm{\Omega}},\boldsymbol{\mathrm{r}})sansserif_g ( bold_Ω , bold_r ),

𝛀𝐫𝗀=12s[𝖰~(𝐫),𝗀]ε2k[Λ3,𝗀]+iγΩxΘ(ma𝛀)δ(xxa)[Λ+,𝗀]+iΩxΘ(mb𝛀)δ(xxb)[Λ,𝗀],𝛀subscriptbold-∇𝐫𝗀12subscripts~𝖰𝐫𝗀𝜀2𝑘subscriptsans-serif-Λ3𝗀i𝛾subscriptΩ𝑥Θsubscript𝑚adelimited-∥∥subscript𝛀perpendicular-to𝛿𝑥subscript𝑥asubscriptsans-serif-Λ𝗀isubscriptΩ𝑥Θsubscript𝑚bdelimited-∥∥subscript𝛀perpendicular-to𝛿𝑥subscript𝑥bsubscriptsans-serif-Λ𝗀\begin{split}\boldsymbol{\mathrm{\Omega}}\cdot\boldsymbol{\mathrm{\nabla}}_{% \boldsymbol{\mathrm{r}}}\mathsf{g}=&-\tfrac{1}{2\ell_{\rm s}}[\tilde{\mathsf{Q% }}(\boldsymbol{\mathrm{r}}),\mathsf{g}]-\tfrac{\varepsilon}{2k}[\mathsf{% \Lambda}_{3},\mathsf{g}]\\ &+\mathrm{i}\gamma\Omega_{x}\Theta(m_{\mathrm{a}}-\left\|\boldsymbol{\mathrm{% \Omega}}_{\perp}\right\|)\delta(x-x_{\mathrm{a}})[\mathsf{\Lambda}_{+},\mathsf% {g}]\\ &+\mathrm{i}\Omega_{x}\Theta(m_{\mathrm{b}}-\left\|\boldsymbol{\mathrm{\Omega}% }_{\perp}\right\|)\delta(x-x_{\mathrm{b}})[\mathsf{\Lambda}_{-},\mathsf{g}]\>,% \end{split}start_ROW start_CELL bold_Ω ⋅ bold_∇ start_POSTSUBSCRIPT bold_r end_POSTSUBSCRIPT sansserif_g = end_CELL start_CELL - divide start_ARG 1 end_ARG start_ARG 2 roman_ℓ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT end_ARG [ over~ start_ARG sansserif_Q end_ARG ( bold_r ) , sansserif_g ] - divide start_ARG italic_ε end_ARG start_ARG 2 italic_k end_ARG [ sansserif_Λ start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT , sansserif_g ] end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL + roman_i italic_γ roman_Ω start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT roman_Θ ( italic_m start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT - ∥ bold_Ω start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT ∥ ) italic_δ ( italic_x - italic_x start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT ) [ sansserif_Λ start_POSTSUBSCRIPT + end_POSTSUBSCRIPT , sansserif_g ] end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL + roman_i roman_Ω start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT roman_Θ ( italic_m start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT - ∥ bold_Ω start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT ∥ ) italic_δ ( italic_x - italic_x start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT ) [ sansserif_Λ start_POSTSUBSCRIPT - end_POSTSUBSCRIPT , sansserif_g ] , end_CELL end_ROW (13)

reminiscent of the Eilenberger equation from the superconductivity literature [37, 38, *Larkin1973, *Larkin1975, *Larkin1977, 42, 56, 57, 34]. In Eq. (13), [𝖠,𝖡]=𝖠𝖡𝖡𝖠𝖠𝖡𝖠𝖡𝖡𝖠[\mathsf{A},\mathsf{B}]=\mathsf{A}\mathsf{B}-\mathsf{B}\mathsf{A}[ sansserif_A , sansserif_B ] = sansserif_AB - sansserif_BA is the commutator and 𝖰~(𝐫)~𝖰𝐫\tilde{\mathsf{Q}}(\boldsymbol{\mathrm{r}})over~ start_ARG sansserif_Q end_ARG ( bold_r ) is the normalized field:

𝖰~(𝐫)=def1πν𝖰(𝐫)=d𝛀Sd𝗀(𝛀,𝐫),~𝖰𝐫def1𝜋𝜈𝖰𝐫contour-integrald𝛀subscript𝑆𝑑𝗀𝛀𝐫\tilde{\mathsf{Q}}(\boldsymbol{\mathrm{r}})\overset{\mathrm{def}}{=}\frac{1}{% \pi\nu}\mathsf{Q}(\boldsymbol{\mathrm{r}})=\oint\frac{\mathop{}\!\mathrm{d}% \boldsymbol{\mathrm{\Omega}}}{S_{d}}\mathsf{g}(\boldsymbol{\mathrm{\Omega}},% \boldsymbol{\mathrm{r}})\>,over~ start_ARG sansserif_Q end_ARG ( bold_r ) overroman_def start_ARG = end_ARG divide start_ARG 1 end_ARG start_ARG italic_π italic_ν end_ARG sansserif_Q ( bold_r ) = ∮ divide start_ARG roman_d bold_Ω end_ARG start_ARG italic_S start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT end_ARG sansserif_g ( bold_Ω , bold_r ) , (14)

Sdsubscript𝑆𝑑S_{d}italic_S start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT being the surface area of the unit d𝑑ditalic_d-ball.

Equation (13) must be completed by appropriate boundary conditions that can be extracted from the asymptotic behavior of 𝗀(𝛀,𝐫)𝗀𝛀𝐫\mathsf{g}(\boldsymbol{\mathrm{\Omega}},\boldsymbol{\mathrm{r}})sansserif_g ( bold_Ω , bold_r ) prescribed by Eqs. (11)–(12) for 𝐫𝐫\boldsymbol{\mathrm{r}}\rightarrow\inftybold_r → ∞. We show in Sec. III C of the companion paper [52] that

𝗀out(𝛀,𝐫)=(1g1201),𝗀in(𝛀,𝐫)=(10g211).formulae-sequencesubscript𝗀out𝛀𝐫matrix1subscript𝑔1201subscript𝗀in𝛀𝐫matrix10subscript𝑔211\mathsf{g}_{\rm out}(\boldsymbol{\mathrm{\Omega}},\boldsymbol{\mathrm{r}})=% \begin{pmatrix}1&g_{12}\\ 0&-1\end{pmatrix}\>,\quad\mathsf{g}_{\rm in}(\boldsymbol{\mathrm{\Omega}},% \boldsymbol{\mathrm{r}})=\begin{pmatrix}1&0\\ g_{21}&-1\end{pmatrix}\>.sansserif_g start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT ( bold_Ω , bold_r ) = ( start_ARG start_ROW start_CELL 1 end_CELL start_CELL italic_g start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL - 1 end_CELL end_ROW end_ARG ) , sansserif_g start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ( bold_Ω , bold_r ) = ( start_ARG start_ROW start_CELL 1 end_CELL start_CELL 0 end_CELL end_ROW start_ROW start_CELL italic_g start_POSTSUBSCRIPT 21 end_POSTSUBSCRIPT end_CELL start_CELL - 1 end_CELL end_ROW end_ARG ) . (15)

The indices “in” and “out” denote the incoming and outgoing directions 𝛀𝛀\boldsymbol{\mathrm{\Omega}}bold_Ω on the disordered region, respectively. Note that Eq. (15) fixes only three elements of the matrix 𝗀𝗀\mathsf{g}sansserif_g out of four. The remaining elements g12subscript𝑔12g_{12}italic_g start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT and g21subscript𝑔21g_{21}italic_g start_POSTSUBSCRIPT 21 end_POSTSUBSCRIPT are free. As Eq. (13) satisfies the properties 𝐫[𝗀(𝛀,𝐫)2]=𝟬subscriptbold-∇𝐫𝗀superscript𝛀𝐫20\boldsymbol{\mathrm{\nabla}}_{\boldsymbol{\mathrm{r}}}\left[\mathsf{g}(% \boldsymbol{\mathrm{\Omega}},\boldsymbol{\mathrm{r}})^{2}\right]=\boldsymbol{% \mathrm{\mathsf{0}}}bold_∇ start_POSTSUBSCRIPT bold_r end_POSTSUBSCRIPT [ sansserif_g ( bold_Ω , bold_r ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ] = bold_sansserif_0 and 𝐫Tr𝗀(𝛀,𝐫)=𝟎subscriptbold-∇𝐫Tr𝗀𝛀𝐫0\boldsymbol{\mathrm{\nabla}}_{\boldsymbol{\mathrm{r}}}\operatorname{Tr}\mathsf% {g}(\boldsymbol{\mathrm{\Omega}},\boldsymbol{\mathrm{r}})=\boldsymbol{\mathrm{% 0}}bold_∇ start_POSTSUBSCRIPT bold_r end_POSTSUBSCRIPT roman_Tr sansserif_g ( bold_Ω , bold_r ) = bold_0, the boundary conditions (15) impose the following important constraints:

𝗀(𝛀,𝐫)2=𝟣,Tr𝗀(𝛀,𝐫)=0,𝛀,𝐫.formulae-sequence𝗀superscript𝛀𝐫21Tr𝗀𝛀𝐫0for-all𝛀𝐫\mathsf{g}(\boldsymbol{\mathrm{\Omega}},\boldsymbol{\mathrm{r}})^{2}=\mathsf{1% }\>,\quad\operatorname{Tr}\mathsf{g}(\boldsymbol{\mathrm{\Omega}},\boldsymbol{% \mathrm{r}})=0\>,\quad\forall\,\boldsymbol{\mathrm{\Omega}},\boldsymbol{% \mathrm{r}}\>.sansserif_g ( bold_Ω , bold_r ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = sansserif_1 , roman_Tr sansserif_g ( bold_Ω , bold_r ) = 0 , ∀ bold_Ω , bold_r . (16)

It is worth noting that the property 𝗀2=𝟣superscript𝗀21\mathsf{g}^{2}=\mathsf{1}sansserif_g start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = sansserif_1 is more general than the normalization 𝖰~2=𝟣superscript~𝖰21\tilde{\mathsf{Q}}^{2}=\mathsf{1}over~ start_ARG sansserif_Q end_ARG start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = sansserif_1 in nonlinear sigma models, which holds only in the isotropic limit due to the relationship given in Eq. (14). Once the system of equations (13)–(15) is solved, the generating function (4) can be calculated with:

F(γ)=i𝛀mad𝛀Ωx2Vd1Tr[𝗀(𝛀,xa)Λ+],𝐹𝛾isubscriptnormsubscript𝛀perpendicular-tosubscript𝑚ad𝛀subscriptΩ𝑥2subscript𝑉𝑑1Tr𝗀𝛀subscript𝑥asubscriptsans-serif-ΛF(\gamma)=\mathrm{i}\int_{\left\|\boldsymbol{\mathrm{\Omega}}_{\perp}\right\|% \leq m_{\mathrm{a}}}\frac{\mathop{}\!\mathrm{d}\boldsymbol{\mathrm{\Omega}}\>% \Omega_{x}}{2V_{d-1}}\operatorname{Tr}[\mathsf{g}(\boldsymbol{\mathrm{\Omega}}% ,x_{\mathrm{a}})\mathsf{\Lambda}_{+}]\>,italic_F ( italic_γ ) = roman_i ∫ start_POSTSUBSCRIPT ∥ bold_Ω start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT ∥ ≤ italic_m start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT end_POSTSUBSCRIPT divide start_ARG roman_d bold_Ω roman_Ω start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT end_ARG start_ARG 2 italic_V start_POSTSUBSCRIPT italic_d - 1 end_POSTSUBSCRIPT end_ARG roman_Tr [ sansserif_g ( bold_Ω , italic_x start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT ) sansserif_Λ start_POSTSUBSCRIPT + end_POSTSUBSCRIPT ] , (17)

Vdsubscript𝑉𝑑V_{d}italic_V start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT being the volume of the unit d𝑑ditalic_d-ball.

Now, we numerically investigate the validity of the theory (13)–(17) in a two-dimensional waveguide (d=2𝑑2d=2italic_d = 2) of transverse width W𝑊Witalic_W. We assume that the transverse boundary conditions are periodic [χn(y+W)=χn(y)subscript𝜒𝑛𝑦𝑊subscript𝜒𝑛𝑦\chi_{n}(y+W)=\chi_{n}(y)italic_χ start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ( italic_y + italic_W ) = italic_χ start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ( italic_y )]. This assumption aims at making the field 𝖰(x,y)𝖰𝑥𝑦\mathsf{Q}(x,y)sansserif_Q ( italic_x , italic_y ) independent of y𝑦yitalic_y in order to simplify the numerical procedures. Of course, the theory remains perfectly valid for Dirichlet boundary conditions if we let 𝖰(x,y)𝖰𝑥𝑦\mathsf{Q}(x,y)sansserif_Q ( italic_x , italic_y ) vary in y𝑦yitalic_y to vanish at the boundaries. We verified that the impact of boundary conditions is marginal in quasiballistic and diffusive regimes. The two equations (13) and (14) are solved iteratively until self-consistency is achieved. The integration of Eq. (13) along x𝑥xitalic_x under the constraints (16) can be carried out by means of the parametrization 𝗀(𝛀,𝐫)=𝖬(𝛀,𝐫)Λ3𝖬(𝛀,𝐫)1𝗀𝛀𝐫𝖬𝛀𝐫subscriptsans-serif-Λ3𝖬superscript𝛀𝐫1\mathsf{g}(\boldsymbol{\mathrm{\Omega}},\boldsymbol{\mathrm{r}})=\mathsf{M}(% \boldsymbol{\mathrm{\Omega}},\boldsymbol{\mathrm{r}})\mathsf{\Lambda}_{3}% \mathsf{M}(\boldsymbol{\mathrm{\Omega}},\boldsymbol{\mathrm{r}})^{-1}sansserif_g ( bold_Ω , bold_r ) = sansserif_M ( bold_Ω , bold_r ) sansserif_Λ start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT sansserif_M ( bold_Ω , bold_r ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT where 𝖬𝖬\mathsf{M}sansserif_M is an arbitrary matrix with two independent parameters. See the companion paper [52] for further details on the numerical integration of Eq. (13).

Our theoretical predictions are compared to transmission eigenvalue distributions computed numerically via the recursive Green’s function method, applied to the microscopic wave equation associated with Eq. (2), and averaged over disorder realizations.

Refer to caption
Figure 1: Transmission eigenvalue distribution for a disordered waveguide with W/λ=50.5𝑊𝜆50.5W/\lambda=50.5italic_W / italic_λ = 50.5, L/W=1𝐿𝑊1L/W=1italic_L / italic_W = 1, and ma=mb=1subscript𝑚asubscript𝑚b1m_{\mathrm{a}}=m_{\mathrm{b}}=1italic_m start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT = italic_m start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT = 1, at different optical thicknesses L/s𝐿subscriptsL/\ell_{\rm s}italic_L / roman_ℓ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT. Numerical distributions (dots) are computed from the microscopic wave equation and averaged over 5×1045superscript1045\times 10^{4}5 × 10 start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT disorder realizations, and compared to the predictions from the field theory (13)–(17) (solid) and the DMPK theory (dashed).

As can be seen in Fig. 1, the theory (13)–(17) (solid lines) agrees with numerical distributions (dots), including in the quasiballistic regime (Lsless-than-or-similar-to𝐿subscriptsL\lesssim\ell_{\rm s}italic_L ≲ roman_ℓ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT). The marked deviations of the DMPK theory [16, 17, 7] (dashed lines) in the region T<1𝑇1T<1italic_T < 1 arise solely from the approximate nature of the isotropy hypothesis. This hypothesis is only valid in diffusive regime (Lsmuch-greater-than𝐿subscriptsL\gg\ell_{\rm s}italic_L ≫ roman_ℓ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT) where the bimodal law applies (see Fig. 1 for L/s=5𝐿subscripts5L/\ell_{\rm s}=5italic_L / roman_ℓ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT = 5).

It is worth noting that the isotropy hypothesis can be implemented in the matrix transport equation (13) by setting

𝗀(𝛀,𝐫)=𝖰~(𝐫)+𝛀𝗝~(𝐫),𝗀𝛀𝐫~𝖰𝐫𝛀~𝗝𝐫\mathsf{g}(\boldsymbol{\mathrm{\Omega}},\boldsymbol{\mathrm{r}})=\tilde{% \mathsf{Q}}(\boldsymbol{\mathrm{r}})+\boldsymbol{\mathrm{\Omega}}\cdot\tilde{% \boldsymbol{\mathrm{\mathsf{J}}}}(\boldsymbol{\mathrm{r}})\>,sansserif_g ( bold_Ω , bold_r ) = over~ start_ARG sansserif_Q end_ARG ( bold_r ) + bold_Ω ⋅ over~ start_ARG bold_sansserif_J end_ARG ( bold_r ) , (18)

where 𝗝~(𝐫)~𝗝𝐫\tilde{\boldsymbol{\mathrm{\mathsf{J}}}}(\boldsymbol{\mathrm{r}})over~ start_ARG bold_sansserif_J end_ARG ( bold_r ) is a matrix current (a vector of matrices), and assuming that 𝗝~𝖰~much-less-thannorm~𝗝norm~𝖰\|\tilde{\boldsymbol{\mathrm{\mathsf{J}}}}\|\ll\|\tilde{\mathsf{Q}}\|∥ over~ start_ARG bold_sansserif_J end_ARG ∥ ≪ ∥ over~ start_ARG sansserif_Q end_ARG ∥, as done in the diffusion approximation. This approach leads to a Usadel-type equation [42, 33, 34], a nonlinear matrix diffusion equation closed for 𝖰~(𝐫)~𝖰𝐫\tilde{\mathsf{Q}}(\boldsymbol{\mathrm{r}})over~ start_ARG sansserif_Q end_ARG ( bold_r ). This equation can also be derived from the saddle-point equation of a nonlinear sigma model [14, 34, 57, 15] and matches the DMPK solution in the absence of absorption [43, 44, 45].

Additionally, Eq. (13) can be solved in the quasiballistic regime (Lsless-than-or-similar-to𝐿subscriptsL\lesssim\ell_{\rm s}italic_L ≲ roman_ℓ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT), where the field 𝖰~(𝐫)~𝖰𝐫\tilde{\mathsf{Q}}(\boldsymbol{\mathrm{r}})over~ start_ARG sansserif_Q end_ARG ( bold_r ) is nearly constant within the disordered region. This approximation serves as the dual to the isotropy hypothesis in DMPK theory: it approximates position space while treating direction space exactly, whereas DMPK does the reverse. As shown in Sec. III D of the companion paper [52], this approximation leads to a self-consistent equation for f(𝛀)𝑓𝛀f(\boldsymbol{\mathrm{\Omega}})italic_f ( bold_Ω ),

f(𝛀)=1(1f¯)1σtanh(Lσ2sΩx)1+(1+γ1γf¯)1σtanh(Lσ2sΩx),𝑓𝛀11¯𝑓1𝜎𝐿𝜎2subscriptssubscriptΩ𝑥11𝛾1𝛾¯𝑓1𝜎𝐿𝜎2subscriptssubscriptΩ𝑥f(\boldsymbol{\mathrm{\Omega}})=\frac{1-\left(1-\bar{f}\right)\frac{1}{\sigma}% \tanh(\frac{L\sigma}{2\ell_{\rm s}\Omega_{x}})}{1+\left(1+\frac{\gamma}{1-% \gamma}\bar{f}\right)\frac{1}{\sigma}\tanh(\frac{L\sigma}{2\ell_{\rm s}\Omega_% {x}})}\>,italic_f ( bold_Ω ) = divide start_ARG 1 - ( 1 - over¯ start_ARG italic_f end_ARG ) divide start_ARG 1 end_ARG start_ARG italic_σ end_ARG roman_tanh ( divide start_ARG italic_L italic_σ end_ARG start_ARG 2 roman_ℓ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT roman_Ω start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT end_ARG ) end_ARG start_ARG 1 + ( 1 + divide start_ARG italic_γ end_ARG start_ARG 1 - italic_γ end_ARG over¯ start_ARG italic_f end_ARG ) divide start_ARG 1 end_ARG start_ARG italic_σ end_ARG roman_tanh ( divide start_ARG italic_L italic_σ end_ARG start_ARG 2 roman_ℓ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT roman_Ω start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT end_ARG ) end_ARG , (19)

where σ𝜎\sigmaitalic_σ and the directional average f¯¯𝑓\bar{f}over¯ start_ARG italic_f end_ARG are defined by

σ=1+γf¯21γ,f¯=Ωx>0d𝛀12Sdf(𝛀).formulae-sequence𝜎1𝛾superscript¯𝑓21𝛾¯𝑓subscriptsubscriptΩ𝑥0d𝛀12subscript𝑆𝑑𝑓𝛀\sigma=\sqrt{1+\frac{\gamma\bar{f}^{2}}{1-\gamma}}\>,\quad\bar{f}=\int_{\Omega% _{x}>0}\frac{\mathop{}\!\mathrm{d}{\boldsymbol{\mathrm{\Omega}}}}{\frac{1}{2}S% _{d}}\>f(\boldsymbol{\mathrm{\Omega}})\>.italic_σ = square-root start_ARG 1 + divide start_ARG italic_γ over¯ start_ARG italic_f end_ARG start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 1 - italic_γ end_ARG end_ARG , over¯ start_ARG italic_f end_ARG = ∫ start_POSTSUBSCRIPT roman_Ω start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT > 0 end_POSTSUBSCRIPT divide start_ARG roman_d bold_Ω end_ARG start_ARG divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_S start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT end_ARG italic_f ( bold_Ω ) . (20)

The sought generating function is then given by inserting 𝗀(𝛀)=2i1γf(𝛀)Θ(Ωx)Λ𝗀𝛀2i1𝛾𝑓𝛀ΘsubscriptΩ𝑥subscriptsans-serif-Λ\mathsf{g}(\boldsymbol{\mathrm{\Omega}})=\tfrac{-2\mathrm{i}}{1-\gamma}f(% \boldsymbol{\mathrm{\Omega}})\Theta(\Omega_{x})\mathsf{\Lambda}_{-}sansserif_g ( bold_Ω ) = divide start_ARG - 2 roman_i end_ARG start_ARG 1 - italic_γ end_ARG italic_f ( bold_Ω ) roman_Θ ( roman_Ω start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT ) sansserif_Λ start_POSTSUBSCRIPT - end_POSTSUBSCRIPT into Eq. (17). The solution of Eqs. (19)–(20) is significantly more accurate than the DMPK theory in the quasiballistic regime (see Figs. 4–5 of the companion paper [52]).

Apart from its accuracy in the quasiballistic regime, another advantage of Eq. (13) is its ability to directly incorporate the effect of absorption in the medium by redefining the parameter ε𝜀\varepsilonitalic_ε as ε=k/a𝜀𝑘subscripta\varepsilon=k/\ell_{\rm a}italic_ε = italic_k / roman_ℓ start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT, where asubscripta\ell_{\rm a}roman_ℓ start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT denotes the ballistic absorption length. In Fig. 2 we compare the predictions of the full model (13) (solid lines) and its isotropic version based on Eq. (18) (dashed lines) with results from numerical simulations (dots).

Refer to caption
Figure 2: Transmission eigenvalue distribution for an absorbing disordered waveguide with W/λ=50.5𝑊𝜆50.5W/\lambda=50.5italic_W / italic_λ = 50.5, L/W=1𝐿𝑊1L/W=1italic_L / italic_W = 1, L/s=5𝐿subscripts5L/\ell_{\rm s}=5italic_L / roman_ℓ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT = 5, and ma=mb=1subscript𝑚asubscript𝑚b1m_{\mathrm{a}}=m_{\mathrm{b}}=1italic_m start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT = italic_m start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT = 1, shown for different absorption lengths asubscripta\ell_{\rm a}roman_ℓ start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT. Numerical distributions (dots) are compared with predictions from the anisotropic field theory (13)–(17) (solid) and its isotropic version based on Eq. (18) (dashed). The inset displays the distribution for L/s=L/a=0.2𝐿subscripts𝐿subscripta0.2L/\ell_{\rm s}=L/\ell_{\rm a}=0.2italic_L / roman_ℓ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT = italic_L / roman_ℓ start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT = 0.2, with the corresponding numerical distribution shown as dotted line.

Interestingly, we find that absorption tends to break the isotropy of 𝗀𝗀\mathsf{g}sansserif_g, even in the diffusive regime (L/s=5𝐿subscripts5L/\ell_{\rm s}=5italic_L / roman_ℓ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT = 5 in Fig. 2), owing to the variation in average path lengths taken by incoming waves based on their angle of incidence. This deviation disappears for L/s10much-greater-than𝐿subscripts10L/\ell_{\rm s}\gg 10italic_L / roman_ℓ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT ≫ 10. In the quasiballistic regime, lobes become visible (see inset of Fig. 2), resulting from higher transverse modes which, due to their extended paths in the absorbing medium, contribute to lower transmission. This lobe-splitting effect cannot be captured by an isotropic theory. We also emphasize that no theoretical prediction has been available until now for the experimentally relevant regime of moderate absorption, aLξasubscripta𝐿much-less-thansubscript𝜉a\ell_{\rm a}\leq L\ll\xi_{\rm a}roman_ℓ start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT ≤ italic_L ≪ italic_ξ start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT, where ξa=sa/2subscript𝜉asubscriptssubscripta2\xi_{\rm a}=\sqrt{\ell_{\rm s}\ell_{\rm a}/2}italic_ξ start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT = square-root start_ARG roman_ℓ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT roman_ℓ start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT / 2 end_ARG denotes the diffusive absorption length. The approximation proposed in Ref. [46] applies only in the strongly absorbing regime, Lξamuch-greater-than𝐿subscript𝜉aL\gg\xi_{\rm a}italic_L ≫ italic_ξ start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT, where transmission is minimal.

Finally, we consider an experimental situation where only modes within the input port’s numerical aperture masubscript𝑚am_{\mathrm{a}}italic_m start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT (i.e., 𝛀manormsubscript𝛀perpendicular-tosubscript𝑚a\left\|\boldsymbol{\mathrm{\Omega}}_{\perp}\right\|\leq m_{\mathrm{a}}∥ bold_Ω start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT ∥ ≤ italic_m start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT) are injected. At the outport port, we assume that all waveguide modes are accessible (mb=1subscript𝑚b1m_{\mathrm{b}}=1italic_m start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT = 1), making the transmission matrix 𝗍𝗍\mathsf{t}sansserif_t rectangular with more rows than columns.

Refer to caption
Figure 3: Transmission eigenvalue distribution for a disordered waveguide with W/λ=50.5𝑊𝜆50.5W/\lambda=50.5italic_W / italic_λ = 50.5, L/W=1𝐿𝑊1L/W=1italic_L / italic_W = 1, L/s=5𝐿subscripts5L/\ell_{\rm s}=5italic_L / roman_ℓ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT = 5, and L/a=0𝐿subscripta0L/\ell_{\rm a}=0italic_L / roman_ℓ start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT = 0, shown for different input port numerical apertures masubscript𝑚am_{\mathrm{a}}italic_m start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT. The output port has maximum numerical aperture (mb=1subscript𝑚b1m_{\mathrm{b}}=1italic_m start_POSTSUBSCRIPT roman_b end_POSTSUBSCRIPT = 1). Numerical distributions (dots) are compared with the predictions from the field theory (13)–(17) (solid) and FRM theory [49] with T¯=1/(1+2Lπs)0.239¯𝑇112𝐿𝜋subscriptssimilar-to-or-equals0.239\bar{T}=1/(1+\tfrac{2L}{\pi\ell_{\rm s}})\simeq 0.239over¯ start_ARG italic_T end_ARG = 1 / ( 1 + divide start_ARG 2 italic_L end_ARG start_ARG italic_π roman_ℓ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT end_ARG ) ≃ 0.239 (dashed).

The distributions obtained in this case are compared in Fig. 3 with the predictions of filtered random matrix (FRM) theory [49], in a regime of moderate optical thickness, where the FRM model is known to be imperfect. Once again, our new formalism provides a better fit to the observed distributions. In principle, Eq. (13) could also be extended to account for the effects of lateral diffusion in open geometries, avoiding the need for nontrivial renormalization of the FRM parameters in terms of long-range correlations [50, 51].

In conclusion, this Letter introduces an anisotropic field framework for the transmission eigenvalue distribution in disordered media. The core of this approach is the self-consistent matrix transport equation (13), which allows us to predict the correct distribution in the quasiballistic regime—a novel achievement, building on prior work [33, 57, 34] inspired by superconductivity [37, 38, *Larkin1973, *Larkin1975, *Larkin1977, 42]. Our model also rigorously accommodates medium absorption and practical measurement setups. Because of its flexibility, this framework opens up new avenues: it can describe incident beams with finite spatial extent [50, 51], statistics of observables beyond transmission, such as internal field intensity [47, 48], and is adaptable to non-rectilinear geometries [58]. Furthermore, due to its microscopic foundation, it can be extended to capture effects such as anisotropic scattering in biological tissues or the impact of correlated disorder [59]. We anticipate significant breakthroughs in understanding wave propagation in complex media through this approach.

Acknowledgements.
Acknowledgments—The authors thank Romain Pierrat for early-stage discussions and Romain Rescanières for valuable discussions and additional numerical simulations. This research has been supported by the ANR project MARS_light under reference ANR-19-CE30-0026 and by the program “Investissements d’Avenir” launched by the French Government.

References