Energy density distribution of shaped waves inside scattering media mapped onto a complete set of diffusion modes

We show that the spatial distribution of the energy density of optimally shaped waves inside a scattering medium can be described by considering only a few of the lowest eigenfunctions of the diffusion equation. Taking into account only the fundamental eigenfunction, the total internal energy inside the sample is underestimated by only 2%. The spatial distribution of the shaped energy density is very similar to the fundamental eigenfunction, up to a cosine distance of about 0.01. We obtained the energy density inside a quasi-1D disordered waveguide by numerical calculation of the joined scattering matrix. Computing the transmission-averaged energy density over all transmission channels yields the ensemble averaged energy density of shaped waves. From the averaged energy density obtained, we reconstruct its spatial distribution using the eigenfunctions of the diffusion equation. The results from our study have exciting applications in controlled biomedical imaging, efficient light harvesting in solar cells, enhanced energy conversion in solid-state lighting, and low threshold random lasers. © 2016 Optical Society of America OCIS codes: (000.0000) General. References and links 1. J. Goodman, Statistical Optics. Wiley, New York, 2000. 2. A. Ishimaru, Wave propagation and scattering in random media. Academic Press, New York, 1978. 3. M. C. W. van Rossum and T. M. Niewenhuizen, “Multiple scattering of classical waves,” Rev. Mod. Phys., vol. 71, pp. 313–371, 1999. 4. E. Akkermans and G. Montambaux, Mesoscopic physics of electrons and photons. Cambridge University Press, Cambridge, 2007. 5. I. Freund, “Looking through walls and around corners,” Physica A: Statistical Mechanics and its Applications, vol. 168, no. 1, pp. 49 – 65, 1990. 6. I. M. Vellekoop and A. P. Mosk, “Focusing coherent light through opaque strongly scattering media,” Opt. Lett., vol. 32, pp. 2309–2311, 2007. 7. I. M. Vellekoop and A. P. Mosk, “Universal optimal transmission of light through disordered materials,” Phys. Rev. Lett., vol. 101, p. 120601, 2008. 8. S. M. Popoff, A. Goetschy, S. F. Liew, A. D. Stone, and H. Cao, “Coherent control of total transmission of light through disordered media,” Phys. Rev. Lett., vol. 112, p. 133903, Apr 2014. 9. M. Davy, Z. Shi, and A. Z. Genack, “Focusing through random media: Eigenchannel participation number and intensity correlation,” Phys. Rev. B, vol. 85, p. 035105, Jan 2012. ar X iv :1 60 2. 03 92 1v 1 [ ph ys ic s. op tic s] 1 1 Fe b 20 16 10. A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, “Controlling waves in space and time for imaging and focusing in complex media,” Nature Photon., vol. 6, pp. 283–292, 2012. 11. I. M. Vellekoop, “Feedback-based wavefront shaping,” Opt. Express, vol. 23, pp. 12189–12206, May 2015. 12. A. Derode, P. Roux, and M. Fink, “Robust acoustic time reversal with high-order multiple scattering,” Phys. Rev. Lett., vol. 75, pp. 4206–4209, Dec 1995. 13. G. Lerosey, J. de Rosny, A. Tourin, A. Derode, G. Montaldo, and M. Fink, “Time reversal of electromagnetic waves,” Phys. Rev. Lett., vol. 92, no. 19, p. 193904, 2004. 14. G. Lerosey, J. de Rosny, A. Tourin, and M. Fink, “Focusing beyond the diffraction limit with far-field time reversal,” Science, vol. 315, pp. 1120–1122, 2007. 15. E. N. Leith and J. Upatnieks, “Holographic imagery through diffusing media,” J. Opt. Soc. Am., vol. 56, pp. 523– 523, Apr 1966. 16. D. R. Dowling and D. R. Jackson, “Narrow band performance of phase conjugate arrays in dynamic random media,” The Journal of the Acoustical Society of America, vol. 91, no. 6, 1992. 17. Z. Yaqoob, D. Psaltis, M. S. Feld, and C. Yang, “Optical phase conjugation for turbidity suppression in biological samples,” Nat. Photonics, vol. 2, pp. 110–115, 2008. 18. M. Kim, Y. Choi, C. Yoon, W. Choi, J. Kim, Q.-H. Park, and W. Choi, “Maximal energy transport through disordered media with the implementation of transmission eigenchannels,” Nat. Photonics, vol. 6, p. 581, 2012. 19. S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, “Measuring the transmission matrix in optics: An approach to the study and control of light propagation in disordered media,” Phys. Rev. Lett., vol. 104, p. 100601, 2010. 20. S. M. Popoff, G. Lerosey, M. Fink, A. C. Boccara, and S. Gigan, “Controlling light through optical disordered media: transmission matrix approach,” New Journal of Physics, vol. 13, no. 12, p. 123021, 2011. 21. Y. M. Wang, B. Judkewitz, C. A. DiMarzio, and C. Yang, “Deep-tissue focal fluorescence imaging with digitally time-reversed ultrasound-encoded light,” Nat Commun, vol. 3, p. 928, Jun 2012. 22. K. Si, R. Fiolka, and M. Cui, “Fluorescence imaging beyond the ballistic regime by ultrasound-pulse-guided digital phase conjugation,” Nat Photon, vol. 6, pp. 657–661, Oct 2012. 23. J. Bertolotti, E. G. van Putten, C. Blum, A. Lagendijk, W. L. Vos, and A. P. Mosk, “Non-invasive imaging through opaque scattering layers,” Nature, vol. 491, pp. 232–234, Nov 2012. 24. T. Cizmar, M. Mazilu, and K. Dholakia, “In situ wavefront correction and its application to micromanipulation,” Nat Photon, vol. 4, pp. 388–394, Jun 2010. 25. J.-H. Park, C. Park, H. Yu, Y.-H. Cho, and Y. Park, “Dynamic active wave plate using random nanoparticles,” Opt. Express, vol. 20, pp. 17010–17016, Jul 2012. 26. Y. Guan, O. Katz, E. Small, J. Zhou, and Y. Silberberg, “Polarization control of multiply scattered light through random media by wavefront shaping,” Opt. Lett., vol. 37, pp. 4663–4665, Nov 2012. 27. J.-H. Park, C. Park, H. Yu, Y.-H. Cho, and Y. Park, “Active spectral filtering through turbid media,” Opt. Lett., vol. 37, pp. 3261–3263, Aug 2012. 28. E. Small, O. Katz, Y. Guan, and Y. Silberberg, “Spectral control of broadband light through random media by wavefront shaping,” Opt. Lett., vol. 37, pp. 3429–3431, Aug 2012. 29. S. R. Huisman, T. Huisman, T. A. W. Wolterink, A. P. Mosk, and P. Pinkse, “Programmable multiport optical circuits in opaque scattering materials,” Opt. Express, vol. 23, p. 3102–3116, 02/2015 2015. 30. R. Horstmeyer, B. Judkewitz, I. M. Vellekoop, S. Assawaworrarit, and C. Yang, “Physical key-protected one-time pad,” Scientific Reports, vol. 3, pp. 3543 EP –, Dec 2013. Article. 31. S. A. Goorden, M. Horstmann, A. P. Mosk, B. Skoric, and P. Pinkse, “Quantum-secure authentication of a physical unclonable key,” Optica, vol. 1, pp. 421–424, 2014. 32. M. Krames, O. Shchekin, R. Mueller-Mach, G. O. Mueller, L. Zhou, G. Harbers, and M. Craford, “Status and future of high-power light-emitting diodes for solid-state lighting,” Display Technology, Journal of, vol. 3, pp. 160– 175, June 2007. 33. J. Phillips, M. Coltrin, M. Crawford, A. Fischer, M. Krames, R. Mueller-Mach, G. Mueller, Y. Ohno, L. Rohwer, J. Simmons, and J. Tsao, “Research challenges to ultra-efficient inorganic solid-state lighting,” Laser ‘&’ Photonics Reviews, vol. 1, no. 4, pp. 307–333, 2007. 34. V. Y. F. Leung, A. Lagendijk, T. W. Tukker, A. P. Mosk, W. L. IJzerman, and W. L. Vos, “Interplay between multiple scattering, emission, and absorption of light in the phosphor of a white light-emitting diode,” Opt. Express, vol. 22, pp. 8190–8204, 04/2014 2014. 35. T. Ogi, A. B. D. Nandiyanto, K. Okino, F. Iskandar, W.-N. Wang, E. Tanabe, and K. Okuyama, “Towards better phosphor design: Effect of sio2 nanoparticles on photoluminescence enhancement of yag:ce,” ECS Journal of Solid State Science and Technology, vol. 2, no. 5, pp. R91–R95, 2013. 36. M. Meretska, A. Lagendijk, H. Thyrrestrup, A. P. Mosk, W. L. IJzerman, and W. L. Vos, “How to distinguish elastically scattered light from stokes shifted light for solid-state lighting?,” arXiv:1511.00467, 2016. 37. J. A. Levitt and W. H. Weber, “Materials for luminescent greenhouse solar collectors,” Appl. Opt., vol. 16, pp. 2684–2689, Oct 1977. 38. A. Polman and H. Atwater Nat. Mater., vol. 11, p. 174, 2012. 39. F. T. Si, D. Y. Kim, R. Santbergen, H. Tan, R. A. C. M. M. van Swaaij, A. H. M. Smets, O. Isabella, and M. Zeman, “Quadruple-junction thin-film silicon-based solar cells with high open-circuit voltage,” Applied Physics Letters, vol. 105, no. 6, 2014. 40. O. Yizhar, L. E. Fenno, T. J. Davidson, M. Mogri, and K. Deisseroth, “Optogenetics in neural systems,” Neuron, vol. 71, no. 1, pp. 9 – 34, 2011. 41. W. Choi, A. P. Mosk, Q.-H. Park, and W. Choi, “Transmission eigenchannels in a disordered medium,” Phys. Rev. B, vol. 83, p. 134207, Apr 2011. 42. M. Davy, Z. Shi, J. Park, C. Tian, and A. Z. Genack Nat. Communications., vol. 6, p. 6893, 2015. 43. S. F. Liew and H. Cao, “Modification of light transmission channels by inhomogeneous absorption in random media,” Opt. Express, vol. 23, pp. 11043–11053, May 2015. 44. B. Gérardin, J. Laurent, A. Derode, C. Prada, and A. Aubry, “Full transmission and reflection of waves propagating through a maze of disorder,” Phys. Rev. Lett., vol. 113, p. 173901, Oct 2014. 45. O. S. Ojambati, H. Yılmaz, A. Lagendijk, A. P. Mosk, and W. L. Vos, “Selective coupling of optical energy into the fundamental diffusion mode of a scattering medium,” arXiv:1505.08103v2. 46. D. Y. K. Ko and J. C. Inkson, “Matrix method for tunneling in heterostructures: Resonant tunneling in multilayer systems,” Phys. Rev. B, vol. 38, pp. 9945–9951, Nov 1988. 47. J. B. Pendry and A. MacKinnon, “Calculation of photon dispersion relations,” Phys. Rev. Lett., vol. 69, pp. 2772– 2775, Nov 1992. 48. C. W. J. Beenakker, “Random-matrix theory of quantum transport,” Rev. Mod. Phys., vol. 69, pp. 731–808, 1997. 49. P. M. Morse and H. Feshbach, Methods of theoretical physics. New York, McGraw-Hill., 1953. 50. A. Lagendijk, R. Vreeker, and P. de Vries, “Influence of internal reflection on diffusive transport in strongly scattering media,” Phys. Lett. A, vol. 136, no. 1–2, pp. 81–88, 1989.


