Coupling of energy into the fundamental diffusion mode of a complex nanophotonic medium

We demonstrate experimentally that optical wavefront shaping increases light coupling into the fundamental diffusion mode of a scattering medium. The total energy density inside a scattering medium of zinc oxide nanoparticles was probed by exciting fluorescent spheres that were randomly positioned in the medium and collecting the fluorescent power. We optimized the incident wavefront to obtain a bright focus at the back surface of the sample and found that the concomitant fluorescent power is enhanced compared to a non-optimized incident wavefront. The observed enhancement increases with sample thickness. Based on diffusion theory, we derive a model wherein the distribution of the energy density of wavefront-shaped light is dominated by the fundamental diffusion mode. Our model agrees remarkably well with our experiments, notably since the model has no freely adjustable parameters.


Introduction
Diffusion is a process that leads to a uniform spreading of matter or energy as a result of randomness. Numerous physical phenomena are described by diffusion, ranging from colloidal particles to currencies to animal motion [1][2][3][4][5][6][7][8]. Many wave phenomena are also well described by diffusion: upon averaging over the disorder, interference is washed out [5,9,10]. Hence, waves become diffuse after a distance of the order of one transport mean free path ℓ, i.e., the distance over which the memory of the incident light direction is lost [5]. Of particular relevance to diffusing waves is the energy density that is defined as the amount of energy in a wave field that is stored in a given volume [3]. Due to diffusion, the ensemble-averaged energy density of the wave W d has a characteristic shape shown in figure 1(a) [3,7,[11][12][13][14]. The gradient of the energy density at the exit surface of the scattering medium is related to the transport of energy, which yields the wave-equivalent of the well-known Ohm'(s) law » ℓ T L, with L the thickness of the medium. A fundamental question addressed here is whether it is feasible to couple energy into a diffusion eigensolution, which has a greater energy density than that of the usual non-optimized energy density W d , as shown in figure 1(a). In case of light, an enhanced energy density inside the scattering medium is important for applications such in white LEDs [15][16][17][18], solar cells [19][20][21], random lasers [22][23][24], and biomedical optics [25].
The ensemble-averaged energy density W is described by the diffusion equation as , where D is the diffusion constant. For a real and finite medium, such as the widely studied slab geometry, there is translational invariance in the r = ( ) x y , plane [12]. The steady-state solution ( ) W z d along the direction of light propagation z is given by [26], where m is the eigensolution index, C m is the relevant coefficient (see (B.9) in appendix B.2), and = + + L L z z ex e1 e2 is the effective thickness of the sample, where z e1 and z e2 are the extrapolation lengths at the front and back surfaces of the sample, respectively. Figure 1(b) shows the first three eigensolutions. The fundamental eigensolution m=1 is the only physical solution with a positive energy density along the sample depth z, as shown in figure 1(b). When a plane wave is incident on a scattering medium, energy is coupled into all eigensolutions that add up to the usual energy density W d shown in figure 1(a) [27]. Recently, it has been demonstrated that waves in complex media can be exquisitely controlled by shaping the spatial phase of the incident waves, notably by feedback-based wavefront shaping [28][29][30][31][32][33][34], time reversal [35,36], phase conjugation [37,38], and transmission matrix-based control [39][40][41]. In wavefront shaping, the spatial phase of the field incident on the scattering medium is controlled to enhance the intensity several hundred times in a single speckle spot (area l p A 2 2 ) at the back surface of the sample [42]. With this method, the total transmission T through a scattering medium can be made to differ from Ohm'(s) law [30,33]. Therefore, it is relevant to investigate whether the energy density inside a scattering medium can be controlled. In particular, we wonder if wavefront shaping couples energy into the diffusion fundamental mode, which has a higher energy density. To date, only numerical calculations [43][44][45] and a single-realization elastic wave experiment [46] have addressed the distribution of the energy density inside two-dimensional (2D) scattering media. Control of the optical energy density inside a three-dimensional (3D) scattering medium has so far not been explored.
In this paper, we experimentally demonstrate the coupling of light predominately into the fundamental diffusion mode by wavefront shaping. Since it is difficult to probe waves inside a 3D material, we use embedded fluorescent nanospheres as reporters. The incident shaped waves excite the fluorescent spheres, whose fluorescent intensity is proportional to the local energy density. By averaging over many fluorescent spheres distributed throughout the samples, we probe the total internal energy of the 3D scattering medium. We observe that the total optical energy is enhanced when the wavefront of the incident light is optimized. The enhancement increases with sample thickness. To interpret our results, we propose a parameter-free model wherein the energy density of wavefront-shaped light is described by the fundamental eigensolution of the diffusion equation, which agrees well with the experimental observations.

Samples and preparation
In our experiments, we studied scattering media that are slabs of zinc oxide (ZnO) nanoparticles (Sigma-Aldrich Zinc Oxide 205532, average grain size of 200 nm). Using a professional airbrush (Harder & Steenbeck Airbrush Evolution Silverline, 123003), we sprayed an aqueous suspension of 2.5 g of the ZnO nanoparticles and a suspension of the fluorescent particles, which are dye-doped polystyrene spheres (Thermo Scientific R50 Fluoro-Max Dyed Red, diameter = d 50 nm). The concentration of the fluorescent spheres is m -1.5 m 3 in 7.5 ml of deionized water. The two suspensions were thoroughly mixed and then spray painted on a cleaned microscope glass slide of thickness 0.17 mm with lateral dimensions of 25 mm×50 mm. In order to obtain a The parameters used in these plots are: well-defined thickness of our samples, we automated the spray-painting process by using a LEGO mindstorm kit. The samples have thicknesses L ranging from 2 to 22 mm, as a result of varying the spraying time. Inside the samples, the fluorescent spheres are randomly dispersed (see figure 2) to probe the total energy density inside the ZnO scattering medium. The fluorescent spheres are excited by incident light with wavelength l = 561 nm i and emit fluorescent light at a different wavelength l = 612 nm f . In order to ensure that the distribution of the energy density at l i inside the scattering medium is not perturbed by the absorption from the probing fluorescent spheres [45,47], we use fluorescent spheres with low absorption probability of 2´-10 % 4 [48]. The absorption probability is the ratio of the absorption cross section to the physical cross section of the fluorescent spheres. Moreover, the concentration of the fluorescent spheres in the samples is as low as m -0.35 m 3 , which results into an albedo » -a 1 10 6 averaged over all the fluorescent spheres [49]. The total fluorescent power emitted from the sample fluctuates by 10% from different positions on the sample, which shows that the fluorescent particles are fairly homogenously distributed.

Experimental set-up
The experimental set-up is shown in figure 3. The coherent light source is a diode-pumped solid state laser (Cobolt Jive, 100 mW continuous wave output at 561 nm). We control the laser power using a combination of a half-wave plate (l 2) and a polarizing beam splitter (PBS). A beam expander with magnification of 3 expands the beam and lens L1 ( f=11 mm) focuses the beam onto the front surface of a single mode fiber (SMF) (Thorlabs P3-460B-FC-1). The SMF spatially filters the laser beam in order to prevent pointing instability of the beam. A fiber collimator C (Kirchoff 60FC-L-0-M75-26) collimates the spatially cleaned beam and it passes through an aperture A1 to illuminate a phase-only liquid-crystal spatial light modulator (SLM) (Holoeye PLUTO-VIS-014-C). The diameter of the beam on the SLM is 8 mm. We used the piece-wise sequential algorithm described in [29] to find an optimized incident wavefront. The image of the SLM is demagnified and imaged to the pupil of a microscope objective MO1 (Nikon: infinity corrected, 100×, NA=0.9) through a telescope consisting of lenses L2 ( f=200 mm) and L3 ( f=100 mm). An aperture A2 filters out the higher diffraction orders of the SLM'(s) pixels and transmits only the 0th order. The demagnified image of the SLM is focused onto the surface of the sample by the microscope objective MO1 to a diffraction-limited spot of 620 nm, when wavefront of light is flat. The back surface of the sample is imaged onto a charge-coupled device (CCD) camera (Allied Vision Technologies Stingray F-145) and an electron multiplying CCD (EMCCD) camera (Andor iXon Ultra 897) by an oil-immersion microscope objective MO2 (Nikon: Infinity corrected, 60×, NA=1.49) and an achromatic lens L4 ( f=500 mm). A dichroic beam splitter DB (Semrock Di02-R561-25×36) reflects the excitation green beam and transmits the emitted fluorescent intensity from the probes. The dichroic beam splitter has a transmission of 93% from 578.4 to 1200 nm. A combination of a notch filter (Semrock NF03-561E-25) and a single-band bandpass filter F2 (Semrock FF01-620/52-25) blocks the remnant excitation green light. The green light is detected by the CCD camera. A polarizer P blocks the orthogonal polarization that is not controlled by the SLM. The intensity of the laser beam is reduced by a neutral density filter F1 with an optical density of 2 to prevent saturating the CCD camera. Setup to probe the total energy density inside a scattering medium. The scattering medium is a slab of disordered ZnO particles in air. The medium is illuminated with a shaped incident wavefront such that incident light at l i (green intensity) is optimized in a speckle spot at the back surface of the sample (see inset). The scattering medium is lightly doped with randomly positioned fluorescent spheres to probe the energy density inside the sample. The total fluorescent power emitted by the fluorescent spheres at l f (red intensity) is detected by the EMCCD (see inset).

Measurement procedure
To obtain ensemble-averaged data that can be interpreted with theory, we performed wavefront shaping experiments at 100 to 130 different positions on each sample. At each position, we determine the fluorescent where P f o is the total fluorescent power with an optimized phase pattern on the SLM, and P f n is the total fluorescent power for a reference incident wavefront. As a reference we use a phase pattern stored from an optimization on a very different part of the sample. This reference field experiences the same diffraction efficiency and surface reflection as our optimized field, in contrast to a plane wave.

Fidelity of optimization
The incident optimized field in any wavefront shaping experiment is a superposition of the ideally optimized incident field, and a noise field that is by definition orthogonal to the ideal field [30]. The fidelity g | | 2 quantifies the overlap between the experimentally obtained incident field and the ideal incident field that optimally couples light to the target spot [30]. A perfectly optimized field has g = | | 1 2 , while all imperfections such as limited spatial control, temporal decoherence or modulation noise, inevitably reduce g | | 2 . Here we neglect a correction of the order of N 1 , where N is the number of transmitting channels » N 10 5 , see [50]. The measured fluorescent power also consists of a linear combination of the ideal incident field term and a noise term. Since the noise field is unrelated to the sample properties, its response is on average equal to that of an unoptimized incident field with the same power. Thus the fluorescent power enhancement h f necessarily depends on fidelity g | | 2 as Here h f p is the fluorescent power enhancement at perfect fidelity g = (| | ) 1 2 due to a perfectly shaped incident field. The second term on the right-hand side of (3) is the contribution from the noise fields. To measure the fidelity of each generated field, we exploit the fact that g | | 2 is well approximated by the ratio of the power in the target spot (figure 2, inset) to the average total transmitted power for an unoptimized incident field [30]. This approximation holds up to the level of fluctuation of the total transmittance, which is about 2% [51]. In practice, determination of fidelity is limited by experimental noise estimated to be about 7%. Factors such as inhomogeneous sample thickness, experimental noise, and instability in environmental conditions limit g | | 2 , and result in variations thereof [52,53]. Although variations are not desired, they have the advantage of yielding a range of g | | 2 to investigate. In each experiment at a different position on the sample, we determined both the fidelity g | | 2 and the fluorescent power enhancement h f . The back surface of the sample is imaged onto the EMCCD camera, which detects the fluorescent intensity and a CCD camera, which detects the excitation green laser light. The dichroic beam splitter separates the fluorescent intensity from the excitation beam. l 2 is a half-wave plate, PBS is a polarizing beam splitter, BD is a beam dump, BE is a beam expander, L1, L2, L3, and L4 are lenses with a focal distances of 11 mm, 200 mm, 100 mm, and 500 mm respectively, SMF is a single mode fiber, C is a collimator, A1 and A2 are apertures, M1 and M2 are mirrors, MO1 and MO2 are microscope objectives, DM is a dichroic beam splitter, F1 and F2 are filters, and P is a polarizer. | | 2 has increased to about 0.035. (For the 23 μm thick sample, the fidelity is limited due to the large number of transmitting channels, of the order of 10 5 , that is much greater than the number of degrees of control on the SLM, of the order of 10 3 .) This fluorescent enhancement implies that the total energy density for optimized incident wavefronts is higher than for unoptimized incident wavefronts. According to (3), h f is directly proportional to g | | 2 with a slope a h = -1 f p and an intercept of 1. Physically, the intercept means that for an

Experimental result
). Using linear least-squares, we therefore make a linear regression of the data with only the slope as adjustable parameter. For this sample, we obtain a slope a =  3.6 0.48, where the error margin gives the 95% confidence level. The residuals of the data are randomly distributed about the model, see figure 4(b), with more than 60% of the data within one standard deviation as expected for normally distributed errors [54].    (4). The uncertainty in the fluorescent power enhancement dh f p is plotted as the two green dashed curves, which were calculated using equation (5). The black dashed-dotted curve is for an invariant distribution of energy density along the sample depth. Parameters to obtain the theoretical curves are given in table 1.
In figure 5, we plot results as a function of sample thickness. The fluorescent power enhancement at perfect fidelity h f p increases with sample thickness L, which means that wavefront shaping causes more energy to be stored in the bulk of the medium. The horizontal error bars denote the standard deviation of the measurement of the sample thickness at different positions on the sample, while the vertical error bars denote uncertainty in the determination of h f p . For perfect fidelity, the total fluorescent power inside a 22.8 μm thick sample is 4.6 times greater than for non-optimized light.

