# Mesoscopic light transport in nonlinear disordered media

Alfonso Nardi, Andrea Morandi, Romain Pierrat, Arthur Goetschy, Xuanchen Li, Frank Scheffold, Rachel Grange

► To cite this version:

Alfonso Nardi, Andrea Morandi, Romain Pierrat, Arthur Goetschy, Xuanchen Li, et al.. Mesoscopic light transport in nonlinear disordered media. Physical Review Research, 2025, 7 (2), pp.023323. 10.1103/ztyj-rc63. hal-05382977

HAL Id: hal-05382977

https://hal.science/hal-05382977v1

Submitted on 26 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Mesoscopic light transport in nonlinear disordered media

Alfonso Nardi,¹ Andrea Morandi,¹ Romain Pierrat,² Arthur Goetschy,²
Xuanchen Li,¹ Frank Scheffold,³ and Rachel Grange¹

¹ETH Zurich, Department of Physics, Institute for Quantum Electronics,
Optical Nanomaterial Group, 8093 Zurich, Switzerland

²ESPCI Paris, PSL University, CNRS, Institut Langevin, 75005 Paris, France

³Department of Physics, University of Fribourg, 1700 Fribourg, Switzerland

(Dated: September 10, 2024)

Nonlinear disordered media uniquely combine multiple scattering and second-harmonic genera-
tion. Here, we investigate the statistical properties of the nonlinear light generated within such
media. We report super-Rayleigh statistics of the second-harmonic speckle intensity, and demon-
strate that it is caused by the mesoscopic correlations arising in extreme scattering conditions. The
measured conductance is the lowest ever observed in an isotropically scattering 3D medium, with
applications in broadband second-harmonic generation, wavefront shaping in nonlinear disordered
media, and photonic computing.

When a coherent wave interacts with a disordered
medium, it generates a complex interference phenomenon
that results in the formation of speckle patterns. These
speckles, which are randomly distributed diffraction-
limited grains, are ubiquitous across waves of various
origins, including electromagnetic [1], acoustic [2], and
matter waves [3]. Despite their different physical ori-
gins, speckles exhibit universal statistical properties that
are referred to as Rayleigh statistics. This universality
is a consequence of the very general conditions under
which Rayleigh statistics emerge. Specifically, the only
requirement is that the field arises from the interference
of a large number of uncorrelated waves whose phases
are uniformly distributed over a 2π range [1, 4]. Non-
Rayleigh statistics can be achieved by using a spatial
light modulator to induce correlations between partial
waves and redistribute the intensity among the speckle
pattern's grains, while preserving ergodicity [5-7].

Conversely, observing deviations from Rayleigh statis-
tics in multiple-scattering media under conventional laser
illumination is significantly more challenging. Devia-
tions occur only when strong scattering introduces meso-
scopic correlations between partial waves, breaking er-
godicity [8–12]. These correlations are associated with
the existence of a finite number *g* of open transmission
channels in disordered media [13, 14]. When *g* (also
known as *dimensionless conductance*) is moderate, the
central limit theorem does not apply, causing the to-
tal intensity in the speckle pattern to fluctuate signifi-
cantly from one disorder configuration to another. Con-
sequently, the intensity in each speckle grain exhibits
super-Rayleigh statistics [9, 15]. This peculiar regime
has been reported in three-dimensional systems only in
nanowire mats [16] and, more weakly, in isotropic ZnO
scattering media [17].

While mesoscopic transport in linear scattering media
is well understood, the same cannot be said for nonlin-
ear disordered media, a class of materials that has gar-
nered significant interest recently. These materials are

composed of nanodomains that exhibit a second-order
susceptibility tensor. The nonlinear characteristics en-
able each nanodomain to exhibit electro-optic effects and
generate second-harmonic (SH) waves when illuminated
with a fundamental beam [18]. The interference of the
waves generated by the nanodomains results in efficient
emission of nonlinear light without stringent conditions
on the polarization and wavelength of the fundamental
light, contrary to bulk crystals [19]. This phenomenon
has been extensively studied in the framework of ran-
dom quasi-phase-matching [20, 21], where the disorder
is used to achieve efficient and broadband SH genera-
tion [22–26]. Researchers investigated the fundamental
properties of diffusion and weak localization in this class
of media [27–30], as well as the effect of scatterers dis-
placement [31, 32]. Notably, nonlinear disordered media
are also emerging as a prominent platform for photonics
processing. The characterization of the scattering tensor
that defines their nonlinear input-output response allows
the use of nonlinear disordered media for encryption and
as all-optical logic gates [33]. In addition, these media en-
ables the implementation of large-scale nonlinear optical
operators for photonic machine learning [34]. However,
simplifying the analysis by assuming Rayleigh statistics,
these applications might overlook the richer physics aris-
ing from strong scattering effects.

In this Letter, we examine the statistical properties of SH light generated within a strongly scattering, nonlinear disordered medium. Specifically, we measure the intensity fluctuations of both the fundamental and SH light under various illumination conditions. The histograms of speckle intensity and total transmission show deviations from Rayleigh statistics, particularly evident for SH light. Since SH light can be generated throughout the medium and cannot be described by a linear transmission matrix, an open question is whether these observations can be interpreted in terms of an effective number *g* of open channels for SH light. Here, we demonstrate that our results align well with mesoscopic transport theory
<figure><img src="image_2.png" /><figcaption>FIG. 1. Experimental setup and typical speckle patterns. (a) $\lambda_F$ = 976 nm: pulsed laser; $\lambda_{SH}$ = 488 nm; L: lens; NLDM: nonlinear disordered medium; SF: spectral filter; CAM: camera. (b-c) Camera images of typical speckle patterns at $\lambda_F$ when the sample is illuminated with (b) out-of-focus and (c) in-focus fundamental light. (d) Camera image of a typical speckle pattern from SH light generated within the medium with in-focus fundamental illumination.</figcaption></figure>

based on a reduced $g$, and we explain it by analyzing the location of SH generation within the medium and the wavelength dependence of the transport mean free path of light. This outcome has relevant applications in broadband SH generation, wavefront shaping [35, 36], and photonic computing [34].

Our experimental setup is shown in Fig. 1a. A linearly polarized Ti:Sa pulsed laser (wavelength $\lambda_F$ = 976 nm, 80 MHz repetition rate, 100 fs pulse duration) is focused by an aspheric lens (NA = 0.5) into a nonlinear disordered medium (disordered assembly of $LiNbO_3$ nanoparticles of thickness $L \approx 10$ µm, see Supplemental Material, Sec. S1 for detailed information [37]). Both the lens and the disordered sample are mounted on motorized translation stages, to accurately control the position of the focal plane, as well as the position of the sample in the plane perpendicular to the beam propagation direction. By moving the focal plane, we can control the size of the beam at the input facet of the disordered medium, a crucial parameter for the analysis of mesoscopic transport. Changing the position of the sample allows the light traveling through the medium to interact with different realizations of disorder. The light scattered and upconverted by the nonlinear disordered medium is collected by an objective (NA = 0.75), and a tube lens (focal length 200 mm) is used to image the output facet of the medium onto a scientific CMOS camera. Finally, we use a linear polarizer to select a single polarization state, and a spectral filter (bandpass region 360 – 580 nm) to