Introduction
There are numerous scattering media that are either of natural origin -such as biological tissue, fog, cloud, and teeth -or man-made -such as paint, diffuser glass, and phosphor inclusion in white LEDs. These materials multiply scatter waves in all directions, due to the random spatial variations of the refractive indices. Random interference occurs between the multiply scattered waves that is well-known as speckles [1]. Upon averaging over realizations of scatterers, the speckles average out and the resulting average energy density is well described by diffusion theory [2][3][4].
Although these methods are useful tools in controlling the intensity at the interfaces of a scattering medium, the energy density distribution of the shaped waves inside the scattering medium is still unknown. In the case of light, the knowledge of the energy density distribution is important for applications such as enhanced energy conversion in white LEDs [32][33][34][35][36], efficient light harvesting in solar cells [37][38][39], and controlled illumination in biomedical imaging [40]. To date, only numerical calculations of scalar waves [41][42][43] and a single-realization elastic wave experiment [44] have addressed the energy density distribution of shaped waves inside twodimensional (2D) scattering media. However, none of these studies provide an analytical model for the energy density distribution. For high-transmission channels, Davy et al [42] described the energy density using a parabolic function. In Ref. [45], our team reported the first measurement of the total energy density of light inside a three-dimensional (3D) scattering medium and reported an enhanced total energy density for shaped waves. The results were interpreted with a heuristic model that the total energy density of shaped waves can be described with only the fundamental eigenfunction of the diffusion equation, which agreed well with experimental observation. However, the exact distribution of the spatial energy density for shaped waves still remains unknown.
In this paper, we study the energy density of shaped waves inside a scattering medium by using the eigenfunctions of the diffusion equation to reconstruct the energy density distribution. Here, we refer to shaped waves to mean perfect phase conjugation of the transmitted waves, which is equivalent to shaping to optimally focus light to a single speckle spot. The diffusion eigenfunctions form a complete set and therefore describe the ensemble average energy density within the domain of the sample boundary. By numerical calculation of concatenated scattering matrices [46][47][48], we obtain the intensity transmission coefficients τ n of nth transmission channel and its energy density W n (z, τ) along the depth z of a quasi-1D disordered waveguide. A quasi-1D disordered waveguide provides a tractable platform to study the essential physics. By calculating the transmission-averaged energy density for all transmission channels, we obtained the energy density of shaped waves [7,11]. We show that only the first seven diffusion eigenfunctions are sufficient to reconstruct the spatial energy density distribution of the shaped waves. Taking into account only the fundamental eigenfunction of the diffusion equation, the total internal energy inside the sample is underestimated by only 2%. The spatial energy density is very similar to the fundamental eigenfunction, with a cosine distance of 0.01. Furthermore, we are able to reconstruct the energy density distribution of both high-and low-transmission channels, using a few M eigenfunctions of the diffusion equation, e.g. M = 16 for channel with transmission τ = 0.1. The reconstruction of the distribution of the energy densities that we found here can be extended to three dimensional (3D) scattering samples as well as to other forms of classical and quantum waves.