Theoretical model
To interpret our experimental results, we employ diffusion theory (see appendix B for a full description of the model). For light with an optimized incident wavefront, the distribution of the energy density W o inside the medium is a priori unknown. The energy density within the boundary domain of the sample must be a linear combination of the diffusion eigensolutions, since the diffusion eigensolutions form a complete set. In the case of optimized light, the energy density W o must also be a linear combination of diffusion eigensolutions. The diffusion eigensolutions that describe optimized light are expected to have the highest positive contribution to the total transmission, since wavefront shaping is known to increase the total transmission [30,33]. Here, we model the energy density of optimized light W o as the fundamental eigensolution (m=1) of the diffusion equation for three reasons. Firstly, the contribution to the total transmission 1 of the first six eigensolutions in figure 6 shows that the fundamental eigensolution (m=1) contributes most to the total transmission. The m=1 contribution is even greater than the total transmission itself, which is a summation of all eigensolutions. Secondly, the fundamental eigensolution is the only physical solution with a positive energy density along the sample depth z, as shown in figure 1(b). Thirdly, a symmetric function peaked near the sample center, similar to the fundamental eigensolution, has been obtained from numerical calculation of wave-front shaped light [43], theoretical and numerical calculations [44,45], as well as a single-realization experiment [46] of high transmitting channels. Hence we neglect small corrections due to coupling to other modes than the fundamental eigensolution.
The prefactor of the energy density of wavefront-shaped light in our model is determined by the known total transmission of 2/3 for perfectly wavefront-shaped light [30,55]. The energy density of optimized light W o thus obtained is much larger than the diffusive energy density W d , as shown in figure 1(a).  . 1 The total transmission is the ratio of the total flux at the back surface of the sample (z=L) to the incident flux and is equal to i , z inj is the injection depth at which the incident light becomes diffuse, which accounts for the angular distribution of the incident shaped wavefront [56]. With an incident numerical aperture of 0.9, we find l = ℓ ( ) z 0.86 inj i . The fluorescent power enhancement h f p is determined by factors such as transport mean free path ℓ and the extrapolation lengths z e1 and z e2 , which are determined experimentally (see table 1). It therefore follows that the evaluation of h f p also has an uncertainty dh f p . We calculate dh f p using To compare our model to our experimental results, we plot in figure 5 the analytic model for h f p versus sample thickness L. Within the uncertainty of the samples' parameters, our model agrees remarkably well with our experimental results, especially in view of the fact that the model has no adjustable parameters. The excellent agreement between our model and our experimental results confirms that energy has been mostly coupled into the fundamental (m=1) diffusion mode.
In our experiments we obtain the fluorescent power enhancement h f p rather than the energy density enhancement h ed that is defined as , 6 ed f where ( ) O l L represents higher orders, see appendix B. The relation in (6) shows that the total fluorescent power depends on the total energy density inside the medium. Hence, the observed increase of the fluorescence is indeed a probe of the increased energy density.