remove the fundamental light when measuring the SH signal ($\lambda_{SH} = \lambda_F/2 = 488$ nm). Typical speckle patterns for different illumination conditions are shown in Fig. 1b-d. Figure 1b displays a camera image of the fundamental light when the focal plane is far from the input facet of the medium, i.e., when the illuminating beam size is large. Figures 1c-d report typical speckle patterns for the fundamental light (Fig. 1c) and the SH light generated within the nonlinear disordered medium (Fig. 1d) when the focal plane of the fundamental illumination coincides with the input facet of the sample.

In the following, we rigorously characterize the statistical properties of the transmitted light after the propagation through the nonlinear disordered medium. In particular, we measure the intensity fluctuations, which contains strong signature of mesoscopic transport, and are not affected by absorption [38]. For a given realization of disorder, we use the camera to record the intensity distribution at the output facet of the disordered sample. From the measured speckle pattern, we obtain the transmission coefficients $T_{ab}$, which relate a speckle spot $b$ with an incoming wave $a$ (kept fixed for each statistical dataset). Summing over the output modes, we obtain the total transmission as $T_a = \sum_b T_{ab}$. To collect a statistically meaningful set of measurements, we measure the intensity distributions of $10^4$ different realizations of disorder, obtained by moving the sample in the plane perpendicular to the propagation direction of the illuminating light. From $T_{ab}$ and $T_a$, we then extract the most relevant quantities for our analysis, which are the normalized speckle intensity $s_{ab} = T_{ab}/\langle T_{ab} \rangle$ and the normalized total transmission $s_a = T_a/\langle T_a \rangle$, where the angle brackets denote the average over the ensemble of random configurations (see Supplemental Material, Sec. S2 for details about the data analysis [37]).

The results of the analysis for various illumination conditions are reported in Fig. 2. We first examine the case of fundamental light with out-of-focus illumination (Fig. 2, gray). The large size of the input beam ensures that we are in a regime of large number of open channels. In this case, the normalized speckle intensity obeys the Rayleigh law, meaning that the probability distribution $P(s_{ab})$ is a negative exponential (gray solid line in Fig. 2a). Conversely, the probability distribution of the normalized total transmission $P(s_a)$ follows a Gaussian distribution (Fig. 2b, gray). This is a consequence of the central limit theorem, as the numerous scattering paths that are summed to obtain the total transmission are uncorrelated.

Reducing the size of the illuminating beam by aligning the focal plane with the input facet of the medium (referred to as *in-focus*) increases the probability of two paths crossing during light propagation inside a disordered medium. When two paths cross, there is a small but non-zero probability that they will become correlated, giving rise to long-range correlations between dis-
<figure><img src="image_3.png" /></figure>

FIG. 2. Intensity fluctuations statistics of fundamental light with out-of-focus (OF Fund.) and in-focus (IF Fund.) illumination, and of SH light generated within the medium with in-focus fundamental illumination (IF SH). (a) Histograms of normalized speckle intensity. Solid lines represent a negative exponential $P(s_{ab}) = -\exp(s_{ab})$ (gray), and the mesoscopic transport model for normalized speckle intensity $P_g(s_{ab})$ with conductance values $g = 13$ (red) and $g = 6.5$ (blue). (b) Histograms of normalized total transmission. Solid lines indicate a Gaussian fit (gray), and the mesoscopic transport model for normalized total transmission $P_g(s_a)$, using the same conductance values as in (a), $g = 13$ (red) and $g = 6.5$ (blue).

<figure><img src="image_4.png" /></figure>

FIG. 3. Intensity fluctuations statistics of SH light generated within the nonlinear disordered medium for varying distance $\Delta z$ between the focal plane of the input lens and the entry facet of the nonlinear disordered medium ($\Delta z = \{-20~\mu\text{m}, -10~\mu\text{m}, 0\}$). (a) Histograms of normalized speckle intensity. Solid lines represent a negative exponential $P(s_{ab}) = -\exp(s_{ab})$ (gray), and the mesoscopic transport model $P_g(s_{ab})$, with conductance values $g = 16$ (violet), $g = 9.5$ (blue) and $g = 6.5$ (green). (b) Histograms of normalized total transmission. Solid lines are the mesoscopic probability density functions $P_g(s_a)$ with the same conductance values as in (a), $g = 16$ (violet), $g = 9.5$ (blue), $g = 6.5$ (green).

tant speckle grains [10, 39]. This, in turn, lowers the dimensionless conductance $g$, which quantifies the number of open channels, i.e., the number of channels that carry the majority of the light [13, 14, 40]. The introduced correlations change the statistical properties of the scattered light, enhancing in particular the intensity fluctuations [41]. Mesoscopic transport theory predicts modified probability distributions, referred to as $P_g(s_{ab})$ and $P_g(s_a)$, containing the dimensionless conductance $g$ as a single parameter (see Ref. [15] and Supplemental Material, Sec. S3 [37]).

The normalized speckle intensity and total transmission histograms for fundamental light with in-focus illumination is shown in red in Fig. 2a and b, respectively. Considering the speckle intensity, the slight deviation from Rayleigh distribution at large $s_{ab}$ values is consistent with previously reported observations in isotropic scattering materials [17]. Conversely, the total transmission histogram shows a clear deviation from the Gaussian distribution predicted by the uncorrelated wave model. Fitting the measured data with the mesoscopic transport model of Ref. [15] demonstrates good agreement (Fig. 2, red), resulting in a conductance value of $g = 13$.

With the same in-focus position of the input lens, we

measure the intensity statistics of the SH light generated within the nonlinear disordered medium. The histogram of the normalized speckle intensity $P(s_{ab})$, shown in blue in Fig. 2a, exhibits a more pronounced deviation from the Rayleigh distribution. Consequently, fitting the data with the mesoscopic model $P_g(s_{ab})$ yields a smaller conductance ($g = 6.5$). The normalized total transmission also demonstrates increased variance (Fig. 2b, in blue). Using the same conductance value of $g = 6.5$, derived from the fit of the speckle intensity histogram, we achieve a strong agreement between the theoretical model $P_g(s_a)$ and the experimental data. This is particularly evident for large $s_a$ values, where the data points closely follows a negative exponential of the form $P(s_a) \propto e^{-gs_a}$.