Eigenfunctions of the diffusion equation
The diffusion equation describes the ensemble averaged energy density W (r,t) of multiply scattered waves as a function of position r and time t inside a scattering medium [2,3], where D is the diffusion constant. Anticipating the fact that the directions (x, y) are decoupled from the z-direction we decompose the energy density as Since we are considering a slab that is a quasi-1D scattering medium, we assume translational invariance in the perpendicular directions (x, y). The sample boundaries are at z = 0 and z = L, where L is the sample thickness. This symmetry allows us to solve W ⊥ (x, y,t) by using its two-dimensional spatial Fourier transform W ⊥ (q ⊥ ,t). Solving the (x, y)-part of Eq. (1) gives Next we solve the remaining 1D equation We will turn this into a Sturm-Liouville eigenvalue problem by introducing the Ansatz W (z,t) = e −γt W (z) The values were obtained by solving Eq. 8. The extrapolation lengths are z e1 = z e2 = z e = (π/4) (for a 2D sample [4]) and L = 12.1 .
As is well-known from linear algebra [49] solutions of Sturm-Liouville equations are orthogonal and form a complete set. Applying boundary condition will lead to a discrete set of eigenvalues Γ m . The general eigenfunction for eigenvalue Γ m can be cast in the form As outlined in detail in Ref. [50] the wavevector κ m and phase η m will be determined from solving Eq. 5 with the appropriate boundary conditions. The paramater A m takes care of the normalization. The boundary conditions are derived from the sample parameters: the mean free path , the left and right extrapolation lengths z e1 , z e2 and the sample thickness L.
The allowed wave vectors κ m are obtained by solving the implicit equation Following Ref. [50], the phase factor η m fulfills the following equation and the normalization factor A m is given by The relation between eigenvalue and wavevector is Γ m = Dκ 2 m . For the exact solution of Eq. 8, see Appendix. The first ten wavevectors obtained from the solution are shown in Table 1 for a particular set of boundary conditions. In Fig. 1, we plot the first three eigenfunctions using Eq. 7 with the same set of boundary conditions.