Summary
We have measured the increase of the total energy density inside a scattering medium that is the result of a shaped incident wavefront. The results agree quantitatively with a model that considers energy to be dominantly coupled to the fundamental mode of the diffusion equation. Our results are of broad relevance since they apply to other wave control methods in scattering media, such as time reversal, phase conjugation, and transmission matrix-based control, as well as to other types of waves such as microwaves, acoustic waves, elastic waves, surface waves, and electron waves. We expect our results to be relevant for applications that require enhanced total optical energy density, such as efficient light harvesting in solar cells especially in near infrared where silicon has low absorption; for enhanced energy conversion in white LEDs, to reduce the quantity of expensive phosphor; for low threshold and higher output of random lasers; and in homogeneous excitation of probes in biological tissues. Last but not least, it will be fruitful to investigate possible relationships between the fundamental diffusion eigensolution and the universal diffusion time as obtained in [7, 59, 60]. Table 1. Values of the transport mean free path ℓ, extrapolation lengths z e1 and z e2 for wavelengths of the incident and fluorescent intensities. The transport mean free path was determined from total transmission measurement of similar samples (see [57]). Following [58], we derived the extrapolation lengths from the average reflectivities at the interfaces for a scattering medium with an effective refractive index n eff . The effective refractive index n eff was determined from measurement of the angle-resolved escape function on similar samples [57].

