Fluctuations and correlations of transmission eigenchannels in diffusive media
FFluctuations and correlations of transmission eigenchannels in diffusive media
Nicholas Bender , ∗ Alexey Yamilov , ∗ , † Hasan Yılmaz , and Hui Cao ‡ Department of Applied Physics, Yale University, New Haven, Connecticut 06520, USAPhysics Department, Missouri University of Science & Technology, Rolla, Missouri 65409, USA (Dated: April 28, 2020)Selective excitation of a diffusive systems transmission eigenchannels enables manipulation of itsinternal energy distribution. The fluctuations and correlations of the eigenchannels spatial profiles,however, remain unexplored so far. Here we show that the depth profiles of high-transmissioneigenchannels exhibit low realization-to-realization fluctuations. Furthermore, our experimental andnumerical studies reveal the existence of inter-channel correlations, which are significant for low-transmission eigenchannels. Because high-transmission eigenchannels are robust and independentfrom other eigenchannels, they can reliably deliver energy deep inside turbid media.
In recent years, extensive studies of coherent wavetransport in multiple-scattering media have been con-ducted with light, microwaves, and acoustic waves [1, 2].One of the goals of this research is overcoming the limita-tions imposed by incoherent diffusion: thereby enablingenergy delivery deep inside a turbid medium. In any lin-ear scattering system with static disorder, the coherentwave transport is described by the field transmission ma-trix, t , which maps the incident waves to the transmittedwaves [3]. Performing a singular value decomposition onthe transmission matrix gives t = U √ τ V † ; where thecolumns of V and U are the input and output wavefrontsof the transmission eigenchannels, and the diagonal ele-ments of τ are the transmission eigenvalues [4–8].Transmission eigenchannels form a unique basis for t because they diagonalize the matrix. Each eigenchannel, V α , independently propagates through the system de-scribed by t with a transmittance given by its correspond-ing eigenvalue τ α . Any incoming waves can be decom-posed into a linear superposition of eigenchannels. Be-cause transmission eigenchannels are specified by the dis-order configuration of a system, they vary from system-to-system. It is unknown whether eigenchannels varyin a coordinated manner between statistically equivalentdisordered systems, and if so, how?Transmission eigenchannels are more than basis func-tions. In addition to having different transmittances [9–15], each eigenchannel also has a distinct energy distri-bution inside diffusive systems [16–24]. For example, aneigenchannel with τ (cid:39) ∗ These two authors contributed equally. † [email protected] ‡ [email protected] a single sample where the eigenchannel profile deviatesfrom the ensemble-averaged eigenchannel profile. Know-ing how significant this deviation is, from realization-to-realization, is essential in many practical applications.Here, we experimentally and numerically investigateboth the fluctuations and correlations of transmissioneigenchannel depth profiles in optical diffusive sys-tems. We fabricate two-dimensional disordered waveg-uides, excite their individual eigenchannels, and mea-sure the eigenchannel field profiles within the waveg-uides. The high-transmission eigenchannels exhibit smallrealization-to-realization fluctuations in their depth pro-files, demonstrating a robustness when compared to ei-ther low-transmission eigenchannels or random inputs.Furthermore, different eigenchannels are correlated intheir fluctuations from realization-to-realization. Thecorrelations are weaker for higher-transmission eigen-channels, indicating they are more independent thanlower-transmission eigenchannels. Our results uncoverthe second-order statistical properties of transmissioneigenchannels, which are promising for applications indeep tissue imaging and light delivery.To directly study the depth profiles of transmissioneigenchannels within a diffusive system, we fabricatetwo-dimensional (2D) waveguide structures on a silicon-on-insulator wafer with electron beam lithography andplasma etching [33]. As shown in Fig. 1(a), 100 nm-diameter holes are randomly etched into the waveguides,which have photonic crystal sidewalls to reflect light [34].At the wavelength of our probe light, λ = 1 . µ m,the transport mean free path, (cid:96) t = 3 . µ m, is muchshorter than the disordered region length, L = 50 µ m, ineach waveguide. Therefore, the light undergoes multiplescattering and diffusive transport through each waveg-uide [33]. Light scatters out-of-plane from the randomholes, providing a direct probe of the light inside thedisordered region. This process can be modeled as aneffective loss, and accounted for in the diffusive dissi-pation length: ξ a = 28 µ m. The waveguides are each15 µ m wide, supporting N = 55 propagating modes at λ = 1 . µ m. Before entering one of the diffusive waveg-uides, light is injected via the edge of the wafer into aridge waveguide. Due to the large refractive index mis-match between silicon and air, only low-order waveguide a r X i v : . [ c ond - m a t . d i s - nn ] A p r 𝑚𝑎𝑥 𝒕 slm→end 𝒕𝒕 slm→buff 𝒕 buff→end = 𝒕 slm→end ∙ 𝒕 slm→buff−1 (a) −𝜋 𝜋 (c) Measured Intensity PatternReconstructed Phase Pattern
End (b)
Buffer Region
Waveguide Structure
FIG. 1. Waveguide structure and full-field measurement. Acomposite SEM image of a diffusive waveguide is shown in(a). The matrix mapping the field in the buffer region to theend region, t buff → end , is related to the matrices t slm → buff and t slm → end . In (b) the 2D intensity pattern of a measured high-transmission eigenchannel is shown. Using our interferometricsetup, we can reconstruct the phase of the light field (c). In (b-c) the edges of the diffusive region are marked by the verticaldashed lines. modes are excited at the interface. Before the disorderedregion, the waveguide width is tapered from 300 µ m to 15 µ m in order to convert the lower-order modes to higher-order ones. The taper enables us to access all waveguidemodes incident on the disordered region [22].To measure the light field inside individual diffusivewaveguides, we use an interferometric setup describedin [33]. In our setup, the monochromatic light from awavelength-tunable laser source is split into two beams.One beam is modulated by a spatial light modulator(SLM) and then injected into one of the waveguides viathe edge of the wafer. The other beam is used as a ref-erence beam. It is spatially overlapped with the out-of-plane scattered light from the diffusive waveguide: onthe CCD camera chip. The CCD camera records theresulting interference pattern, from which the complexfield profile across the diffusive waveguide is obtained, asshown in Fig. 1 (b,c).By sequentially applying an orthogonal set of phasepatterns to the 128 SLM macropixels, and measuring thefield within the sample, we acquire a matrix that mapsthe field from the SLM to the field inside the disorderedwaveguide t slm → int . This operator encompasses informa-tion about the light transport inside the waveguide andthe light propagation from the SLM to the waveguide. Toseparate these, we need access to the field incident on thedisordered region of the waveguide. We obtain this infor-mation by adding an auxiliary weakly-scattering regionin front of the diffusive region called the ‘buffer’ region:as shown in Fig. 1 (a). From the light scattered out-of-
22 eigenchannel profiles
High-transmission channel
Low-transmission channel (a) (b)(c)
Depth z ( m m) I n t en s i t y I ( z ) α I n t en s i t y I ( z ) n t en s i t y I ( z ) ExperimentSimulationSimulation
Depth z ( m m)Depth z ( m m) FIG. 2. Depth profiles of transmission eigenchanels. High( α = 1) and low ( α = 20) transmission eigenchannel profilesare presented in (a,b) while the 22 measured eigenchannelprofiles are juxtaposed in (c). The experimentally measuredprofiles (blue lines) agree well with the profiles calculated fromnumerical simulations using the transmission matrix t (blackdashed lines) and the matrix t buff → end (red lines). plane from the buffer, we recover the field right in front ofthe strongly-scattering region. The length of the bufferregion is 25 µ m, which is shorter than its 32 µ m-lengthtransport mean free path. Therefore, light only experi-ences single scattering in the buffer, and as a result, thediffusive wave transport in the original disordered regionis not appreciably altered.With access to the field inside the buffer, we can con-struct the matrix relating the field on the SLM to thebuffer, t slm → buff . From t slm → int , we can also constructthe matrix, t slm → end , which maps the field from the SLMto a region near the end of the diffusive waveguide. Withthese we calculate the matrix which maps the field fromthe buffer to the end, t buff → end = t slm → end t − → buff , us-ing Moore-Penrose matrix inversion. Although t buff → end is not the field transmission matrix t , the depth profilesof its eigenchannels match those of transmission eigen-channels in our numerical simulation (see Fig. 2 and dis-cussion below). Therefore, t buff → end can be used as anexperimental proxy for the field transmission matrix, t ,of the diffusive waveguide.To excite a single eigenchannel, we first perform a sin-gular value decomposition on t buff → end to obtain the fielddistribution in the buffer corresponding to one eigenchan-nel. Then we multiply the field profile in the bufferwith t − → buff to calculate the requisite SLM phase-modulation pattern. By displaying this pattern on theSLM, we excite a single eigenchannel of the diffusivewaveguide. We record the spatial intensity profile ofeach eigenchannel within the diffusive waveguide. Fromthis measurement, we obtain the eigenchannel depth pro-file ˜ I ( z ) associated with each measurement by summingthe intensity over the waveguide cross-section. For eachdepth profile, the measured intensity profile ˜ I ( z ) is nor-malized to I ( z ) = ˜ I ( z ) / [(1 /L ) (cid:82) L ˜ I ( z (cid:48) ) dz (cid:48) ].In Figures 2(a) and (b), the experimentally-measureddepth profiles of a high-transmission ( α = 1) and a low-transmission ( α = 20) eigenchannel are juxtaposed. Thehigh-transmission eigenchannel in (a) has an arch-shapedenergy-density distribution which spans the depth ofthe diffusive region. In (b), the energy-density distri-bution of the low-transmission eigenchannel rapidly de-cays with depth. The experimentally measured profilesmatch the corresponding depth profiles generated fromnumerical simulations of both t and t buff → end ; confirmingthat we excite individual eigenchannels in our measure-ments. Furthermore, the agreement between the eigen-channels of t buff → end and t , confirms that the depth pro-files of t buff → end have a one-to-one correspondence withthe eigenchannels of t .In total, we measure 50 eigenchannel profiles fora single experimental system realization. Each pro-file matches one of the ensemble-averaged profiles of t buff → end generated numerically without any fitting pa-rameters [33]. The experimental profiles are indexed by α E and the numerical profiles by α . Measurement noisecauses multiple α E to be mapped to one α , and this limits α ≤
22. Fig. 2(c) shows the depth profiles for all 22 re-covered eigenchannels, which agree well with the numer-ical simulations [33]. The transmittance of the measuredeigenchannels varies from τ (cid:39) .
43 to τ (cid:39) . × − ,with a mean value of (cid:104) τ α (cid:105) = 0 . (cid:104) I α ( z ) (cid:105) , and the realization-specific deviation, δI α ( z ) = I α ( z ) − (cid:104) I α ( z ) (cid:105) . From this,the total fluctuation of each eigenchannel profile is quan-tified by ˜ C α = (1 /L ) (cid:82) L (cid:104) [ δI α ( z )] (cid:105) dz , where (cid:104) ... (cid:105) repre-sents ensemble averaging. Fig. 3(a) shows that the totalfluctuation of each eigenchannel profile increases mono-tonically as a function of eigenchannel index. The un-certainty of ˜ C α -due to the finite number of ensemblesin our experiment- is estimated from simulations to be ±
25% the value of ˜ C α : which is smaller than the over-all change of ˜ C α with α . Hence, the depth profiles ofhigh-transmission eigenchannels fluctuate less than theprofiles generated by random illumination patterns (indi-cated by the green dashed line); while lower-transmissioneigenchannels fluctuate more.Now we look into the position-dependent fluctuationsof individual eigenchannel profiles about their ensem-ble average, var[ I α ( z )] = (cid:104) [ δI α ( z )] (cid:105) as a function of Eigenchannel index
Experiment
Depth z ( m m) E i gen c hanne l i nde x -3 -2 -1 Simulation
0 10 20 30 40 505101520 E i gen c hanne l i nde x Depth z ( m m) -3 -2 -1 Depth z ( m m) -3 -2 -1
10 0 10 20 30 40 50 (a) (b) (e) (d)
Eigenchannel fluctuationsDepth-resolved fluctuations α = 1 α = 20 Experiment
Simulation
Simulation
Simulation
0 10 20 30 40 505101520 E i gen c hanne l i nde x Depth z ( m m) -3 -2 -1 (c) (f) Depth z ( m m) Relative variance
Random Inputs C α ~ v a r αα v a r α α = 1 α = 20 FIG. 3. Eigenchannel fluctuations. In (a), the spatially-averaged depth-profile fluctuations of the eigenchannels, ˜ C α ,increase monotonically with the channel index α . The greendashed line indicates the experimentally observed fluctuationsfor random incident wavefronts: 0.59. In (b), experimentallyobserved depth-resolved intensity fluctuations, var[ I α ( z )], ofhigh ( α = 1) and low ( α = 20) transmission eigenchan-nels (circles) are closely reproduced by the numerical sim-ulations of transmission eigenchannels from t buff → end (solidlines) and t (dashed lines). In (c), var[ I α ( z )] is dividedby (cid:104) I α ( z ) (cid:105) for the high/low-transmission eigenchannels of t . In (d-e), the experimentally-observed and numerically-calculated depth-resolved intensity fluctuations for individualeigenchannels show how var[ I α ( z )] evolves with α . depth, z . Fig. 3(b) reveals distinct differences in thedepth dependence of high and low-transmission eigen-channels. While var[ I α ( z )] is nearly flat for the high-transmission eigenchannel, it features a fast drop with z for the low-transmission eigenchannel. Figs. 3(d-f) are2D plots of var[ I α ( z )] for all 22 eigenchannels: calcu-lated using experimental data, as well as simulations of t buff → end , and t . As the transmittance decreases, themaximum of var[ I α ( z )] moves towards the front surfaceof the diffusive region. The decrease in the variancewith depth results from the decay of the mean inten-sity with depth: (cid:104) I α ( z ) (cid:105) . However, the relative intensityfluctuations of the low-transmission eigenchannels, char-acterized by var[ I α ( z )] / (cid:104) I α ( z ) (cid:105) , actually increase withdepth as shown in Fig. 3(c) for α = 20. In contrast, therelative intensity fluctuation of high-transmission eigen- -3 -2 -1 V a r i a n c e C o v a r i a n c e -3 -2 -1 -3 -2 -1 (a)c)( (b)(d) Variance and covarianceSimulationExperiment Simulation E i gen c hanne l i nde x Eigenchannel index Eigenchannel index Eigenchannel index Eigenchannel index E i gen c hanne l i nde x E i gen c hanne l i nde x FIG. 4. Inter-channel correlations. The covariance ˜ C αβ between any two pairs of eigenchannels, α and β , are calcu-lated from experimental data (a) and numerical simulations(b,c). The cumulative covariance (cid:80) β (cid:54) = α ˜ C αβ exceeds the vari-ance ˜ C αα in (d). The blue symbols represent experimentaldata and red lines represent numerical simulations based on t buff → end . channels are uniform with depth and small: for examplevar[ I ( z )] / (cid:104) I ( z ) (cid:105) < .
04 for all z .The experimentally observed fluctuations of individ-ual transmission eigenchannels are quantitatively repro-duced by the numerical simulations of t buff → end and t in Figs. 3(a,b,d-f). The excellent agreement betweenexperimental and numerical results confirms that eigen-channel fluctuations depend on their transmittance. Thehigher the transmittance, the lower the fluctuations.This means that high-transmission eigenchannels have arobust and consistent depth profile: irrespective of thedisorder configuration of a system.Finally, we investigate the cross-correlations betweendifferent transmission eigenchannels. For any given dis-order realization, eigenchannels are an orthogonal set offunctions at the front and back surfaces of the medium.While eigenchannels differ from realization-to-retaliation,their orthogonality implies that the differences in theirfield profiles must be correlated from realization-to-realization. This does not mean, however, that the inten-sity fluctuations of their profiles inside the sample shouldbe correlated. To study cross-correlations in the eigen-channels’ intensity-fluctuations across the sample, we in-troduce the covariance ˜ C αβ = (cid:104) δI α ( z ) δI β ( z ) (cid:105) z , where (cid:104) ... (cid:105) z describes both ensemble averaging and depth av-eraging. For α = β , ˜ C αα reduces to the variance ˜ C α which describes the eigenchannel fluctuations. Fig. 4(a-c) shows the experimental and numerical re-sults of ˜ C αβ for all α and β . The non-vanishing off-diagonal elements of ˜ C αβ ( α (cid:54) = β ) reveal coordinatedchanges in the eigenchannels’ depth-profiles. Betweendifferent pairings of eigenchannels, the correlations dif-fer. The larger the difference in the transmittances of apair, the weaker the correlation of their depth profile fluc-tuations. Furthermore, lower-transmission eigenchannelstend to correlate more with other low-transmission eigen-channels than higher-transmission eigenchannels do withother high-transmission eigenchannels. Quantitativelywe can describe the correlation of a single eigenchannelto all others by the cumulative covariance (cid:80) β (cid:54) = α ˜ C αβ . Asshown in Fig. 4(d), the cumulative covariance increaseswith α , indicating higher-transmission eigenchannels aremore independent from other eigenchannels than lower-transmission eigenchannels. Moreover, the cumulativecovariance exceeds the variance ˜ C αα = ˜ C α by a factor of2. Hence, the total cross-correlation for a single eigen-channel is stronger than its own fluctuation.In summary, we investigated realization-to-realizationfluctuations and correlations of the spatial structures oftransmission eigenchannels inside 2D diffusive systems.High-transmission eigenchannels exhibit low fluctuationsin their depth profiles (cross-section-averaged intensityvs. depth), when compared to low-transmission eigen-channels. While the intensity fluctuations are nearlyconstant -as a function of depth- for high-transmissioneigenchannels, they have a strong depth dependence forlow-transmission eigenchannels. Moreover, the changesin eigenchannel profiles from realization-to-realizationare correlated. The correlations are stronger betweeneigenchannels with similar transmittances. A lower-transmission eigenchannel tends to coordinate more withother eigenchannels than a higher-transmission eigen-channel. The overall correlation of any eigenchannelto all other eigenchannels exceeds its own fluctuations.Our findings of the second-order statistical propertiesof transmission eigenchannels are general and applica-ble to other types of waves such as microwaves, acousticwaves and matter waves. In practical applications, spa-tial wavefront shaping techniques enable the excitationof high-transmission eigenchannels which we now knowcan reliably be used to deliver energy deep inside a sin-gle diffusive sample: because of their consistency androbustness. ACKNOWLEDGMENTS
This work is supported partly by the Office of NavalResearch (ONR) under Grant No. N00014-20-1-2197,and by the National Science Foundation under GrantNos. DMR-1905465 and DMR-1905442. [1] A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink,Nat. Photon. , 283 (2012).[2] S. Rotter and S. Gigan, Rev. Mod. Phys. , 015005(2017).[3] C. W. Beenakker, Rev. Mod. Phys. , 731 (1997).[4] O. N. Dorokhov, Solid State Commun. , 381 (1984).[5] Y. Imry, Europhys. Lett. , 249 (1986).[6] P. Mello, P. Pereyra, and N. Kumar, Ann. Phys. ,290 (1988).[7] J. B. Pendry, A. MacKinnon, and A. B. Pretre, Phys-ica A , 400 (1990).[8] Y. V. Nazarov, Phys. Rev. Lett. , 2129 (1996).[9] I. M. Vellekoop and A. P. Mosk, Phys. Rev. Lett. ,120601 (2008).[10] M. Kim, Y. Choi, C. Yoon, W. Choi, J. Kim, Q.-H. Park,and W. Choi, Nat. Photon. , 581 (2012).[11] H. Yu, T. R. Hillman, W. Choi, J. O. Lee, M. S. Feld,R. R. Dasari, and Y. K. Park, Phys. Rev. Lett. ,153902 (2013).[12] M. Kim, W. Choi, C. Yoon, G. H. Kim, and W. Choi,Opt. Lett. , 2994 (2013).[13] S. M. Popoff, A. Goetschy, S. F. Liew, A. D. Stone, andH. Cao, Phys. Rev. Lett. , 133903 (2014).[14] J. Bosch, S. A. Goorden, and A. P. Mosk, Opt. Express , 26472 (2016).[15] C. W. Hsu, S. F. Liew, A. Goetschy, H. Cao, and A. D.Stone, Nat. Phys. , 497 (2017).[16] W. Choi, A. P. Mosk, Q. H. Park, and W. Choi,Phys. Rev. B , 134207 (2011).[17] B. G´erardin, J. Laurent, A. Derode, C. Prada, andA. Aubry, Phys. Rev. Lett. , 173901 (2014).[18] A. Pe˜na, A. Girschik, F. Libisch, S. Rotter, and A. Cha-banov, Nat. Commun. , 1 (2014).[19] M. Davy, Z. Shi, J. Park, C. Tian, and A. Z. Genack,Nat. Commun. , 6893 (2015).[20] R. Sarma, A. Yamilov, S. F. Liew, M. Guy, and H. Cao,Phys. Rev. B , 214206 (2015).[21] O. S. Ojambati, H. Yılmaz, A. Lagendijk, A. P. Mosk,and W. L. Vos, New J. Phys. , 043032 (2016).[22] R. Sarma, A. G. Yamilov, S. Petrenko, Y. Bromberg,and H. Cao, Phys. Rev. Lett. , 086803 (2016).[23] P. Hong, O. S. Ojambati, A. Lagendijk, A. P. Mosk, andW. L. Vos, Optica , 844 (2018).[24] H. Yılmaz, C. W. Hsu, A. Yamilov, and H. Cao,Nat. Photon. , 352 (2019).[25] S. F. Liew, S. M. Popoff, S. W. Sheehan, A. Goetschy,C. A. Schmuttenmaer, A. D. Stone, and H. Cao, ACSPhoton. , 449 (2016).[26] A. M. Paniagua-Diaz, A. Ghita, T. Vettenburg, N. Stone,and J. Bertolotti, Opt. Express , 33565 (2018).[27] M. Kim, W. Choi, Y. Choi, C. Yoon, and W. Choi, Opt.Express , 12648 (2015).[28] S. F. Liew, S. M. Popoff, A. P. Mosk, W. L. Vos, andH. Cao, Phys. Rev. B , 224202 (2014).[29] S. F. Liew and H. Cao, Opt. Express , 11043 (2015).[30] O. S. Ojambati, A. P. Mosk, I. M. Vellekoop, A. La-gendijk, and W. L. Vos, Opt. Express , 18525 (2016).[31] M. Koirala, R. Sarma, H. Cao, and A. Yamilov, Phys.Rev. B , 054209 (2017).[32] P. Fang, C. Tian, L. Zhao, Y. P. Bliokh, V. Freilikher,and F. Nori, Phys. Rev. B , 094202 (2019). [33] See Supplementary Materials at http://link.aps.org/ fordetailed descriptions of the sample design and fabrica-tion, the experiment, and numerical simulations, whichincludes Refs. [22, 34–39].[34] A. G. Yamilov, R. Sarma, B. Redding, B. Payne, H. Noh,and H. Cao, Phys. Rev. Lett. , 023904 (2014).[35] N. Bender, H. Ylmaz, Y. Bromberg, and H. Cao, APLPhoton. , 110806 (2019).[36] A. A. Lisyansky and D. Livdan, Phys. Rev. B , 14157(1993).[37] C. W. Groth, M. Wimmer, A. R. Akhmerov, andX. Waintal, New J. Phys. , 063065 (2014).[38] A. Yamilov, S. Petrenko, R. Sarma, and H. Cao, Phys.Rev. B , 100201(R) (2016).[39] R. Sarma, A. Yamilov, P. Neupane, B. Shapiro, andH. Cao, Phys. Rev. B , 014203 (2014). Supplementary:
I. SAMPLE DESIGN AND FABRICATION (a)
Waveguide Structure
Hole Size
500 nm500 nm
Photonic Crystal (c) (d)
Diffusive RegionBuffer
Taper
Buffer vs. Diffusive Region (b)
FIG. 5. Scanning electron microscope (SEM) images of a2D diffusive waveguide. In (a) we show a composite SEMimage which outlines the structures we etch into a silicon-on-insulator wafer when fabricating our structures. A SEMimage of the interface between the buffer and diffusive regionsis marked by the blue dashed line in (b). Close-up images ofthe photonic crystal sidewall and randomly-distributed holesare shown in (c) and (d).
Figure 5 shows a schematic of our two-dimensional(2D) disordered waveguide structures. The major com-ponents are the tapered waveguide, the buffer region, andthe diffusive region. The air holes (diameter = 100 nm),which induce light scattering in the buffer and diffusiveregions, are randomly distributed with a minimum (edge-to-edge) distance of 50 nm. The diffusive region has 5250holes, which results in an air filling fraction in the Si of5 . . µ m, length =15 mm). It then enters a tapered waveguide (taperingangle = 15 ◦ ). The tapered waveguide width decreasesgradually from 300 µ m to 15 µ m. The tapering results in waveguide mode coupling and conversion [22]. To avoidlight leakage, the tapered waveguide has photonic crystalsidewalls. II. OPTICAL SETUP
Beam
SplitterPolarizer MirrorLens ( 𝒇 ) MirrorLens ( 𝒇 ) x100NA 0.7 O b j . Obj.2
Beam
Splitter 𝝀/𝟐
PlateLens ( 𝒇 )IR CCD Sample Laser
SLM SlitShutter x100 NA . FIG. 6. A depiction of our experimental setup. Monochro-matic light from our laser is linearly polarized and split intotwo beams. One beam illuminates the phase modulating sur-face of a spatial light modulator (SLM), while the other isused as a reference beam. The SLM is used to control the in-put wavefront in our diffusive waveguide structures. A beamsplitter merges the light collected from the top of our samplewith the reference beam on an IR CCD. The focal length ofthe three lenses used in this setup are: f = 400 mm, f = 75mm, and f = 100 mm. Fig. 6 is a schematic of the experimental setup.Continuous-wave (CW) output from a tunable laser(Keysight 81960A) -operating around 1554 nm- is lin-early polarized and split into two beams. One beam il-luminates the phase modulating surface of a phase-onlySLM (Hamamatsu LCoS X10468), while the other is usedas a reference beam. A one-dimensional (1D) phase-modulation pattern is displayed on the SLM, consistingof 128 macropixels. Each macropixel consists of 4 × f = 400 mm and f = 75 mm, we image thefield on the SLM plane onto the back focal plane of along-working-distance objective ( Obj. 1 ) (Mitutoyo MPlan APO NIR HR100 × , Numerical Aperture = 0.7). Toprevent the unmodulated light from entering the objec-tive lens, we display a binary diffraction grating withineach macropixel to shift the modulated light away fromthe unmodulated light in the focal plane of the f lens.Using a slit in the same focal plane, we block everythingexcept the phase-modulated light in the first diffractionorder. Before the the f lens, we insert a half-wave ( λ/ Obj. 1 and illuminated with the Fourier transform of the phase-modulation pattern displayed on the SLM. From the topof the wafer, a second long-working-distance objective(
Obj. 2 ) (Mitutoyo M Plan APO NIR HR100 × ) col-lects light scattered out of plane of the waveguide. Weuse a third lens with a focal length of f = 100 mm to-gether with Obj. 2 to magnify the sample image by ×
50. With a second beam splitter, we combine the lightcollected from the sample and the reference beam. Theirinterference patterns are recorded with an IR CCD cam-era (Allied Vision Goldeye G-032 Cool).
III. INTERFEROMETRIC MEASUREMENT
With the interferometric setup described in the lastsection, we can measure the field distribution of lightscattered out-of-plane from within the diffusive waveg-uide: for any phase-modulation pattern displayed on theSLM. To do this, we first measure the 2D intensity dis-tribution of the scattered light by blocking the referencebeam with a shutter. Then using the reference beam inour setup, we retrieve the phase profile of the scatteredlight with a four-phase measurement (wherein the globalphase of the pattern on the SLM is modulated four timesin increments of π/ t slm → buff and t slm → end , which map the field from the SLM surface tothe buffer and to the far end of the disordered waveguide,respectively. To construct the matrix relating the field inthe buffer region to the field near the end of the diffusivewaveguide, t buff → end , we define the field-mapping matrixbetween the two regions t buff → end ≡ t slm → end t − → buff .To calculate the inverse of t slm → buff , we use Moore-Penrose matrix inversion. In this operation we only takethe inverse of the 55 highest singular values of t slm → buff ,and set the inverse of the remaining singular values tozero. This restriction is imposed because our diffusivewaveguide only has 55 transmission eigenchannels. IV. DETERMINATION OF TRANSPORTPARAMETERS
Diffusive wave propagation in a scattering mediumwith loss is determined by two parameters: the transportmean free path (cid:96) t and the diffusive dissipation length ξ a . In a 2D system, the latter can be expressed as ξ a = (cid:112) (cid:96) t (cid:96) a /
2, where (cid:96) a is the ballistic dissipation length.To determine ξ a and (cid:96) t in the diffusive region of the 2Dwaveguide, we first measure the cross-section-averagedintensity I ( z ) as a function of depth z for multiple ran- Depth z ( m m) I n t en s i t y I ( z ) TheoryExperiment l = 3.2 m m t ξ = 28 m m a FIG. 7. Determining the transport mean free path (cid:96) t andthe diffusive dissipation length ξ a of the diffusive waveguidesby fitting the experimentally-measured average depth profilefor random incident wavefronts, (cid:104) I ( z ) (cid:105) , (blue) to theoreticalpredictions from the diffusion equation (red). dom input wavefronts. We then ensemble average thedata, (cid:104) I ( z ) (cid:105) , and fit the theoretically-predicted depthprofile -based on the diffusive equation- to it. The the-oretical (cid:104) I ( z ) (cid:105) is found by convolving the incident bal-listic intensity I exp[ − z/(cid:96) s ] ( (cid:96) s represents the scatteringmean free path, and (cid:96) s ≈ (cid:96) t in our case), which actsas the source, and the Green’s function of the diffusionequation [36]: G ( z, z (cid:48) ) = (cid:26) P ( z ) P ( L − z (cid:48) ) , z < z (cid:48) P ( z (cid:48) ) P ( L − z ) , z > z (cid:48) (1)where P ( z ) = sinh( z/ξ a ) + z /ξ a cosh( z/ξ a ), and z =( π/ × (cid:96) t is the so-called extrapolation length.We compute the difference between the experimentaland theoretical (cid:104) I ( z ) (cid:105) for different values of ξ a and (cid:96) t ,and identify the minimum difference at ξ a = 28 µ m and (cid:96) t = 3 . µ m. Fig. 7 shows an excellent agreement betweenthe measured (cid:104) I ( z ) (cid:105) and the theoretical prediction.In the buffer region, the air hole density is 10 timeslower than in the diffusive region. Thus, the transportmean free path is 10 times longer, (cid:96) buff t = 32 µ m. Theloss, caused by out-of-plane scattering from the air holes,is also 10 times weaker, thus the ballistic dissipationlength (cid:96) a is 10 times longer. This leads to a tenfold in-crease in the diffusive dissipation length: ξ buff a = 280 µ m. V. EIGENCHANNEL PROFILEMEASUREMENT
In total, we measure the transmission eigenchannel in-tensity profiles of 13 independent realizations. We obtainthese measurements from two samples with different ran-dom hole configurations. To generate independent sys-tem realizations from the same random hole configura-tion, we vary the wavelength of the input light beyond thespectral correlation width of the diffusive region: 0.4 nm.Over a wavelength span of 3 nm, we vary the input wave-length of our laser in increments of 0.5 nm. We choose thespecific wavelength range of the measurement -for eachrandom hole configuration- such that the effective dissi-pation in the diffusive region is minimal over the wave-length range and homogeneous. While our waveguidestructure has a width of 15 µ m, we only use the central10 µ m region of the waveguides out-of-plane-scatteredlight when performing our measurements to avoid arti-facts from light scattered out-of-plane from the photoniccrystal boundaries. VI. NUMERICAL SIMULATIONS
In our numerical simulations, we use the recursiveGreen’s function method in the Kwant simulation pack-age [37]. We simulate a two-dimensional (2D) rectangu-lar waveguide geometry, which is defined using a tight-binding model for scalar waves on a square grid. At thewaveguide boundaries, which are reflective, the grid isterminated. The leads are attached to the open endsof the waveguide, allowing for computation of the com-plete scattering S matrix of the system and the wave fieldthroughout the bulk of the system: under an excitationby an arbitrary combination of field amplitudes for thepropagating modes. The width W of the simulated sys-tem is selected so that the number of waveguide modes N matches the number found in the experiment. Once W is chosen, the length of the disordered waveguide isdetermined by the ratio L/W of the waveguides used inthe experiment. Due to the low filling fraction of the airholes in the experimental waveguides, both in the bufferregion and in the main disordered region, we assume thatthe number of propagating modes is equal to N .Scattering is introduced by a randomly (box distribu-tion) fluctuating real-valued on-site ‘energy’ in the tight-binding model, see Refs. [22, 38]. The addition of a pos-itive imaginary constant to the same term simulates theeffect of absorption. In our previous works, we confirmedthat the process of vertical leakage due to the holes in ourdisordered waveguides can be modeled via absorption ina 2D system [22, 34, 39]. The actual material absorptionin our experimental system is negligible. By a properchoice of these parameters, we can match the experimen-tal values for the transport mean free path (cid:96) t and thediffusive dissipation length ξ a .To model the weakly scattering ‘buffer’ region, we re-duce the scattering (the amplitude of the on-site fluctu-ation) so that transport mean free path is reduced by afactor of 10. The latter corresponds to a 10 times re-duction in the areal density of the air holes in the bufferregion. Furthermore, because the out-of-plane scatteringloss is reduced 10 times, the diffusive dissipation lengthis also reduced by the same factor.The buffer region is incorporated into the experimentalwaveguides to measure t buff → end of the diffusive waveg-uide, which is not a direct measurement of the field trans- mission matrix t . We numerically simulate the eigen-chanels of both matrices to confirm their depth profilesare equivalent. The matrix t is obtained from the in-cident and transmitted fields in the left and right leadswithout the buffer. To compute t buff → end , we computethe auxiliary matrices t in → buff and t in → end . The formermatrix relates the incident fields in the left lead to thefields at 2 × N randomly selected points within 10 µ m × µ m region centered in the buffer region (of an area15 µ m × µ m). The chosen points are at least 2.5 µ m separate from each other or any boundary/interface.The second auxiliary matrix t in → end relates the imping-ing fields in the left lead to the fields at 2 × N randomlyselected points within 10 µ m × µ m region at the end ofthe diffusive waveguide. Again all points are at least 2.5 µ m (which is on the order of (cid:96) t ) spaced. In the last step,we compute t buff → end = t in → end t − → buff , where t − → buff is calculated with the Moore-Penrose pseudo-inverse.To calculate the spatial structure of the transmissioneigenchannels, we perform a singular value decomposi-tion on the t matrix, and use the right singular vectorsas input fields in the left lead to excite individual eigen-channels. For the matrix t buff → end , its right singular vec-tors are transformed to the incident fields in the left leadby multiplying t − → buff . To further mimic the phase-onlymodulation of the SLM in the experiment, we only keepthe phases of the incident fields, and set the field mag-nitudes equal. We calculate all eigenchannels for t and t buff → end for an ensemble of 1000 disorder configurationsof the waveguides. The numerical results are presentedin Figs. 2-4 in the main text.To compare the variance ˜ C α and covariance˜ C αβ numerically calculated from t buff → end to theexperimentally-measured ones, we need to account forsome experimental limitations and imperfections. Onone hand, the finite spatial resolution of our detectionoptics effectively enlarges the speckle grain size of thefield measured inside the diffusive waveguide. Thisreduction in the number of speckle grains increases thefluctuations of the cross-section-averaged intensity. Onthe other hand, the combined effects of sample driftduring measurements and the presence of two linearpolarizations in the light scattered out-of-plane from oursample; decrease the fluctuations of the cross-section-averaged intensity. For random incident wavefronts,the spatially-averaged intensity variance of our exper-imental measurements is var[ I ( z )] = 0 .
59, comparedto var[ I ( z )] = 0 .
64 from the numerical simulations of t buff → end . For all eigenchannels, we re-scale the nu-merical var[ I α ( z )] and ˜ C αβ by the multiplicative factor0 . / .
64, in order to compare them to the experimentalvalues. While we applied the re-scaling factor to thefluctuations and correlations calculated from numericalsimulations of t buff → end , we did not apply it to theresults from simulations of t . VII. IDENTIFICATION OF EXPERIMENTALEIGENCHANNELS
10 20 30 40 50
Experimental eigenchannel index E N u m e r i c a l e i gen c hanne l i nde x (b) T r an s m i ss i on e i gen v a l ue -5 -10 -15 Numerical eigenchannel index (a)
FIG. 8. Calculated transmission eigenvalues, as a function ofeigenchannel index α , are shown in (a). In (b), we show themapping between the experimentally-measured eigenchannelprofiles with index α E and the first 22 (in the order of de-creasing transmittance) eigenchannels with index α found inthe numerical simulations based on the t buff → end matrix. In this section, we analyze the normalized eigenchan-nel profiles measured in the experiment (cid:104) I α E ( z ) (cid:105) andthe numerical simulations (cid:104) I α ( z ) (cid:105) . For each experimen-tal channel with an index of α E ∈ [1 ... (cid:82) L ( (cid:104) I α E ( z ) (cid:105) − (cid:104) I α ( z ) (cid:105) ) dz . Wedo not use any channel-specific adjustments/fits in thisidentification. This process gives the mapping of α E to α , shown in Fig. 8. A few experimental eigenchannelsare redundant, particularly in the range α ∈ [6 ... α >
22 are observed experimen-tally. We attribute this to the finite signal-to-noise ratioin the experimental data. The eigenchannels with α >
22 have a transmittance less than ∼ . (cid:104) ... (cid:105) for the α -th channel thatcorresponds to multiple α E ’s include both disorder con-figuration and { α E } averaging. VIII. EIGENCHANNEL VARIANCE
10 20 30 40 500.050.10.150.2 ~ ~ F l u c t ua t i on s Effect of normalization on fluctuations
Eigenchannel index
FIG. 9. We compare the variance of eigenchannel pro-files calculated with the normalized intensity (purple solidline) to the intensity variance normalized by the mean inten-sity (cid:104) var[ ˜ I α ( z )] / (cid:104) ˜ I α ( z ) (cid:105) (cid:105) z squared (green dashed line). Theyshow similar growth with the eigenchannel index α . In the main text, we present the realization-to-realization fluctuations of the egenchannels normalizedintensity profiles. For each eigenchannel, the mea-sured intensity profile ˜ I ( z ) is normalized to I ( z ) =˜ I ( z ) / [(1 /L ) (cid:82) L ˜ I ( z (cid:48) ) dz (cid:48) ]. Using a different normaliza-tion procedure, we check the effect of our normaliza-tion on the eigenchannel fluctuations using numericalsimulations of t . For an eigenchannel α , the variancevar[ ˜ I α ( z )] = (cid:104) δ ˜ I α ( z ) (cid:105) of the unnormalized intensity fluc-tuation δ ˜ I α ( z ) = ˜ I α ( z ) − (cid:104) ˜ I α ( z ) (cid:105) can be normalized bydividing the square of the mean intensity (cid:104) ˜ I α ( z ) (cid:105) atthe same depth z . Then this ratio var[ ˜ I α ( z )] / (cid:104) ˜ I α ( z ) (cid:105) can be averaged over all z . In Fig. 9, we compare thisquantity to the variance of the normalized intensity pro-file, ˜ C α , calculated in the main text. Both exhibit anincrease with the eigenchannel index αα