Reconstructing the energy density with eigenfunctions of the diffusion equation
The diffusion equation is of the Sturm-Liouville type so its eigenfunctions form a complete, orthogonal, set [49]. Therefore, any distribution of energy density within the domain of the sample boundaries can be decomposed into a sum over a finite number M of eigenfunctions for a specific set of coefficients. Such a decomposition is only useful if the number of needed modes is small. We express the ensemble averaged energy density W (z) inside the scattering medium in terms of the normalized eigenfunctions W m (z) where d m are coefficients. We calculate the overlap integral I p of the product of W (z) and W p (z) and substitute Eq. 11 into Eq. 12 to obtain We use orthonormality of functions W p (z) and W m (z), which is expressed as where δ mp is the Kronecker delta. Therefore, Eq. 13 becomes We now can write the decomposition (11) as We obtained W (z) for shaped waves and different transmission eigenchannels, which we obtained from the simulation described in Section 3. The simulation results are compared to the reconstructed energy density W re (z) in Section 4. In order to quantify the overlap between the reconstructed function W re (z) and the numerical data W (z), we use cosine distance COSD Ref[], which is defined as N s is the number of points in the numerical data. COSD varies between 1 and 0 and tends towards 0 as the reconstructed function fully describes the numerical data. As a figure-of-merit of a good reconstruction, we choose the eigenfunction, which has COSD ≈ 10 −4 .