Incident intensity
Fluorescent intensity @ l i =561 nm @ l f =612 nm In order to determine s lc , we measured the total fluorescent power emanating from the sample for 12 h. We determine the laser and camera noise s lc from the standard deviation of the fluorescent power with respect to the mean power and we obtain s = 0.03 lc . The mesoscopic intensity fluctuation s m is given by the fluctuation of the overlap integral between the illuminating speckle intensity and the spatial distribution of the fluorescent spheres in the scattering medium. We estimated s m using a simple model. For a given number of spheres N sphere we generated a random array with N sphere non-zero values which represent our spheres embedded inside the scattering medium. To simulate the speckle intensity inside the scattering medium, we generated a second random array consisting of numbers with a Rayleigh distribution. We multiply the two random arrays and take the summation of the resulting array, which gives the simulated total fluorescence power P f s . We repeated the procedure to simulate 100 realizations and calculated the standard deviation of the normalized total fluorescence power á ñ P P , as mentioned in the main text.
A.1.1. Error in the estimation of the fidelity g | | 2 . In the main manuscript, we determine the fidelity g = á ñ | | where P o is the power in the optimized focus and á ñ P n is the ensemble-averaged power transmitted before optimization. Therefore, the error in estimating the fidelity s g where s C2 is due to expected fluctuation in the total transmitted intensity, i.e. the C 2 intensity correlation [12], s s is due to the roughness of the sample, and s c is the fraction of intensity cropped out by the optics. We estimated s C2 as where N is the number of transmitting channels and T tot is the total transmission [56].  figure 3(b), about 60% of the data points are within one standard deviation from the model, which means is expected since the data points are Gaussian distributed.