To confirm the agreement with the mesoscopic transport model, we characterized the intensity fluctuations of the SH light generated within the medium at different input lens positions ($L_1$ in Fig. 1), thereby varying the fundamental spot size at the entrance of the nonlinear medium. The results are shown in Fig. 3. We present histograms for three illumination conditions, based on the distance $\Delta z$ of the focal plane to the input facet of the medium: $\Delta z = -20~\mu\text{m}$, $\Delta z = -10~\mu\text{m}$, $\Delta z = 0$. As pre-
dicted by mesoscopic transport theory [39], the deviation of the speckle intensity histograms (Fig. 3a) from a negative exponential (gray solid line) becomes more evident as the focal plane of the illumination approaches the input facet of the sample. Similarly, in the histograms of the normalized total transmission (Fig. 2b), the exponential decay present at large $s_a$ becomes more pronounced for shorter $\Delta z$. The fitted probability distributions (violet, blue and green solid lines in Fig. 3) show excellent agreement for both the speckle intensity and the total transmission histograms, using the same conductance values for the expressions of $P_g(s_{ab})$ and $P_g(s_a)$ ($g = 16$ for $\Delta z = -20$ µm, $g = 9.5$ for $\Delta z = -10$ µm and $g = 6.5$ for $\Delta z = 0$).

The reduced conductance measured with the SH light is explained by the shorter transport mean free path $\ell_t$ associated to the SH wavelength ($\ell_t \approx 160$ nm at 450 nm, compared to $\ell_t \approx 700$ nm at 950 nm; see Supplemental Material, Sec. S4 for the experimental characterization of $\ell_t$ [37]). The transport mean free path is relevant for the characterization of the mesoscopic effects, because it limits the minimum conductance to a value $g_{min}(\lambda) \approx (2\pi\ell_t/\lambda)^2$ [39, 42]. Considering the measured $\ell_t$, we obtain $g_{min}(950$ nm) $\approx 11.3$ and $g_{min}(450$ nm) $\approx 2.6$ [37]. Notably, by accounting for the beam width dependence of $g$ with a finite SH illumination width of 600 nm, consistent with our experimental parameters, we estimate a conductance of $g = 6.5$ (see Supplemental Material, Sec. S4 [37]).

The question remains on how the statistics can be so accurately described with mesoscopic transport theory, given that the generation of the SH light throughout the medium prevents the definition of a linear transmission matrix [33]. The explanation lies in the fact that, due to the tight focusing, the SH generation efficiency is maximized at the entrance of the medium. Indeed, theoretical calculations show that the SH signal is mostly generated at $z < 500$ nm (see Fig. 4a, on the left; theoretical and simulation details in Supplemental Material, Sec. S5 [37]). This means that, given that the nanoparticles' sizes are distributed between 100 nm and 400 nm [37], the SH light is generated within the very first layers, and then undergoes linear scattering for the rest of the thickness $L \approx 10$ µm of the material (Fig. 4a, on the right). Therefore, under tightly focused illumination, we can assume that the system behaves similarly to a linear scattering medium. To prove it, we compare the statistical properties of SH light generated within the nonlinear medium, and by an external 200 µm thick BaB$_2$O$_4$ crystal (BBO, experimental setup depicted in Fig. 4b). Although the light generated within the nonlinear disordered medium (Fig. 4a, right) penetrates slightly deeper than the light generated by the BBO (Fig. 4c, right), the sources of both concentrate at depth $z \ll L \approx 10$ µm (Fig. 4a,c, left). In addition, the measured histograms of $s_{ab}$ for SH light from both the nonlinear medium and the

<figure><img src="image_5.png" /><figcaption>FIG. 4. Comparison of SH light generated within the nonlinear disordered medium (NLDM), or externally, via a BBO crystal. (a, c) Calculated mean intensity profiles within the disordered medium of the sources (left) and the resulting diffuse halos (right) for SH light generated (a) within the NLDM or (c) with a BBO crystal [37]. (b) Experimental setup. Light from a pulsed laser is upconverted with a BBO crystal ($\lambda_F = 976$ nm, $\lambda_{SH} = 488$ nm). A dichroic mirror (DM) filters out the fundamental light, while the SH is focused by the lens $L_1$ onto the NLDM. The scattered light is collected by an objective $L_2$ and measured with a camera (CAM). (d) Histograms of normalized speckle intensity $P(s_{ab})$ for SH light generated with a BBO crystal (cross, yellow) or by the NLDM (circles, blue).</figcaption></figure>

BBO (shown in Fig. 4d in yellow and blue, respectively) closely match. The similarity is expected, as the conductance is limited by the transport mean free path and the wavelength, which are identical in the two cases. The excellent agreement between the two histograms confirms that the conductance values measured for SH light generated in the medium are determined by the scattering properties of the sample and can be accurately described by mesoscopic transport theory.

In summary, we have experimentally measured the intensity fluctuations of fundamental and SH light generated by a nonlinear disordered medium. We observed distinct deviations from the predictions of the uncorrelated wave model in the histograms of SH speckle intensity and total transmission. The measured deviations are well described by mesoscopic transport model, marking the first observation of mesoscopic transport in a nonlinear disordered medium to our knowledge. Moreover, the measured conductance of $g = 6.5$ is the lowest ever reported in an isotropically scattering 3D medium. We motivated the description of the nonlinear input-output response in terms of effective number $g$ of open channels by observing that the SH light is mainly generated at the
input of the nonlinear disordered medium. Subsequently, the SH light undergoes linear scattering, allowing us to approximate the nonlinear disordered sample as a linear scattering system. We confirmed this reasoning by comparing the speckle intensity of SH light generated either within the medium, or externally with a BBO crystal, obtaining very similar fluctuations statistics. Our results show SH speckles with intensities more than 20 times stronger than the average, occurring with a probability two orders of magnitude higher than that of uncorrelated light. The increased probability of generating higher SH intensity condensed in few speckles finds useful applications in broadband SH generation from nonlinear disordered media. Additionally, achieving a low number of open transmission modes has also significant implications for wavefront shaping, as the long-range correlations enhance the achievable control [35, 36]. The deviations from the uncorrelated wave model indicate that, when the fundamental light is strongly focused to achieve a higher SH signal, the wave propagation through the complex nonlinear medium cannot be described by a fully random nonlinear input-output response [33]. It is therefore essential to consider this effect for photonic computing applications that aim to exploit randomness to achieve large-scale nonlinear operations [34].

We acknowledge Maria Teresa Buscaglia and Vincenzo Buscaglia from the Institute of Condensed Matter Chemistry and Technologies for Energy (National Research Council, via De Marini 6, 16149 Genoa, Italy) for the synthesis of the LiNbO₃ nanocubes. This research was supported by the Swiss National Science Foundation SNSF via Consolidator Grant 213713 (A.N., A.M., R.G.) and via project grant 188494 (F.S.). R.P. and A.G. acknowledge the support of the program Investissements d'Avenir launched by the French Government. A.N. and A.M. contributed equally to this work.

[1] J. W. Goodman, *Speckle phenomena in optics: theory and applications*, 2nd ed. (The International Society for Optical Engineering, Bellingham, 2020).

[2] R. Wagner, S. Smith, J. Sandrik, and H. Lopez, Statistics of Speckle in Ultrasound B-Scans, IEEE Transactions on Sonics and Ultrasonics **30**, 156 (1983).