Numerical setup
We perform simulations of transport of monochromatic scalar waves through a waveguide containing scatterers. The waveguide is modeled as a quasi-one-dimensional (quasi-1D) system possessing N m = 100 transversal modes, with elastic scattering on bulk impurities. This is implemented as a 2D waveguide with 100 transversal modes, containing scatterers at random positions with a density of 0.2 scatterer per λ 2 , λ is the wavelength. The distribution of the scattering potential P of scatterers in one waveguide configuration is shown in Fig. 2. In order to calculate the energy density inside the waveguide we employ a recursive S-matrix formalism, where the waveguide is first divided into sub-wavelength slices, each containing only few scatterers. For each slice we define the S-matrix in a mode representation as follows [48], where R L and R R are the left and right reflection matrix respectively of dimension N m × N m , T is the left-to-right transmission matrix. The right-to-left transmission matrix in this reciprocal system is the transpose of T . The S-matrix of each slice, of dimensions 2N m ×2N m , is calculated in the approximation that no recurrent scattering takes place within the thin slice. Next, the S-matrices of the slices are joined, including contributions of recurrent scattering, using the composition rule for S-matrices [46]. This numerically stable procedure yields the S-matrix of the entire waveguide and of subsections of it. From the transmission matrix we extract the singular values and vectors corresponding to all channels by singular value decomposition. The energy density inside the waveguide, for a given incident field, is calculated as follows: For a given axial coordinate z c the S matrices S 1 ,S 2 of the waveguide sections 0 < z < z c and z c < z < L respectively are calculated. The field in the z = z c plane is then found by calculating recurrent scattering diagrams: Here, R L2 is the left reflection submatrix of S 2 , R R1 is the right reflection submatrix of S 1 , and T 1 is the left to right transmission submatrix of S 1 . Subsequently for every z s the time-averaged energy density ε 0 |E c | 2 is integrated over the cross section to obtain the projected energy density W (z). The calculation is repeated for 8000 independently generated random waveguides from which the transmission channels are extracted. To obtain plots of the average energy density as a function of transmission, we average the energy densities of channels in a narrow band of transmission values.

Sample characterization
In order to characterize the numerical samples, we plot in Fig. 3a the probability distribution P(τ) as a function of the transmission τ. The probability distribution P(τ) obtained from the simulation is bi-modal: there is a high probability for transmitting channels with transmission τ close to zero and one. In Fig. 3a, we compare the Dorokhov-Mello-Pereyra-Kumar (DMPK) distribution [48] with our numerical result. The probability distribution of transmission for a scattering medium is expected to follow the DMPK distribution. The DMPK distribution agrees well with our numerical result. Furthermore, in order to confirm that our numerical samples are in the diffusive regime, we plot in Fig. 3b the equally-weighted ensemble averaged energy density versus reduced sample depth. The ensemble averaged energy density shows a linear decrease from the front surface towards the end surface of the sample, in agreement with the prediction of diffusion theory for a diffusive sample, see Fig. 3b.