Appendix B. Theory
Here, we present the theory to model the spatial distribution of wavefront-shaped light inside a scattering medium. We model both the fluorescence enhancement and the energy density enhancement.

B.1. Assumptions
We first state the three (3) assumptions that we made. Firstly, we assume that our scattering medium can be described by a slab geometry with infinite lateral dimensions in the x and y directions and one boundary at z=0 and z=L, since the beam propagates through a very small area on the sample compared to the whole extent of the sample. For a 10 μm thick sample, the diameter of the diffuse spot at the back surface of the sample is about 20 μm, much smaller than the 40 mm extent of the whole sample. Therefore, we can safely assume that the boundaries in both x and y directions have negligible effect on light diffusion inside the medium. The boundaries along z, however, have an effect on light propagation. We therefore solve the one-dimensional diffusion equation for only the z direction [61].
Secondly, we assume that absorption is negligible in the samples. This assumption is valid since the albedo » a 1. However, weak absorption is present in our samples due to the presence of the low density of probe spheres. Therefore, the assumption that there is no absorption in the medium is an extremely mild limitation to our model. If necessary, the diffusion equation can be extended to account for absorption, see [62].
Thirdly, we assume that multiple scattering events occur inside the medium before light exits the sample, such that our samples can be modelled using diffusion theory. Therefore where ℓ is the transport mean free path. Also, we assume that where z z , e1 e2 are the extrapolation lengths on the front and back side of the sample respectively. The extrapolation lengths account for reflections on each side of the sample [63]. The extrapolation lengths are obtained from the effective refractive index of the scattering medium and the refractive index of the surrounding medium [63]. The inequality in (B.2) implies that the extrapolation lengths are small compared to the sample thickness, which is the case for our samples with thickness L ranging from ℓ 3 to ℓ 37 , where m = ℓ 0.6 m. We note that in [64], it has been shown that diffusion theory is robust to describe samples even only ℓ 2 thick to an accuracy of 2%.

