Universal structure of transmission eigenchannels inside opaque media

As the desire to explore opaque materials is ordinarily frustrated by multiple scattering of waves, attention has focused on the transmission matrix of the wave field. This matrix gives the fullest account of transmission and conductance and enables the control of the transmitted flux; however, it cannot address the fundamental issue of the spatial profile of eigenchannels of the transmission matrix inside the sample. Here we obtain a universal expression for the average disposition of energy of transmission eigenchannels within random diffusive systems in terms of auxiliary localization lengths determined by the corresponding transmission eigenvalues. The spatial profile of each eigenchannel is shown to be a solution of a generalized diffusion equation. These results reveal the rich structure of transmission eigenchannels and enable the control of the energy distribution inside random media.

The electronic conductance in units of the quantum of conductance is equivalent to the optical transmittance, which is the sum over all flux transmission coefficients, T ¼ P N a;b t ba j j 2 , where the t ba are elements of the TM between N incident and outgoing channels, a and b, respectively 23 . These channels may be the transverse modes of empty waveguides or momentum states of ideal leads connecting to a disordered sample at a given frequency or energy. The eigenchannels or natural channels of transmission are linear combinations of these channels, which can be obtained from the singular value decomposition of the TM, t ¼ P N n¼1 u n ffiffiffiffi ffi t n p v þ n . Here u n and v n are unit vectors comprising the n th transmission eigenchannel and t n are the corresponding transmission eigenvalues 23 . The transmittance may be expressed as the sum over the N transmission eigenvalues, T ¼ P N n¼1 t n . Moreover, many key statistical properties of transport through random media such as the fluctuations and correlations of conductance and transmission may be described in terms of the statistics of the transmission eigenvalues [2][3][4][5][6][7][8][9][10] .
Notwithstanding the success of the TM in describing the statistics of transmission through disordered systems, it cannot shed light on the disposition of energy inside the sample. Recent simulation [24][25][26][27][28] and measurements 22 in single realizations of random samples suggest that the peak of the energy density profile moves towards the centre of the sample and the total energy of the corresponding eigenchannel increases as the transmission eigenvalue t increases. This is consistent with the measurements of the composite phase derivative of transmission eigenchannels, which show that the integrated energy inside the sample increases with t (ref. 29). The possible universality of the structure of each transmission eigenchannels within scattering media remains an unexplored question, which is the subject of this work. In contrast, the average energy distribution over all N eigenchannels can be obtained by solving the diffusion equation yielding a linear falloff of energy inside a diffusive sample governed by Fick's law for particle diffusion. For localized waves, the energy distribution inside the sample is governed by a generalized diffusion equation and deviates dramatically from a linear fall. This deviation can be explained in terms of a positiondependent diffusion coefficient [30][31][32][33][34][35][36][37] .
We seek to discover a universal expression for the average over a random ensemble of samples of the spatial distribution of the energy density averaged over the transverse dimensions for eigenchannels with specified transmission eigenvalue t. This would extend our understanding of the universality in wave propagation from the boundaries to the interior of the sample. Aside from its fundamental importance, a universal description of the scaling of energy density in eigenchannel allows for the control of the energy density profile within the sample. This could be exploited to obtain depth profiles of random media in measurements that are sensitive to optical absorption, emission or nonlinearity. The possibility of depositing energy well below the surface would lengthen the residence time of the emitted photons in active random systems and thus lower the lasing threshold of amplifying diffusive media below that for traditionally pumped random lasers. The pump threshold for random lasers is traditionally high because the residence times of emitted photons are relatively short as a consequence of the shallow penetration of the pump laser due to multiple scattering [38][39][40] .
We show that the average of the energy density profile within the sample of an individual eigenchannel with transmission t is closely related to the generalized diffusion equation 30-37. The energy density of the perfectly transmitting eigenchannel with t ¼ 1 is essentially the sum of a spatially uniform background equal to the energy density of the incident and outgoing wave and the product of this background energy density and the probability density for the wave to return to the cross section at given position within the sample. For to1, we find that the energy density can be expressed as a product of the profile of the fully transmitting eigenchannel and a function governed only by the auxiliary localization length, which was previously used to parameterize the corresponding transmission eigenvalue 1,2 .