[3] R. G. Dall, S. S. Hodgman, A. G. Manning, M. T. Johnson, K. G. H. Baldwin, and A. G. Truscott, Observation of atomic speckle and Hanbury Brown–Twiss correlations in guided matter waves, Nature Communications **2**, 291 (2011).

[4] R. Carminati and J. C. Schotland, *Principles of Scattering and Transport of Light* (Cambridge University Press, Cambridge, 2021).

[5] Y. Bromberg and H. Cao, Generating Non-Rayleigh Speckles with Tailored Intensity Statistics, Physical Review Letters **112**, 213904 (2014).

[6] R. Liu, B. Qing, S. Zhao, P. Zhang, H. Gao, S. Chen, and F. Li, Generation of Non-Rayleigh Nondiffracting Speckles, Physical Review Letters **127**, 180601 (2021).

[7] S. Han, N. Bender, and H. Cao, Tailoring 3D Speckle Statistics, Physical Review Letters **130**, 093802 (2023).

[8] S. Feng, C. Kane, P. A. Lee, and A. D. Stone, Correlations and fluctuations of coherent wave transmission through disordered media, Phys. Rev. Lett. **61**, 834 (1988).

[9] N. Garcia and A. Z. Genack, Crossover to strong intensity correlation for microwave radiation in random media, Phys. Rev. Lett. **63**, 1678 (1989).

[10] M. P. van Albada, J. F. de Boer, and A. Lagendijk, Observation of long-range intensity correlation in the transport of coherent light through a random medium, Phys. Rev. Lett. **64**, 2787 (1990).

[11] F. Scheffold and G. Maret, Universal conductance fluctuations of light, Phys. Rev. Lett. **81**, 5800 (1998).

[12] S. Balog, P. Zakharov, F. Scheffold, and S. E. Skipetrov, Photocount statistics in mesoscopic optics, Phys. Rev. Lett. **97**, 103901 (2006).

[13] O. N. Dorokhov, On the coexistence of localized and extended electronic states in the metallic phase, Solid State Communications **51**, 381 (1984).

[14] Y. Imry, Active Transmission Channels and Universal Conductance Fluctuations, Europhysics Letters **1**, 249 (1986).

[15] T. M. Nieuwenhuizen and M. C. W. van Rossum, Intensity Distributions of Waves Transmitted through a Multiple Scattering Medium, Physical Review Letters **74**, 2674 (1995).

[16] T. Strudley, T. Zehender, C. Blejean, E. P. A. M. Bakkers, and O. L. Muskens, Mesoscopic light transport by very strong collective multiple scattering in nanowire mats, Nature Photonics **7**, 413 (2013).

[17] T. Strudley, D. Akbulut, W. L. Vos, A. Lagendijk, A. P. Mosk, and O. L. Muskens, Observation of intensity statistics of light transmitted through 3D random media, Optics Letters **39**, 6347 (2014).

[18] F. Timpu, J. Sendra, C. Renaut, L. Lang, M. Timofeeva, M. T. Buscaglia, V. Buscaglia, and R. Grange, Lithium Niobate Nanocubes as Linear and Nonlinear Ultraviolet Mie Resonators, ACS Photonics **6**, 545 (2019).

[19] R. W. Boyd, *Nonlinear Optics* (Academic Press, Burlington, 2008).

[20] M. Baudrier-Raybaut, R. Haïdar, P. Kupecek, P. Lemasson, and E. Rosencher, Random quasi-phase-matching in bulk polycrystalline isotropic nonlinear materials, Nature **432**, 374 (2004).

[21] E. V. Makeev and S. E. Skipetrov, Second harmonic generation in suspensions of spherical particles, Optics Communications **224**, 139 (2003).

[22] R. Fischer, S. M. Saltiel, D. N. Neshev, W. Krolikowski, and Y. S. Kivshar, Broadband femtosecond frequency doubling in random media, Applied Physics Letters **89**, 191105 (2006).

[23] Y. Qiao, F. Ye, Y. Zheng, and X. Chen, Cavity-enhanced second-harmonic generation in strongly scattering nonlinear media, Physical Review A **99**, 043844 (2019).

[24] R. Savo, A. Morandi, J. S. Müller, F. Kaufmann, F. Timpu, M. Reig Escalé, M. Zanini, L. Isa, and R. Grange, Broadband Mie driven random quasi-phase-matching, Nature Photonics **14**, 740 (2020).

[25] J. S. Müller, A. Morandi, R. Grange, and R. Savo, Modeling of Random Quasi-Phase-Matching in Birefringent
Disordered Media, Physical Review Applied **15**, 064070 (2021).

[26] A. Morandi, R. Savo, J. S. Müller, S. Reichen, and R. Grange, Multiple Scattering and Random Quasi-Phase-Matching in Disordered Assemblies of $LiNbO_3$ Nanocubes, ACS Photonics **9**, 1882 (2022).

[27] V. M. Agranovich and V. E. Kravtsov, Effects of weak localization of photons in nonlinear optics: Second harmonic generation, Physics Letters A **131**, 378 (1988).

[28] K. M. Yoo, S. Lee, Y. Takiguchi, and R. R. Alfano, Search for the effect of weak photon localization in second-harmonic waves generated in a disordered anisotropic nonlinear medium, Optics Letters **14**, 800 (1989).

[29] S. Faez, P. M. Johnson, D. A. Mazurenko, and A. Lagendijk, Experimental observation of second-harmonic generation and diffusion inside random media, JOSA B **26**, 235 (2009).

[30] C. I. Valencia and E. R. Méndez, Weak localization effects in the second-harmonic light scattered by random systems of particles, Optics Communications **282**, 1706 (2009).

[31] R. Samanta and S. Mujumdar, Intensity-dependent speckle contrast of second harmonic light in a nonlinear disordered medium, Applied Optics **59**, 11266 (2020).

[32] R. Samanta, R. Pierrat, R. Carminati, and S. Mujumdar, Speckle Decorrelation in Fundamental and Second-Harmonic Light Scattered from Nonlinear Disorder, Physical Review Applied **18**, 054047 (2022).

[33] J. Moon, Y.-C. Cho, S. Kang, M. Jang, and W. Choi, Measuring the scattering tensor of a disordered nonlinear medium, Nature Physics **19**, 1709 (2023).

[34] H. Wang, J. Hu, A. Morandi, A. Nardi, F. Xia, X. Li, R. Savo, Q. Liu, R. Grange, and S. Gigan, Large-scale photonic computing with nonlinear disordered media, Nature Computational Science **4**, 429 (2024).

[35] C. W. Hsu, S. F. Liew, A. Goetschy, H. Cao, and A. D. Stone, Correlation-enhanced control of wave focusing in disordered media, Nature Physics **13**, 497 (2017).