B.2. Spatial distribution of the energy density of wavefront-shaped light
The spatial distribution of the energy density of wavefront-shaped light inside the medium is a-priori unknown. As stated in the main manuscript, we postulate that wavefront shaping predominately couples energy into the fundamental eigensolution of the diffusion equation. The one-dimensional diffusion equation is given by where ( ) W z t , is the ensemble-averaged energy density and D is the diffusion constant that is equal to with v E the energy velocity [65]. We use the boundary condition [12] and the source is taken to be the result of the exponential decay of the incident coherent beam where I 0 is the incident intensity and z inj is the injection depth at which the incident light becomes diffuse that accounts for the angular distribution of the incident shaped wavefront [56,66]. We solve (B.3) using the method of separation of variables [26], to obtain the energy density We define the diffusion coefficient We plot the first 3 eigensolutions in figure 1(b) in the main manuscript. The series in (B.8) can be analytically evaluated [67] to obtain In figure 1(a) in the main manuscript, we plot the unoptimized energy density W d using the same parameters as in figure 1(b).

B.3. Contribution of eigensolutions to the total transmission
We have derived the spatial distribution of energy densities W(z) of the diffusion eigensolutions. From the distribution, we want to obtain the contribution of each eigensolution to the total transmission and then find the eigensolution with the greatest contribution to the total transmission. The total transmission is the ratio of the total flux at the back side of the sample (z=L) to the incident flux [3] and it can be expressed as Using similar parameters used to plot figure 1 except for L=5 mm, we obtain the contributions to the total transmission as presented in figure 5 in the main manuscript. The fundamental eigensolution m=1 has the greatest contribution to the total transmission, even more than the total transmission. We therefore hypothesize that the energy density distribution of optimized light is identical to the fundamental eigensolution of the diffusion equation, which is given by

B.4. Enhancement of fluorescent power
We have derived the spatial distribution of the energy density of optimized W o and unoptimized light W d . In our experiments, we probe an integration over these distributions by using fluorescent spheres positioned inside the sample. In this section, we model the enhancement of the fluorescent power measured at the back surface of the sample. We consider the small fluorescent spheres as point emitters located at a depth z 0 inside the scattering medium. The diffuse emission from the fluorescent spheres is also described by the diffusion equation but with a Dirac delta function as the source term [50,68]. At steady-state, the solution is