Results and Discussion
Energy density profiles. The energy density profiles averaged over 500 configurations drawn from a random ensemble with N ¼ 66 and L/x ¼ 0.05 for transmission eigenvalues t ¼ 1, 0.5, 0.1 and 0.001 are shown in linear and semilog plots in Fig. 1a,b. The localization length x is determined from oT4 ¼ px/2L (refs 8,9), where oy4 indicates the average over an ensemble of random samples. For ease of presentation, we normalize the energy density W t (x) so that it is equal to the transmission coefficient t at the output surface, W t (L) ¼ t. This is accomplished by dividing the energy density by the average longitudinal component of the wave velocity at the output surface, v þ . In general, W t (x) may be associated with the sum of the forward and backward flux divided by v þ . The normalized energy density at the incident surface is The linear falloff of the average of the energy density over all eigenchannels is in accord with the diffusion theory. (d) W t (x) for a sample with half the width and length as the sample in a. Note that, though value of L/x is the same as for a, the form of W t (x) differs with lower values of the peaks in the shorter sample. equal to 2 À t. For the Lambertian probability distribution for waves transmitted through random two-dimensional (2D) media without internal reflection at the boundaries, v þ ¼ pc/4, where c is the speed of the wave in the medium. A single peak in W t (x) is seen to increase and shift towards the middle of the sample as t increases. The average of the energy density profiles over all the transmission eigenchannels, WðxÞ ¼ N À 1 P N n¼1 W tn ðxÞ, is seen in Fig. 1c to fall linearly as expected 25 . Simulations for a shorter and wider sample with L/x ¼ 0.03 are shown in Fig. 1d. The peak values of the energy density are lower in this sample for all values of t. The details of the simulations are given in the first part of the Methods section.
Perfectly transmitting eigenchannel. We first consider the energy density profile within the sample of the perfectly transmitting eigenchannel. This is seen in Fig. 1 to be of the form, with F 1 (x) a symmetric function peaked in the middle of the sample, which vanishes at the longitudinal boundaries of the sample. This is reminiscent of the probability density for a wave to return to a cross section at depth x in an open random medium. This probability density in the case where the surface internal reflection vanishes can be obtained from the correlation function, Y(x, x 0 ), which is RR dydy 0 o Gðx; y; x 0 ; y 0 Þ j j 2 4 up to an inessential overall factor that depends upon the wave frequency and the average density of states 35 is the ensemble average of local energy density averaged over the cross section at x due to a unit plane source located at x 0 . By using exact microscopic methods, it was found 31 propagating in open random media. We emphasize that, as in conventional diffusion, D(x) is an intrinsic observable independent of the sources and is determined by microscopic parameters of the system such as disorder strength and incident wavelength. The explicit expression of D(x) can be obtained by either exact microscopic theories [31][32][33] or approximate self-consistent treatments 30 .
With D(x) given explicitly by either analytical or numerical methods, the generalized diffusion equation leads to a generalization of Fick's law giving the average flux, À DðxÞ d dx WðxÞ, and gives Y(x, x 0 ) once the appropriate boundary conditions are implemented. We solve the generalized diffusion equation with the boundary conditions, being the generalized probability density of return to the cross section at x (see the second part of the Methods section).
In the diffusive limit, D(x) reduces to the Boltzmann diffusion coefficient D 0 and equation (2) is the return probability density in the diffusive limit [30][31][32][33] . In 2D, the diffusion coefficient is D 0 ¼ cc/2, where c is the transport mean free path. With the energy density rescaled by v þ , we obtain Y 1 (x,x) ¼ px(L À x)/2Lc. Plots of F 1 (x) for three diffusive samples are shown in Fig. 2a. When normalized to their peak values at L/2 in Fig. 2b, these profiles collapse to a single curve, 4x(L À x)/L 2 , showing that F 1 (x) equals Y 1 (x,x) for diffusive waves. The peak value F 1 (L/2) is seen in Fig. 2c to increase linearly with L/c with a slope found from the expression above for F 1 (x), namely F 1 (L/2) ¼ pL/8c.
The relationship between F 1 (x) and Y 1 (x, x) is further explored by comparing these functions for localized waves for which Y 1 (x, x) depends upon L/x. The profiles of F 1 (x/L)/F 1 (L/2) obtained in simulations for samples with N ¼ 10 and in onedimensional (1D) samples (the first part of the Methods section) are seen in Fig. 3a to narrow as L/x increases. Y 1 (x,x) is computed from equation (2) with D(x) determined from the gradient of W(x) found in simulations as shown in Fig. 3b for the same samples as in Fig. 3a. Profiles of Y 1 (x, x)/Y 1 (L/2, L/2) and F 1 (x)/ F 1 (L/2) are seen in Fig. 3a to match for each of the random ensembles studied, thereby establishing the result, x). In Supplementary Note 1 and Supplementary Figures 1 and 2 we show that this result has a counterpart even in single random samples.
Factorization of the energy density profile. We find from simulations for diffusive waves that the energy density profile can be written as a product in which S t (x) is independent of the value of L/x and L/c, as demonstrated in Fig. 4a. This factorization is proved for both diffusive and localized waves by a non-perturbative diagrammatic technique, and the diagrammatic meaning of S t (x) in terms of the interactions caused by dielectric fluctuations within the medium is given in Supplementary Note 2 and Supplementary Figures 3  and 4. In the case of the perfectly transmitting eigenchannel, x) consists of two parts. One is a spatially uniform background energy density. The other, and the boundary conditions to Y( Solving this equation (see the second part of the Methods section) and setting , which is the same as equation (3). Thus, we establish the relation between the factorization and the generalized diffusion equation. We have not found the analytical form for S t (x) for arbitrary to1. Instead, we could conjecture the form by considering the distribution of transmission eigenvalues. Dorokhov 1,2 showed that for quasi 1D scattering media the variation of average conductance with sample length L (ref. 41) could be understood in terms of the correlated scaling of the transmission eigenvalues in terms of a set of auxiliary localization lengths x n , which gives the average transmission eigenvalue of the n th eigenchannel as cosh À 2 (L/x n ). For noN/2, x À 1 n ¼ ð2n À 1Þ=ð2xÞ varies linearly with channel index n (ref. 8). For diffusive waves, this corresponds to BoT4 open eigenchannels with the average transmission eigenvalue 41/e, 1-3 while for localized waves, this implies that the transmission is dominated by a single eigenchannel and the probability distribution of L/x 0 exhibits peaks at L/x n (refs 42,43). Here the auxiliary localization length x 0 is associated with the transmission eigenvalue t ¼ cosh À 2 (L/x 0 ). To find an expression for S t (x), we hypothesize that average intensity profile in each eigenchannel is related to x 0 .
For diffusive waves, S t (x) is governed only by L/x 0 . A natural assumption is that the analytic continuation of S t (x) inside the sample at the boundaries is S t (L) ¼ t and S t (0) ¼ 2 À t so that W t (L) ¼ t and W t (0) ¼ 2 À t. Thus, S t (x) at the boundaries must reduce to S t (L) ¼ cosh À 2 (L/x 0 ) and S t (0) ¼ 2 À cosh À 2 (L/x 0 ). Perhaps the simplest expression consistent with these conditions is S t ðxÞ ¼ 2 cosh 2 ðL À xÞ=x 0 ð Þ cosh 2 L=x 0 ð Þ À t. However, this expression does not agree with the results of simulations shown in Fig. 1. Rather, agreement with simulations is only found once we incorporate an empirical function h(x/L) into the expression above, giving To find the function h(x/L), we solve equation (5) with S t (x) found in simulations for a given value of t. Surprisingly, we find that h(x/L) obtained in this way for a single value of t in one diffusive sample gives good agreement for all t in all diffusive samples (Fig. 4a). h(x/L) is plotted in Fig. 4b. W t (x) found from equations (1,3,5) using this form for h(x/L) are plotted as dashed black curves in Fig. 1a,b,d. These are in excellent agreement with the profiles found in simulations for diffusive waves. The universal form of energy density profiles of transmission eigenchannels found here apply to all classical and quantum waves in homogeneously disordered samples. The auxiliary localization lengths for transmission eigenvalues proposed by Dorokhov are seen to determine the average energy density profiles of the transmission eigenchannels for samples of specified length and transmission for random ensembles and so provide approximate profiles in individual samples. Aside from extending the knowledge of the average disposition of energy from the boundary into the bulk, the expression for W t (x) provides a window on eigenchannel dynamics and the density of states. The integral of W t (x) is the delay time for the transmission eigenvalue t and gives the contribution of the eigenchannel to the density of states of the sample 29 . The sum of the integrals of W t (x) over all channels is the density of states of the medium. For waves such as light and sound for which it is possible to control the incident waveform, these results make it possible to tailor the distribution of excitation within a random medium. This may allow depth profiling of random media and make it possible to deliver radiation deep into multiply scattered media.

Methods
Numerical simulations. The profiles of the energy density are found from recursive Green's function simulations [44][45][46] . We consider the propagation of a scalar wave in a quasi 1D strip, which is locally 2D, with length L much greater than the transverse dimension L t . The precise geometry of the cross section, which is linear in our study, does not influence the results since waves are mixed in the transverse direction within the bulk of the quasi 1D sample [1][2][3][4][5][6][7][8][9][10] . The results are expected to remain valid for locally 3D samples. We consider a model of the disordered region stripped of inessential nonuniversal elements. The random sample is index matched to its surroundings with index of refraction of unity so that both the internal and external reflections are negligible. Fluctuations de(x, y) in the position-dependent dielectric constant, e(x, y) ¼ 1 þ de(x, y), are drawn from a rectangular distribution, the width of which determines the strength of disorder and so the average transmittance of the sample. The wave equation r 2 Eðx; yÞ þ k 2 0 eðx; yÞEðx; yÞ ¼ 0 is discretized on a square grid and solved using the recursive Green's function method. The product of the wave number in the leads k 0 and the grid spacing is set to unity.
The Green's function G(r,r 0 ) is calculated between grid points r ¼ (0,y) and r 0 ¼ (x 0 ,y 0 ) on the incident plane x ¼ 0 and in the interior of the sample at depth x 0 , with y being the transverse coordinate. The field transmission coefficient at x 0 between the incoming mode a and outgoing mode b is calculated by projecting the Green's function onto the empty waveguide modes f n (y), where v a is the group velocity of the empty waveguide mode a at the frequency of the wave. The product of t(x 0 ) and v n yields the field at x 0 for different channels due to the n th incoming eigenchannel. In particular, summing the square of the field amplitude at x 0 ¼ L over all N channels gives the transmission coefficient t n . Averages are taken over 500 statistically equivalent samples for each ensemble so that precise comparisons can be made with potential forms for energy density profiles. For all the values of t in all quasi 1D samples, simulations of W t (x) are averaged over a subensemble of eigenchannels with transmission between 0.98 t and 1.02 t.
We also carry out scattering matrix simulations in a 1D disordered system to explore the profiles of W t ¼ 1 (x) and D(x). An electromagnetic plane wave is normally incident on a layered system with fixed number of layers L with alternating indices of refraction between a fixed value n and 1. The layer thickness is drawn from a distribution that is much greater than the wavelength. The mean free path in this system is equal to the localization length, c ¼ x (ref. 47). W t ¼ 1 (x) is computed from the average of the energy density over a subset of samples with transmission t ranging from 0.995 to 1. The boundary condition is correspondingly modified such that Y(z(0), z 0 ) ¼ S t (z 0 )/ v þ . This is the normal diffusion equation and is readily solved to give Yðz; z 0 Þ ¼ S t ðz 0 ÞÂ 1=v þ þ ðz À zð0ÞÞðzðLÞ À z 0 Þ=ðzðLÞ À zð0ÞÞ; zoz 0 ; 1=v þ þ ðz 0 À zð0ÞÞðzðLÞ À zÞ=ðzðLÞ À zð0ÞÞ; z4z 0 : : We note that since D(x)40, the function z(x) is strictly monotonic and is therefore invertible. Taking this into account and substituting the identity, zðx 2 Þ À zðx 1 Þ ¼ R Setting x ¼ x 0 , we find Yðx; xÞ ¼ v À 1 þ S t ðxÞð1 þ Y 1 ðx; xÞÞ with the generalized return probability density given by equation (2).