Introduction

The transmission matrix (TM), whose elements are field transmission coefficients through disordered media, provides a powerful basis for a complete understanding of the statistics of transmission1,2,3,4,5,6,7,8,9,10 and of the degree of control that may be exerted over the transmitted wave by manipulating the incident wavefront11,12,13,14,15,16,17,18,19,20,21,22. Control of the transmitted wave has been demonstrated and characterized in recent years in tight focusing11,12,13,14,15, and in enhanced or suppressed transmission16,17,18,19,20,21,22 of sound11, elastic waves22, light12,13,14,16,18,20,21 and microwave radiation15,19. The striking similarities in the statistics and scaling of transport of classical and electron waves in disordered systems reflect the equivalent descriptions of transmission and conductance in terms of the TM, t (refs 1, 2, 3, 4, 5, 6, 7, 8, 9, 10).

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, , where the tba are elements of the TM between N incident and outgoing channels, a and b, respectively23. 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, . Here un and vn are unit vectors comprising the nth transmission eigenchannel and τn are the corresponding transmission eigenvalues23. The transmittance may be expressed as the sum over the N transmission eigenvalues, . 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 eigenvalues2,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 simulation24,25,26,27,28 and measurements22 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 τ 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 τ (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 position-dependent diffusion coefficient30,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 τ. 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 scattering38,39,40.

We show that the average of the energy density profile within the sample of an individual eigenchannel with transmission τ is closely related to the generalized diffusion equation 30–37. The energy density of the perfectly transmitting eigenchannel with τ=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 τ<1, 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 eigenvalue1,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/ξ=0.05 for transmission eigenvalues τ=1, 0.5, 0.1 and 0.001 are shown in linear and semilog plots in Fig. 1a,b. The localization length ξ is determined from <T>=πξ/2L (refs 8, 9), where <…> indicates the average over an ensemble of random samples. For ease of presentation, we normalize the energy density Wτ(x) so that it is equal to the transmission coefficient τ at the output surface, Wτ(L)=τ. 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τ(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 equal to 2−τ. For the Lambertian probability distribution for waves transmitted through random two-dimensional (2D) media without internal reflection at the boundaries, v+=πc/4, where c is the speed of the wave in the medium. A single peak in Wτ(x) is seen to increase and shift towards the middle of the sample as τ increases. The average of the energy density profiles over all the transmission eigenchannels, , is seen in Fig. 1c to fall linearly as expected25. Simulations for a shorter and wider sample with L/ξ =0.03 are shown in Fig. 1d. The peak values of the energy density are lower in this sample for all values of τ. The details of the simulations are given in the first part of the Methods section.

Figure 1: Energy density profiles.
figure 1

(a) Ensemble averages of the eigenchannel energy density profiles Wτ(x) for eigenvalues τ=1, τ=0.5, τ=0.1 and τ=0.001 for a diffusive sample with L/ξ=0.05. The black dashed curves are plots of results given by equations (1, 3, 5). (b) Semilog plot of data in a showing results for small values of the energy density. (c) Ensemble averages of energy density profiles for all the transmission eigenchannels with the eigenchannel indices n from 1 to N, . The linear falloff of the average of the energy density over all eigenchannels is in accord with the diffusion theory. (d) Wτ(x) for a sample with half the width and length as the sample in a. Note that, though value of L/ξ is the same as for a, the form of Wτ(x) differs with lower values of the peaks in the shorter sample.

Full size image

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 F1(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′), which is up to an inessential overall factor that depends upon the wave frequency and the average density of states35, evaluated at x=x′. Physically, Y(x, x′) is the ensemble average of local energy density averaged over the cross section at x due to a unit plane source located at x′. By using exact microscopic methods, it was found31,32,33 that Y(x, x′) is the solution of the generalized diffusion equation, for waves 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 theories31,32,33 or approximate self-consistent treatments30.

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, , and gives Y(x, x′) once the appropriate boundary conditions are implemented. We solve the generalized diffusion equation with the boundary conditions, Y (x=0, x′)=Y(x=L, x′)=1/v+. Upon setting x=x′ in the solution, we find Y(x, x)=1/v+(1+Y1(x, x)), with

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 D0 and equation (2) yields Y1(x, x)/v+=x(L−x)/D0L. The quadratic form of Y1(x, x) is the return probability density in the diffusive limit30,31,32,33. In 2D, the diffusion coefficient is D0=cℓ/2, where ℓ is the transport mean free path. With the energy density rescaled by v+, we obtain Y1(x,x)=πx(L−x)/2Lℓ. Plots of F1(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)/L2, showing that F1(x) equals Y1(x,x) for diffusive waves. The peak value F1(L/2) is seen in Fig. 2c to increase linearly with L/ℓ with a slope found from the expression above for F1(x), namely F1(L/2)=πL/8ℓ.

Figure 2: Profiles for perfectly transmitting eigenchannels.
figure 2

(a) Ensemble averages of the perfectly transmitting channel F1(x) for the two diffusive samples of Fig. 1 with L/ξ=0.05 and a sample with L/ξ=0.025. (b) F1(x) normalized by its peak value F1(L/2). (c) Scaling of F1(L/2) with L/ℓ for samples with N=66 (blue triangles) and N=33 (black circles). The linear fit to the results of simulations with a slope of π/8=0.39 given by Y1(L/2,L/2) gives plotted as the red curve.

Full size image

The relationship between F1(x) and Y1(x, x) is further explored by comparing these functions for localized waves for which Y1(x, x) depends upon L/ξ. The profiles of F1(x/L)/F1(L/2) obtained in simulations for samples with N=10 and in one-dimensional (1D) samples (the first part of the Methods section) are seen in Fig. 3a to narrow as L/ξ increases. Y1(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 Y1(x, x)/Y1(L/2, L/2) and F1(x)/F1(L/2) are seen in Fig. 3a to match for each of the random ensembles studied, thereby establishing the result, F1(x)=Y1(x,x) and Wτ=1(x) =1+Y1(x, x). In Supplementary Note 1 and Supplementary Figures 1 and 2 we show that this result has a counterpart even in single random samples.

Figure 3: Perfectly transmitting eigenchannel profiles and position-dependent diffusion coefficients.
figure 3

(a) F1(x) normalized by F1(L/2) for multichannel samples (N=10) and for 1D samples (N=1) with values of L/ξ ranging from 0.25 to 11.7. The dashed black curves are plots of Y1(x,x) normalized by Y1(L/2,L/2) found from equation (2) with position-dependent diffusion coefficient D(x) shown in b. (b) Normalized position-dependent diffusion coefficient D(x)/D(0) found from the gradient of W(x) shown in the inset for the same samples as in a.

Full size image

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τ(x) is independent of the value of L/ξ and L/ℓ, 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τ(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, Wτ=1(x)=1+Y1(x, x) consists of two parts. One is a spatially uniform background energy density. The other, Y1(x, x′=x), is associated with the probability of return of waves that have reached the cross section at x, with the source of uniform strength, Sτ=1(x)=1, provided by the background energy density that enters into the right-hand side of the generalized diffusion equation as the coefficient of the Dirac delta function. For eigenchannels with τ<1, the strength of the source term Sτ(x) falls with increasing depth into the sample since waves do not readily penetrate into the bulk. This modifies the generalized diffusion equation to

Figure 4: Factorization of energy density profiles.
figure 4

(a) Wτ(x)/Wτ=1(x) for the three ensembles of samples of Fig. 1 (blue and red curves) and a sample with L/ξ=0.025 (green curve). The black dashed curves are given by equation (5) with h(x/L) shown in b. Curves of simulations of Wτ(x)/Wτ=1(x) for the same values of τ in different samples and for Sτ(x) given by equation (5) overlap. (b) h(x/L) determined from equation (5).

Full size image

and the boundary conditions to Y(x=0, x′)=Y(x=L, x′)=Sτ(x)/v+. Solving this equation (see the second part of the Methods section) and setting x=x′ in the solution gives Wτ(x)=Sτ(x)+Sτ(x)Y1(x, x), 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τ(x) for arbitrary τ<1. Instead, we could conjecture the form by considering the distribution of transmission eigenvalues. Dorokhov1,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 ξn, which gives the average transmission eigenvalue of the nth eigenchannel as cosh−2(L/ξn). For n<N/2, varies linearly with channel index n (ref. 8). For diffusive waves, this corresponds to ∼<T> open eigenchannels with the average transmission eigenvalue >1/e,1,2,3 while for localized waves, this implies that the transmission is dominated by a single eigenchannel and the probability distribution of L/ξ′ exhibits peaks at L/ξn (refs 42, 43). Here the auxiliary localization length ξ′ is associated with the transmission eigenvalue τ=cosh−2(L/ξ′). To find an expression for Sτ(x), we hypothesize that average intensity profile in each eigenchannel is related to ξ′.

For diffusive waves, Sτ(x) is governed only by L/ξ′. A natural assumption is that the analytic continuation of Sτ(x) inside the sample at the boundaries is Sτ(L)=τ and Sτ(0)=2−τ so that Wτ(L)=τ and Wτ(0)=2−τ. Thus, Sτ(x) at the boundaries must reduce to Sτ(L)=cosh−2(L/ξ′) and Sτ(0)=2−cosh−2(L/ξ′). Perhaps the simplest expression consistent with these conditions is . 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τ(x) found in simulations for a given value of τ. Surprisingly, we find that h(x/L) obtained in this way for a single value of τ in one diffusive sample gives good agreement for all τ in all diffusive samples (Fig. 4a). h(x/L) is plotted in Fig. 4b. Wτ(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τ(x) provides a window on eigenchannel dynamics and the density of states. The integral of Wτ(x) is the delay time for the transmission eigenvalue τ and gives the contribution of the eigenchannel to the density of states of the sample29. The sum of the integrals of Wτ(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 simulations44,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 Lt. 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 sample1,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 δɛ(x, y) in the position-dependent dielectric constant, ɛ(x, y)=1+δɛ(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 is discretized on a square grid and solved using the recursive Green’s function method. The product of the wave number in the leads k0 and the grid spacing is set to unity.

The Green’s function G(r,r′) is calculated between grid points r=(0,y) and r′=(x′,y′) on the incident plane x=0 and in the interior of the sample at depth x′, with y being the transverse coordinate. The field transmission coefficient at x′ between the incoming mode a and outgoing mode b is calculated by projecting the Green’s function onto the empty waveguide modes φn(y), , where va is the group velocity of the empty waveguide mode a at the frequency of the wave. The product of t(x′) and vn yields the field at x′ for different channels due to the nth incoming eigenchannel. In particular, summing the square of the field amplitude at x′=L over all N channels gives the transmission coefficient τ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 τ in all quasi 1D samples, simulations of Wτ(x) are averaged over a subensemble of eigenchannels with transmission between 0.98 τ and 1.02 τ.

We also carry out scattering matrix simulations in a 1D disordered system to explore the profiles of Wτ=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, ℓ=ξ (ref. 47). Wτ=1(x) is computed from the average of the energy density over a subset of samples with transmission τ ranging from 0.995 to 1.

Solution of the generalized diffusion equation

Let us multiply both sides of equation (4) by D(x). With the change of variable, which gives dz=dx/D(x), we rewrite equation (4) as

The boundary condition is correspondingly modified such that Y(z(0), z′)=Sτ(z′)/v+. This is the normal diffusion equation and is readily solved to give

We note that since D(x)>0, the function z(x) is strictly monotonic and is therefore invertible. Taking this into account and substituting the identity, , into the expression of Y(z, z′), we may rewrite Y(z, z′) as

Setting x=x′, we find with the generalized return probability density given by equation (2).

Additional information

How to cite this article: Davy, M. et al. Universal structure of transmission eigenchannels inside opaque media. Nat. Commun. 6:6893 doi: 10.1038/ncomms7893 (2015).