[36] N. Bender, A. Yamilov, A. Goetschy, H. Yılmaz, C. W. Hsu, and H. Cao, Depth-targeted energy delivery deep inside scattering media, Nature Physics **18**, 309 (2022).

[37] See Supplemental Material at [URL will be inserted by publisher], which included Refs. [42–46], for the fabrication procedure of the nonlinear disordered sample, for the details of the data analysis, for the characterization of the transport mean free path and the minimum conductance, and for the derivation of the equations that govern the propagation of the linear and SH mean intensities through the disordered slab.

[38] A. A. Chabanov, M. Stoytchev, and A. Z. Genack, Statistical signatures of photon localization, Nature **404**, 850 (2000).

[39] F. Scheffold, W. Härtl, G. Maret, and E. Matijević, Observation of long-range correlations in temporal intensity fluctuations of light, Physical Review B **56**, 10942 (1997).

[40] C. W. J. Beenakker, Random-matrix theory of quantum transport, Reviews of Modern Physics **69**, 731 (1997).

[41] M. C. W. van Rossum and T. M. Nieuwenhuizen, Multiple scattering of classical waves: microscopy, mesoscopy, and diffusion, Rev. Mod. Phys. **71**, 313 (1999).

[42] A. García-Martín, F. Scheffold, M. Nieto-Vesperinas, and J. J. Sáenz, Finite-size effects on intensity correlations in random media, Phys. Rev. Lett. **88**, 143901 (2002).

[43] M. Reufer, L. F. Rojas-Ochoa, S. Eiden, J. J. Sáenz, and F. Scheffold, Transport of light in amorphous photonic materials, Applied Physics Letters **91** (2007).

[44] L. Schertel, L. Siedentop, J.-M. Meijer, P. Keim, C. M. Aegerter, G. J. Aubry, and G. Maret, The structural colors of photonic glasses, Advanced Optical Materials **7**, 1900442 (2019).

[45] D. Contini, F. Martelli, and G. Zaccanti, Photon migration through a turbid slab described by a model based on diffusion approximation. i. theory, Applied optics **36**, 4587 (1997).

[46] V. A. Markel, Introduction to the maxwell garnett approximation: tutorial, J. Opt. Soc. Am. A **33**, 1244 (2016).
Supplemental Material:
Mesoscopic light transport in nonlinear disordered media

Alfonso Nardi,¹, * Andrea Morandi,¹, * Romain Pierrat,²

Arthur Goetschy,² Xuanchen Li,¹ Frank Scheffold,³ and Rachel Grange¹

¹ETH Zurich, Department of Physics, Institute for Quantum Electronics, Optical Nanomaterial Group, 8093 Zurich, Switzerland

²ESPCI Paris, PSL University, CNRS, Institut Langevin, 75005 Paris, France

³Department of Physics, University of Fribourg, 1700 Fribourg, Switzerland

(Dated: September 10, 2024)

This document provides supplementary information for the article "Mesoscopic Light Transport in Nonlinear Disordered Media". In Sec. S1, we present the procedure for fabricating the nonlinear disordered medium. Section S2 reports the detailed data analysis, and Sec. S3 presents the theoretical expressions for the probability density functions according to mesoscopic transport theory [1]. In Sec. S4 we report the measurements of the sample transmittance, yielding the experimental values of the transport mean free path and the minimum conductance. Finally, in Sec. S5, we derive the equations that govern the propagation of the linear and SH mean intensities through the disordered slab.

S1. FABRICATION PROCEDURE

The nonlinear disordered medium used is a slab, assembled starting from crystalline LiNbO₃ nanoparticles, pro-
duced by solvothermal synthesis [2]. The precursor oxides, Nb2O5 (HC Starck, 99.92%) and LiOH (Aldrich, 98%),
are dispersed in a mixture of ethylene glycol and distilled water. Following ultrasonication, the suspension is trans-
ferred into a polytetrafluoroethylene-coated stainless steel acid digestion bomb (model PA4748, volume 120 ml, Parr
Instrument Company) and treated hydrothermally at 250°C for 70 hours. The reaction product is then washed
with water via centrifugation. This chemical synthesis method allows precise control over the size distribution of
the nanoparticles, which range from 100 to 400 nm in size, have a linear refractive index of approximately 2.3, and
exhibit negligible absorption at visible wavelengths. These nanoparticles feature a non-centrosymmetric hexagonal
R3c crystal structure that enables second-harmonic (SH) generation. To form the slabs, we deposit the nanoparticles
by drop deposition and allow solvent evaporation. An aqueous suspension (1wt%) of the LiNbO₃ nanoparticles is
mixed with polyvinyl alcohol in a 60:1 ratio and deposited onto a glass substrate framed with hydrophilic tape. The
sample is placed on a horizontal substrate holder and maintained at 0 °C for 24 hours. An optical microscope image
of the resulting slab is reported in Fig. S1a. The white opaque color is typical of multiple scattering sample that have
a thickness much larger than the transport mean free path. The porous structure composed by randomly oriented
nanoparticles and air gaps is shown in Fig. S1b, measured via scanning electron microscopy (SEM). By controlling
the amount of deposited nanoparticles solution, we can tune the thickness of the slab. We measured the thickness
of the sample used for the experiments in the main text by profilometry. The resulting thickness map is shown in
Fig. S2a.

<figure><img src="image_6.png" /><figcaption>FIG. S1. Nonlinear disordered sample. (a) Optical image (in reflection configuration) of the disordered sample. (b) Scanning electronic microscope image of the slab, showing a typical arrangement of the nanoparticles.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. S2. Thickness map of the sample, measured with a profilometer. (a) Thickness map of the entire sample. (b) Thickness map of the area of 300 × 300 µm², used for the mesoscopic transport analysis. (c) Thickness profile along the blue dashed line.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. S3. Total transmission fluctuations over different positions of the sample, resulting from the thicknesses gradient and fluctuations of the laser power. (a) Fluctuation over the entire dataset, and (b) close up on 1000 steps. The orange line is the result of a moving average over 100 steps.</figcaption></figure>

## S2. DATA ANALYSIS

We collected the data by measuring the speckle patterns ($\tilde{T}_{ab,i}$, where b are the output modes and a the fixed input modes) at $i = 1...N = 100 \times 100$ different sample positions. For each measurement, we moved the sample by 3 µm using a motorized stage. The thickness of the total scanned area of 300 × 300 µm², measured by profilometry, is shown in Fig. S2b. The thickness along the blue dashed line in Fig. S2b is presented in Fig. S2c. Due to the inhomogeneous thickness and the laser fluctuations, the total transmission ($TT_i$, i.e., the sum of all the intensities on the camera pixels) exhibited periodic fluctuations, as shown in Fig. S3. To account for these fluctuations, we used a moving average with a window size of 100 steps to smooth the total transmission. This filtered total transmission $TT_i^{\text{filt.}}$ is shown in orange in Fig. S3. We normalized each speckle pattern by the corresponding value of filtered total transmission (similar approach to Ref. [3])