Energy density of shaped waves
We obtained the energy density of shaped waves W sw as the transmission-weighted ensemble average of the energy density of transmission channels W n (z, τ), as shown in Ref. [7,11]: where N is the number of transmission channels. Using Eq. 21, we obtained the energy density of the shaped waves from the transmission-weighted averaged over the energy density of all transmission channels. We show in Fig. 4a the ensemble averaged energy density distribution of shaped waves. The energy density distribution of the shaped waves increases from the interfaces of the sample and is maximum at z = 4.9 , which is about 20% off the center of the sample (z = 6.1 ). The ensemble averaged energy density of wavefront-shaped light was first obtain in Ref. [41] by solving the Maxwell equation using the finite-difference time-domain (FDTD) method. Despite the different numerical calculation methods, our numerical result is closely similar to the one obtained by Choi et al [41]. The exact peak position of the energy density obtain in [41] is difficult to estimate due to the noise in the data, however, the peak position is close to the center of the sample. We reconstruct the distribution of the energy density of shaped waves from eigenfunctions of the diffusion equation. From the reconstruction, we obtain the cosine distance COSD which we plot in Fig. 4. In case of shaped waves, a summation over the first M = 7 eigenfunctions is remarkably sufficient (COSD ≈ 10 −4 ) to reconstruct the energy density, see Fig. 4a. Only the diffusion fundamental eigensolution m = 1 diffusion eigensolution has COSD ≈ 10 −2 and is not sufficient to describe the energy density profile. The fundamental eigensolution peaks at exactly the center of the sample, in contrast to the energy density of shaped waves.
We further investigated the contribution of the eigenfunctions to the total energy density integrated along the depth of the sample relative to the contribution of the first 100 eigenfunctions. As a reference, we choose the first 100 eigenfunctions because the higher order eigenfunctions has a very negligible contribution. In Fig. 5, we plot the contributions. For shaped waves, the m = 1 contributes 98% of the total energy density inside the sample. Remarkably, our finding about the contribution of m = 1 to the total internal energy density agrees very well with the heuristic model in Ref. [45] that wavefront shaping predominately couples light to the fundamental diffusion eigensolution. The contribution of m = 2 is negligible compared to m = 1 since m = 2 has negative energy density in almost half of the sample depth, which cancels out the positive energy density (see Fig. 1). A similar effect also happens for other even index eigenfunction.

Energy density of transmission channels
The ensemble averaged energy density of open transmitting eigenchannels with a transmission τ in the range 0.99 < τ < 1 is shown in Fig. 6. The energy density shows a symmetric profile, which peaks in the middle of the sample. Our numerical result is similar to the one obtained by Davy et al [42]. Interestingly, only the fundamental diffusion eigenfunction m = 1 is sufficient to reconstruct the energy density of the open channel, see Fig. 4a. The fundamental diffusion eigenfunction m = 1 is expected to describe the open eigenchannels for the following reasons: Firstly, as the open eigenchannel has the highest individual transmission, the fundamental eigenfunction (m = 1) as well contributes most to the total transmission as shown in Ref. [45]. Secondly, the fundamental eigenfunction is the only physical solution with a positive energy density along the sample depth z, see Fig. 1. Interestingly, a simplified mathematical model by Davy et al [42] described the energy density of open channels as a parabolic function. We note that the deviations from a parabolic shape are obvious only near the edges and the shape of the energy density curve agrees better with the cosine shape of the fundamental eigenfunction, see Fig. 6(a). The COSD for the parabola is 10 −3 and 10 −5 for the fundamental eigenfunction.
In order to describe the distribution of energy density of low-transmission channels, we also use the eigenfunctions of the diffusion equation to reconstruct the ensemble averaged energy densities of low-transmission channels obtained from simulation. In Figs. 6 (b -d), we show the ensemble averaged energy densities for transmission eigenchannels with transmissions τ = 0.5, 0.10 and 0.01. For the eigenchannel with transmission τ = 0.5, only M = 7 eigenfunctions of the diffusion equation are sufficient to reconstruct its energy density, M = 16 eigenfunctions for τ = 0.1, and M = 33 is required for a closed channel τ = 0.01. As the transmission τ decreases, the number of eigenfunctions sufficient to reconstruct the energy density increases. The increase in the number M of sufficient eigenfunctions is because the asymmetry of the distribution of energy density increases as the transmission τ reduces, and likewise, the asymmetry of the eigenfunctions increases with the eigensolution index. Therefore, the higher index eigenfunctions contribute more to an asymmetric function. The contribution of the eigenfunctions with index m = 1, 2, 3 and 4 relative to a summation over the first 100 eigenfunctions is shown in Fig. 7. The relative contribution of the fundamental eigenfunction of the diffusion equation m = 1, which is a symmetric function increases as the transmission increases, while the relative contributions of the higher index eigenfunctions, which are asymmetric functions are higher for low-transmission channels.