$$T_{ab,i} = \tilde{T}_{ab,i} \frac{\langle TT_i \rangle}{TT_i^{\text{filt.}}}, \qquad (S1)$$

where the angle brackets stands for the average over the N realizations of disorder. We then extracted the ensemble average $\langle T_{ab} \rangle$ of the speckle patterns, as

$$\langle T_{ab} \rangle = \frac{1}{N} \sum_{i=1}^{N} T_{ab,i}. \qquad (S2)$$

For the rest of the analysis, we considered only the central area (average intensity larger than 75% of the maximum) of the measured images to minimize the effect of the background. For each collected camera image, we subtracted the background and normalized the speckle patterns by the average intensity:

$$s_{ab,i} = \frac{T_{ab,i}}{\langle T_{ab} \rangle}. \qquad (S3)$$
We then collected all N realizations into a single vector, to obtain the histogram of $s_{ab}$. Finally, by summing all the output modes, we derived the normalized total transmission

$$s_{a,i} = \sum_{b} s_{ab,i} \qquad (S4)$$

and collected the N realizations into the same vector $s_a$.

## S3. THEORETICAL EXPRESSIONS OF PROBABILITY DENSITY FUNCTIONS

For clarity, we present the expressions for $P_g(s_a)$ and $P_g(s_{ab})$ used to fit the histograms, as derived in Ref. [1]. The probability density function $P_g(s_a)$ for the normalized total transmission is

$$P_g(s_a) = \int_{-i\infty}^{i\infty} \frac{dx}{2\pi i} \exp[xs_a - \Psi(x)]. \qquad (S5)$$

For broad Gaussian illumination (waist $w$ larger than the sample thickness $L$), it was proved that the function $\Psi(x)$ takes the form

$$\Psi(x) = g \int_{0}^{1} \frac{dy}{y} \left[ \log \left( \sqrt{1 + \frac{xy}{g}} + \sqrt{\frac{xy}{g}} \right) \right]^2, \qquad (S6)$$

where $g = k^2 w^2 \ell_t / 3L$. We note that this definition of $g$, specific for Gaussian beam, leads to the following relation between the second moment of the normalized total transmission and the conductance $g$:

$$\langle s_a^2 \rangle = 1 + \frac{1}{3g}. \qquad (S7)$$

For a Gaussian beam with a waist significantly smaller than the slab thickness, an analytical expression for $P_g(s_a)$ is not available. However, as argued in Ref. [4], the first cumulants are expected to remain nearly unchanged, provided that the parameter $g$ is properly computed for this specific geometry. This is why we use Eq. (S6) to fit our experimental data, with $g$ as the only fitting parameter.

Furthermore, the probability density function $P_g(s_{ab})$ for the normalized speckle intensity is given by

$$P_g(s_{ab}) = \int_0^\infty ds_a P_g(s_a) \frac{e^{-s_{ab}/s_a}}{s_a}. \qquad (S8)$$

The latter expression reveals that the speckle intensity distribution differs from the Rayleigh law, $P(s_{ab}) = e^{-s_{ab}}$, only if the normalized total transmission $s_a$ significantly deviates from its mean $\langle s_a \rangle = 1$.

## S4. TRANSPORT MEAN FREE PATH

We characterized the transport mean free path $\ell_t$ by measuring the scaling of the total transmission with the sample thickness. The relation between thickness and transport mean free path in absence of absorption follows [5, 6]

$$T = \frac{\ell_t + z_0}{L + 2z_0} \qquad (S9)$$

where $T$ is the transmittance, $L$ is the local thickness of the sample, and $z_0$ the extrapolation length. The expression of $z_0$ is

$$z_0 = \frac{2}{3} \left( \frac{1 + R_i}{1 - R_i} \right) \ell_t = \beta \ell_t. \qquad (S10)$$

with $R_i$ the internal reflectivity. The internal reflectivity can be derived from the effective refractive index of the slab $n_{\text{eff}}$. In practice, we use the following approximate expression which is a polynomial fit of the analytical expression [7]

$$\frac{1 + R_i}{1 - R_i} = 504.332889 - 2641.0021 n_{\text{eff}} + 5923.699064 n_{\text{eff}}^2 - 7376.355814 n_{\text{eff}}^3 + \\ + 5507.53041 n_{\text{eff}}^4 - 2463.357945 n_{\text{eff}}^5 + 610.965647 n_{\text{eff}}^6 - 64.8047 n_{\text{eff}}^7. \qquad (S11)$$
<figure><img src="image_9.png" /><figcaption>FIG. S4. Thickness and inverse transmittance maps of the sample used for the characterization of the transport mean free path. (a) Thickness map measured with a profilometer. (b, c) Inverse transmittance measured at (b) 950 nm and (c) 450 nm. The circles represent the area used for the analysis of $\ell_t$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. S5. Experimental setup for the characterization of the transmittance. (a) For the measurement at 950 nm, the laser light is directly focused on the disordered medium. (b) Light at 450 nm is obtained by upconverting the laser with a BBO crystal. Then, the fundamental light is filtered with a dichroic mirror, and the SH is focused on the nonlinear disordered sample. Finally, the scattered light is collected in both cases by an integrating sphere, and directed towards a detector. L: lens; DM: dichroic mirror; NLDM: nonlinear disordered medium; IS: integrating sphere; DET: detector.</figcaption></figure>

We obtain the effective refractive index from the Maxwell-Garnett mixing rule [8]

$$n_{\text{eff}}^2 = \frac{1 + \frac{1+2\text{ff}}{3}(n_{\text{LNO}}^2 - 1)}{1 + \frac{1-\text{ff}}{3}(n_{\text{LNO}}^2 - 1)}, \qquad (\text{S12})$$

where $n_{\text{LNO}}$ is the average refractive index of bulk lithium niobate, and ff is the filling fraction of the nanocubes in the composite. We calculate the filling fraction by comparing the mass $m$ of the sample, and its volume $V$, estimated from the profilometer map (Fig. S4a), using the formula

$$\text{ff} = \frac{m}{V\rho_{\text{LNO}}} \approx 0.55, \qquad (\text{S13})$$

with $\rho_{\text{LNO}}$ the density of lithium niobate. Combining Eqs. (S10-S13) we obtain the value of $\beta$. Given $\beta$, we recast Eq. (S9) to highlight the linear relation between the inverse transmission and the transport mean free path

$$\frac{1}{T} = \frac{1}{\ell_t + z_0} L + \frac{2z_0}{\ell_t + z_0} = \frac{1}{\ell_t(1 + \beta)} L + \frac{2\beta}{1 + \beta}, \qquad (\text{S14})$$

which has $\ell_t$ as the only free parameter. Therefore, we can derive $\ell_t$ from the slope of the relation between the local thickness of the material $L$ and the inverse transmittance.

To this aim, we experimentally measured the thickness map by profilometry (see Fig. S4a), and the transmittance of the disordered sample with the experimental setup sketched in Fig. S5. For the measurement of the transmittance at the fundamental wavelength ($\lambda_F = 950$ nm), the laser light is directly focused on the nonlinear disordered medium (Fig. S5a). For the measurement at the SH wavelength ($\lambda_{\text{SH}} = 450$ nm), we first needed to frequency convert the fundamental beam with a BBO crystal. We then filtered the fundamental light with a dichroic mirror, and we focused the SH beam on the sample (Fig. S5b). In both cases, the scattered light is collected with an integrating sphere, and then guided with a multimode fiber to a powermeter. We subtracted the background intensity (due to stray light reaching the integrating sphere), and we normalized the measured transmittance values with the intensity transmitted through the bare glass coverslip. By moving the sample in the plane orthogonal to the illuminating beam, we obtained
<figure><img src="image_11.png" /><figcaption>FIG. S6. Measured relation between thickness and inverse transmittance at 450 nm and 950 nm. The dots and the errorbars represent the mean values and the standard deviation, respectively. From the slope of the linear fit (solid lines), we extract the transport mean free paths $\ell_t$ according to Eq. (S14).</figcaption></figure>

the 2D map of the sample inverse transmittance at the fundamental and SH wavelength (reported in Figs. S4b and c, respectively). By matching the profilometer and the inverse transmittance map, we obtain the relation between the thickness of the sample and the inverse transmittance (Fig. S6). The areas evaluated for the fit are highlighted with circles in Figs. S4b, c. Note that, due to the lower power available at the SH wavelength, the transmittance measurements at high thicknesses are strongly affected by noise, thus we did not include them in the analysis of the transport mean free path.

We fit the measured data with Eq. (S14), resulting in the following values of transport mean free paths:

* $\ell_t = 700$ nm at $\lambda = 950$ nm,

* $\ell_t = 160$ nm at $\lambda = 450$ nm.

The transport mean free path is relevant for the characterization of the mesoscopic effects, because it limits the minimum measurable conductance. In fact, the light incident on the slab does not contribute to the long-range correlations before it is scattered at least once, linking the intensity fluctuations to the scattering mean free path [9]. Combining Eq. (S7), valid for Gaussian illumination, with the heuristic expression of the intensity fluctuation $C_2 = \langle s_a^2 \rangle - 1$ established in Ref. [10] for small beam diameter, we obtain the following expression of the conductance $g$:

$$g = \frac{1}{3} \frac{1}{\langle s_a^2 \rangle - 1} \approx \frac{2}{3} \left(\frac{8}{9}\right)^2 \frac{[1 + (9/32) w/\ell_t]^2}{1 + (3/16) w/\ell_t} (k\ell_t)^2, \quad (S15)$$

where $w$ is the input beam width. For $w \to 0$ we then obtain the minimum measurable conductance $g_{min}$, which in our case is

* $g_{min}(950$ nm) $\approx 11.3$,

* $g_{min}(450$ nm) $\approx 2.6$.

It is noteworthy that, to obtain a conductance value of $g = 6.5$ in Eq. (S15), we should consider an input width of $w = 600$ nm, which corresponds well to the transverse width of the source intensity profile reported in the main text (Fig. 4a on the left).

## S5. TRANSPORT MODELS IN THE DIFFUSIVE REGIME

To compare the linear propagation of light at frequency $2\omega$ with the propagation of SH light generated within the disordered medium at the same frequency, we compute the mean intensity for both cases, using radiative transport models. The model for SH light, detailed in Ref. [11], has been shown to quantitatively match simulations of the microscopic wave equation that incorporate the SH generation process. Here we solve both propagation models in the
diffusive limit where the transport mean free path is much smaller than the medium thickness, using focused beams for illumination.

We consider a disordered slab translation-invariant in x and y directions and of thickness L in direction z. Scattering is supposed to be isotropic such that the scattering and transport mean-free paths are equals (i.e., $l_s = l_t$ and anisotropy factor $g = 0$) and there is no absorption. When this slab is illuminated by a Gaussian beam from $z < 0$ at a frequency $\omega$, the mean ballistic intensity at depth $z \ge 0$ and frequency $\omega$ takes the following expression:

$$I_b(\mathbf{R}, z, \omega) = I_0 \exp \left[ -\frac{2R^2}{w^2} - \frac{z}{l_s(\omega)} \right], \qquad (S16)$$

where $I_0$ is the intensity of the incident beam of waist $w$, and $\mathbf{R} = (x, y)$ is the transverse coordinate, with $R = \sqrt{x^2 + y^2}$. This ballistic intensity is a source term for the diffuse intensity $I_d$ which obeys the stationary diffusion equation

$$-\frac{l_s(\omega)^2}{3} \nabla^2 I_d(\mathbf{R}, z, \omega) = I_b(\mathbf{R}, z, \omega) \qquad (S17)$$

with the boundary conditions

$$I_d(\mathbf{R}, z = 0, \omega) - z_0(\omega) \frac{\partial I_d}{\partial z}(\mathbf{R}, z = 0, \omega) = 0, \qquad (S18)$$

$$I_d(\mathbf{R}, z = L, \omega) + z_0(\omega) \frac{\partial I_d}{\partial z}(\mathbf{R}, z = L, \omega) = 0, \qquad (S19)$$

where $z_0(\omega) = 2l_s(\omega)/3$ is the extrapolation length. This set of equations can be easily solved applying a Fourier transform with respect to $\mathbf{R}$. In particular, the Green function, solution of

$$-\left(\frac{\partial^2}{\partial z^2} - q^2\right) G_d(\mathbf{q}, z, z', \omega) = \delta(z - z') \qquad (S20)$$

with

$$G_d(\mathbf{q}, z = 0, z', \omega) - z_0(\omega) \frac{\partial G_d}{\partial z}(\mathbf{q}, z = 0, z', \omega) = 0, \qquad (S21)$$

$$G_d(\mathbf{q}, z = L, z', \omega) + z_0(\omega) \frac{\partial G_d}{\partial z}(\mathbf{q}, z = L, z', \omega) = 0, \qquad (S22)$$

is given by