Effect of sample thickness
Here, we investigate the effect of the sample thickness on the number of eigenfunctions sufficient to reconstruct the energy density of shaped waves. We therefore perform simulation of the energy density on different samples with various sample thicknesses: L = 5 , 12.1 , 18.8 and 39.5 . We show in Fig. 8 the COSD versus the eigensolution index for the 4 samples. Interesting, all the samples have a convergence to a value of COSD ≈ 10 −4 at the seventh eigensolution. This convergence of the COSD shows that the summation over the first 7 eigenfunctions of the diffusion equation is sufficient to describe the energy density of shaped waves, irrespective of the sample thickness. There is however a slight deviation for the thickest sample L = 39.5 and this deviation is probably due to artefacts from the simulation due to the large number of waveguide slices.
In the previous sections, we have shown that the diffusion fundamental eigenfunction of the diffusion equation m = 1 dominates the energy density of high-transmission channels. We investigate if this result holds for the different samples thicknesses. In Fig. 8, we show the coefficient of the m = 1 mode relative to the sum of the first 100 eigenfunctions versus the transmission for the 4 different thicknesses. The relative contribution of the m = 1 mode is similar for all samples, irrespective of the thickness and this is quite remarkable. The diffusion fundamental eigensolution dominates the energy density of high-transmission channels for all sample thicknesses.

Conclusion
We have shown that only a few eigenfunctions of the diffusion equation suffice to accurately reconstruct the distribution of the shaped energy density inside a quasi-1D scattering waveguides. In particular, the fundamental eigenfunction of the diffusion equation is very similar to the distribution of energy density of shaped waves. To reconstruct the distribution of open channels only the diffusion fundamental eigenfunction m = 1 is sufficient. In addition, we have shown that a few number of diffusion eigenfunctions reconstruct the energy density of low transmit-ting channels. Our results are relevant for applications that require the precise knowledge of distribution of the energy density inside scattering media. Such applications include as efficient light harvesting in solar cells especially in near infrared where silicon has low absorption; for enhanced energy conversion in white LEDs, which serves to reduce the quantity of expensive phosphor; for low threshold and higher output yield of random lasers; as well as in homogeneous excitation of probes in biological tissues.

Acknowledgment
We thank Pepijn Pinkse, Henri Thyrrestrup, and Hasan Yilmaz for useful discussions. This project is part of the research program of the "Stichting voor Fundamenteel Onderzoek der Materie" (FOM) FOM-program "Stirring of light!", which is part of the "Nederlandse Organisatie voor Wetenschappelijk Onderzoek" (NWO). We acknowledge NWO-Vici, DARPA, and STW.

A. Method of calculating the wave vectors of diffusion eigenfunctions
In Eq. 8, we show an implicit equation, which defines the allowed wave vectors of the diffusion eigenfunctions. The equation is again stated here: We define where and Our goal is to find the wave vectors κ r m that fulfills the condition In Fig. 9, we plot the two functions f 1 and f 2 . The condition in Eq. 26 is fulfilled by the wave vectors κ r m at the intersection of the two functions f 1 and f 2 , as indicated in Fig. 9. We found the wave vectors κ r m by dividing function f into domains with size π/L. In most of the domains, only one root is expected to be in each of the domains, expect the domain where function f 2 has a divergence, which has two roots. In each domain, we use a Simple Bisection method to search the domain for wave vectors that fulfill the desired condition. Each domain is divided into two equal parts. The function f is first evaluated for left part of the domain and then check if there is any root present. If there is no root present, the left domain is further divided into smaller subdomains and each subdomains is searched for the root. We divide the domains until the size of the subdomain is 10 −12 . If there is no root in all the subdomains of the left domain, then, the right domain is searched. This way, all the desired roots are found, iteratively.   Eigenfunction index m Cum. relative contrib.