$$G_d(\mathbf{q}, z, z', \omega) = \begin{cases} \frac{1}{q} \frac{[\sinh(qz) + qz_0(\omega) \cosh(qz)][\sinh\{q(L-z')\} + qz_0(\omega) \cosh\{q(L-z')\}]}{[1 + q^2z_0(\omega)^2] \sinh(qL) + 2qz_0(\omega) \cosh(qL)} & \text{for } z < z', \\ \frac{1}{q} \frac{[\sinh(qz') + qz_0(\omega) \cosh(qz')][\sinh\{q(L-z)\} + qz_0(\omega) \cosh\{q(L-z)\}]}{[1 + q^2z_0(\omega)^2] \sinh(qL) + 2qz_0(\omega) \cosh(qL)} & \text{for } z > z'. \end{cases} \qquad (S23)$$

From this expression, the diffuse intensity is given by

$$I_d(\mathbf{q}, z, \omega) = \frac{3}{l_s(\omega)^2} \int_0^L dz' G_d(\mathbf{q}, z, z', \omega) I_b(q, z', \omega), \qquad (S24)$$

where $I_b(\mathbf{q}, z', \omega)$ is the Fourier transform of $I_b(\mathbf{R}, z', \omega)$. We can write it explicitly as

$$I_d(\mathbf{q}, z, \omega) = I_0 \frac{3\pi w^2}{2l_s(\omega)^2} e^{-q^2 w^2 / 8} \bar{G}_d(\mathbf{q}, z, \omega), \qquad (S25)$$

where the integrated Green's function, $\bar{G}_d(\mathbf{q}, z, \omega) = \int_0^L dz' G_d(\mathbf{q}, z, z', \omega) e^{-z'/l_s(\omega)}$, is given by

$$\bar{G}_d(\mathbf{q}, z, \omega) = \frac{l_s(\omega)^2}{q^2 l_s(\omega)^2 - 1} \left[ e^{-z/l_s(\omega)} - \frac{e^{-L/l_s(\omega)}[2ql_s(\omega) \cosh(qz) + 3 \sinh(qz)] + 10ql_s(\omega) \cosh\{q(L-z)\} + 15 \sinh\{q(L-z)\}}{12ql_s(\omega) \cosh(qL) + [9 + 4q^2 l_s(\omega)^2] \sinh(qL)} \right]. \qquad (S26)$$
Taking the inverse Fourier transform, we obtain

$$I_d(\mathbf{R}, z, \omega) = I_0 \frac{3w^3}{4\ell_s(\omega)^2} \int_0^\infty dq e^{-q^2 w^2 / 8} \bar{G}_d(\mathbf{q}, z, \omega) q J_0(qR), \quad (S27)$$

where $J_0$ is the Bessel function of the first kind of order zero. Finally, the total mean intensity for linear propagation at frequency $\omega$ reads

$$I_{\text{lin}}(\mathbf{R}, z, \omega) = I_b(\mathbf{R}, z, \omega) + I_d(\mathbf{R}, z, \omega). \quad (S28)$$

The source term represented in the left part of Fig. 4c of the main text corresponds to $I_b(\mathbf{R}, z, 2\omega)$, while the total intensity profile shown in the right part is $I_{\text{lin}}(\mathbf{R}, z, 2\omega)$, evaluated at $x = 0$. Parameters of the calculation are chosen close to the experimental values: the waist $w \simeq 0.41 \ \mu\text{m}$ of the incident beam is such that its FWHM, $\sqrt{2\ln(2)}w$, is equal to $(\lambda_F/2)/2NA$ with $\lambda_F = 0.976 \ \mu\text{m}$ and NA = 0.5; the mean free path is $\ell_s(2\omega) = 0.2 \ \mu\text{m}$; the sample thickness is $L = 10 \ \mu\text{m}$.

Regarding SH light propagating at frequency $2\omega$, it is possible to show that the diffuse mean intensity still verifies a diffusion equation, with a source term that involves the square of the mean intensity of light propagating at frequency $\omega$. For isotropic scattering, we have [11]

$$-\frac{\ell_s(2\omega)}{3} \nabla^2 I_d(\mathbf{R}, z, 2\omega) = \mathcal{M}_0(\mathbf{R}, z, \omega), \quad (S29)$$

with

$$\mathcal{M}_0(\mathbf{R}, z, \omega) = \frac{\alpha}{4\pi} I_{\text{lin}}(\mathbf{R}, z, \omega)^2. \quad (S30)$$

Here $\alpha$ is a parameter proportional in particular to the SH susceptibility $\chi^{(2)}$, whose explicit value does not affect the profile of the SH light propagation. With the same boundary conditions as in the fundamental case, the solution is given by

$$I_d(\mathbf{R}, z, 2\omega) = \frac{3}{\ell_s(2\omega)} \int_0^L dz' \iint d\mathbf{R}' G_d(\mathbf{R} - \mathbf{R}', z, z', 2\omega) \mathcal{M}_0(\mathbf{R}', z', \omega), \quad (S31)$$

where $G_d(\mathbf{R}, z, z', \omega)$ is the inverse Fourier transform of the expression given in Eq. (S23). The source term shown in the left part of Fig. 4a of the main text corresponds to Eq. (S30) evaluated at $x = 0$, with $w \simeq 0.83 \ \mu\text{m}$ and $\ell_s(\omega) = 0.7 \ \mu\text{m}$. On the other hand, the right part of Fig. 4a shows the diffusive SH intensity given by Eq. (S31), computed with the values of the parameters mentioned above.

* These authors contributed equally to this work.

[1] T. M. Nieuwenhuizen and M. C. W. van Rossum, Intensity Distributions of Waves Transmitted through a Multiple Scattering Medium, Physical Review Letters 74, 2674 (1995).

[2] F. Timpu, J. Sendra, C. Renaut, L. Lang, M. Timofeeva, M. T. Buscaglia, V. Buscaglia, and R. Grange, Lithium Niobate Nanocubes as Linear and Nonlinear Ultraviolet Mie Resonators, ACS Photonics 6, 545 (2019).

[3] T. Strudley, T. Zehender, C. Blejean, E. P. A. M. Bakkers, and O. L. Muskens, Mesoscopic light transport by very strong collective multiple scattering in nanowire mats, Nature Photonics 7, 413 (2013).

[4] M. C. W. van Rossum and T. M. Nieuwenhuizen, Rev. Mod. Phys. 71, 313 (1999).

[5] M. Reufer, L. F. Rojas-Ochoa, S. Eiden, J. J. Sáenz, and F. Scheffold, Transport of light in amorphous photonic materials, Applied Physics Letters 91 (2007).

[6] L. Schertel, L. Siedentop, J.-M. Meijer, P. Keim, C. M. Aegerter, G. J. Aubry, and G. Maret, The structural colors of photonic glasses, Advanced Optical Materials 7, 1900442 (2019).

[7] D. Contini, F. Martelli, and G. Zaccanti, Photon migration through a turbid slab described by a model based on diffusion approximation. i. theory, Applied optics 36, 4587 (1997).

[8] V. A. Markel, Introduction to the maxwell garnett approximation: tutorial, J. Opt. Soc. Am. A 33, 1244 (2016).

[9] F. Scheffold, W. Härtl, G. Maret, and E. Matijević, Observation of long-range correlations in temporal intensity fluctuations of light, Phys. Rev. B 56, 10942 (1997).

[10] A. García-Martín, F. Scheffold, M. Nieto-Vesperinas, and J. J. Sáenz, Finite-size effects on intensity correlations in random media, Phys. Rev. Lett. 88, 143901 (2002).

[11] R. Samanta, R. Pierrat, R. Carminati, and S. Mujumdar, Speckle Decorrelation in Fundamental and Second-Harmonic Light Scattered from Nonlinear Disorder, Physical Review Applied 18, 054047 (2022).