# Spatio-temporal visualization of light transport in complex photonic structures

Lorenzo Pattelli, Romolo Savo, Matteo Burresi, Diederik S. Wiersma

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

## To cite this version:

Lorenzo Pattelli, Romolo Savo, Matteo Burresi, Diederik S. Wiersma. Spatio-temporal visualization of light transport in complex photonic structures. Light: Science and Applications, 2016, 5, pp.e16090. 10.1038/lsa.2016.90 . hal-01333961

HAL Id: hal-01333961

https://hal.sorbonne-universite.fr/hal-01333961

Submitted on 20 Jun 2016

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

ORIGINAL ARTICLE

Spatio-temporal visualization of light transport in
complex photonic structures

Lorenzo Pattelli<fnref n="1" />, Romolo Savo<fnref n="1" />,*, Matteo Burresi²,† and Diederik S Wiersma<fnref n="1" />,³

Spatio-temporal imaging of light propagation is very important in photonics because it provides the most direct tool available to study the interaction between light and its host environment. Sub-ps time resolution is needed to investigate the fine and complex structural features that characterize disordered and heterogeneous structures, which are responsible for a rich array of transport physics that have not yet been fully explored. A newly developed wide-field imaging system enables us to present a spatio-temporal study on light transport in various disordered media, revealing properties that could not be properly assessed using standard techniques. By extending our investigation to an almost transparent membrane, a configuration that has been difficult to characterize until now, we unveil the peculiar physics exhibited by such thin scattering systems with transport features that go beyond mainstream diffusion modeling, despite the occurrence of multiple scattering.

*Light: Science & Applications* (2016) **5**, e16090; doi:10.1038/lsa.2016.90; published online 6 May 2016

**Keywords:** diffusive approximation; disordered media; light transport; optical gating; wide-field imaging

INTRODUCTION

Light represents a useful and versatile probe to study our world. By tracking its evolution in the domains of time, space, frequency and reciprocal space, we are able to extract a wealth of information regarding the surrounding environment. Further insight is often gained when we are technically able to track this evolution in more than one such domain simultaneously, as demonstrated by a number of recent results from diverse research fields<fnref n="1" />⁻⁸. Of the different approaches, spatio-temporally resolved techniques undoubtedly offer the most straightforward approach, as underlined by the continuously growing interest in fast photography⁹⁻¹¹ and imaging¹²⁻¹⁵ applications. Further development of similar techniques is key to our ability to investigate increasingly more complex media, such as disordered or heterogeneous media, which are ubiquitous in most fields of study, from atmospheric physics to biology and cultural heritage preservation. In this context, one main task of optical investigations is to characterize the light transport process determined by multiple scattering. This characterization is typically achieved through repetitive measurements to retrieve meaningful microscopic transport parameters, such as the transport mean free path *l*ᵢ and the absorption length *l*ₐ, which are directly related to the properties of the investigated medium (that is, the internal structure, composition and/or density), and are key parameters to consider from both the fundamental and the application points of view.

Here, we experimentally combine wide-field spatial imaging and
sub-picosecond (ps) time resolution and show how the acquired

multi-dimensional dataset allows us to characterize previously
inaccessible scattering systems and to improve the accuracy of
measurements in standard situations. The ultrafast time resolution
needed for such applications represents a major experimental chal-
lenge. Established fast imaging techniques involve electronic-based
devices such as streak cameras or sequential scanning over a given
sample in a pump-probe configuration. Despite providing very
accurate results, many of these techniques are limited by low-time
resolutions¹³ or a long acquisition time¹⁴. Conversely, wide-field
images at sub-ps time scales are easily obtained with optical gating
techniques; however, their limited accuracy has limited their use to
qualitative pattern/shadowgram recognition rather than to absolute
intensity measurements across the entire field of view. In other words,
optical gating has been previously exploited as a convenient method to
see through a turbid medium obliterating the complex light transport
arising from its structure, not to study it directly.

In this study, we investigate light transport focusing on three classes
of disordered media that are extremely relevant for practical applica-
tions and yet still very difficult to characterize accurately, namely
samples with limited thickness, large-scale heterogeneities and semi-
transparent membranes. In each of these cases, transport properties
are retrieved by exploiting, for the first time, an all-optical gating
setup, which we developed to offer both wide-field acquisition and
high quantitative accuracy on a sub-ps time scale. Our spatio-temporal
investigation allows us to grasp the evolution of light transport in its
entirety, consequently revealing effects that were either overlooked or

<footnote n="1">¹European Laboratory for Non-linear Spectroscopy (LENS), Università di Firenze, Sesto Fiorentino (FI) 50019, Italy; ²Istituto Nazionale di Ottica (CNR-INO), Firenze (FI) 50125, Italy and ³Dipartimento di Fisica e Astronomia, Università di Firenze, Sesto Fiorentino (FI) 50019, Italy</footnote>

<footnote>*Current affiliation: Laboratoire Kastler Brossel, UMR 8552, CNRS, Ecole Normale Supérieure, Université Pierre et Marie Curie, Collège de France, 24 rue Lhomond, Paris 75005, France.</footnote>

<footnote>†Current affiliation: Gestione SILO Srl, via di Castelpulci 14/d, Scandicci (FI) 50010, Italy.</footnote>

<footnote>Correspondence: L Pattelli, Email: pattelli@lens.unifi.it</footnote>

<footnote>Received 23 September 2015; revised 13 January 2016; accepted 2 February 2016; accepted article preview online 5 February 2016</footnote>
misinterpreted to date, including a peculiar regime occurring in ultra-
thin disordered media that becomes apparent only when the spatial
and temporal domains are studied simultaneously and would not be
detectable with other state-of-the-art techniques. By comparing the
outcome of our analysis with that of the mainstream evaluation model
based on the diffusive approximation (DA), we describe how single-
domain investigations performed to date are susceptible to particularly
deceptive artifacts that can be clearly appreciated in a multi-domain
dataset. At the same time, as widely suggested in the literature, we
discuss how much more robust interpretations, even within the
diffusive framework, can be given in terms of transverse transport,
to which we gain direct access through our technique.

MATERIALS AND METHODS

Experimental setup

The optical scheme employed to resolve the time and space propaga-
tion of light in complex media is based on the cross-correlation gating
technique<sup>16,17</sup>, shown schematically in Figure 1a. Two synchronous,
collinear probe and gate pulses at different wavelengths impinge on a
2-mm-thick β-barium borate (BBO) crystal. When the two beams
overlap both spatially and temporally, a sum-frequency signal is
generated with intensity proportional to the cross-correlation of the
original pulses as a function of the delay-line position. This technique
has been used to fully characterize ultrafast pulse propagation<sup>18-20</sup>,
generally disregarding the spatial distribution of the converted signal
(Figure 1b). Nevertheless, cross-correlation gating techniques offer
several advantages, such as being unaffected by optical chirp (in
contrast to interferometric methods) and being compatible with
collinear schemes because possible spurious second harmonic con-
tributions from gate/probe beams, if any, can still be spectrally
separated. Indeed, a collinear geometry is particularly convenient to

mitigate possible distortion of upconverted images passing through the
BBO, which is a crucial aspect for accurate spatially resolved
experiments<sup>21</sup>. Finally, the gate and probe beams are perfectly
exchangeable, therefore allowing us to easily investigate any sample
at both wavelengths.

Upconversion imaging

To study the dynamics of light emerging from a complex specimen,
we image its exit surface on the BBO crystal. There, the original image
at the probe wavelength is upconverted to a different (signal)
wavelength by sum-frequency generation with a gate pulse. Non-
linear upconversion is a polarization- and wavevector-sensitive pro-
cess, and both of these aspects must be taken into account when
designing an imaging apparatus. In the current study, which addresses
multiple scattering systems, the polarization sensitivity is not relevant
because all polarization channels have the same optical properties.
However, the setup inherently allows to investigate separately different
polarization channels. In the phase-matching condition, the non-
linear crystal converts only those wavevectors laying within a relatively
narrow paraxial range in reciprocal space, which defines the angular
acceptance of the BBO crystal (see Supplementary Information). This
acceptance angle is weighted differently in the x- and y-directions of
the BBO surface (that is, along the image plane, see Figure 1c and 1d).
A double telecentric optical apparatus in a 4f configuration ensures
that the paraxial component emerging from the sample is preserved at
the crystal interface. This allows us to obtain a quantitatively accurate
reproduction of the original intensity distribution at each time step for
each point of the object field. This accuracy condition holds as long as
the transmitted wavevector distribution, however complex, does not
depend on the point and time of exit, as is the case with multiple
scattered light in disordered samples. In this case, the intensity carried

<figure><img src="image_4.png" /><figcaption>**Figure 1** Experimental setup and spatio-temporal resolution. (a) Two linearly polarized fs-pulses are provided by a Ti:Sa laser (probe in red) and a synchronously pumped optical parametric oscillator (gate in gray). The probe beam interacts with a slab sample, which is imaged on a BBO crystal with a 4f lens configuration. A small aperture (FF) allows the *k* vectors to be isotropically filtered, while a dichroic mirror (M2) guarantees signal and gate spatial superposition. A long-pass (F1) and a bandpass filter (F2) block the visible and unconverted light, respectively. (b) A typical cross-correlation time-resolved measurement acquired with a PMT. A fit with the cross-convolution of two sech² pulses returns a FWHM of 174 fs. (c) The spatial information carried by the pulse is retrieved using a CCD camera, whose resolution is influenced by the sum-frequency process. (d) A geometric sketch qualitatively showing how the phase-matching condition results in different angular acceptances in the extraordinary and ordinary directions. A slightly misaligned $\bar{k}_e$ and a largely misaligned $\bar{k}_o$ may result in the same $\Delta\theta$, therefore explaining the sharper resolution in the y-direction.</figcaption></figure>
by any (fixed) subset of wavevectors can be considered to be strictly proportional to the total transmitted intensity.

As seen in Figure 1c, the reduced set of available wavevectors degrades the sharpness of the retrieved image differently along the two axes. To obtain the same resolving power in both the x- and y-directions, the aperture of a Fourier filter (FF) can be reduced to act as a tunable isotropic wavevector filter. Several alternative upconverting configurations have been reported in the literature regarding full-frame imaging applications<sup>6,22–29</sup>, most of which employ a non-linear crystal in the Fourier plane of the imaging optics rather than in the image plane. Both configurations present different advantages, and the preference for their adoption depends on the specific application. In the common Fourier plane configuration, angular acceptance limits upconversion uniformity over the entire field of view; in our case, the image is formed with a cut-off set of *k* vectors (thus with lower resolution) while guaranteeing spatially uniform upconversion efficiency. Therefore, we adopted the latter configuration because our diffuse signal typically does not exhibit any sharp features, and a spatially uniform response represents a critical condition for attaining our quantitative accuracy goal over multiple orders of magnitude. The generated signal is collected by a lens inside a shielded box where a flip mirror allows us to switch between different detection schemes.

We initially performed standard time-resolved measurements using a photomultiplier tube (PMT), that is, collecting intensity without any spatial resolution. Figure 1b shows a typical cross-correlation measurement of our probe pulse passing through a test target together with a fit obtained from the convolution of two squared-hyperbolic secant functions (see Supplementary Information). The measurement does not exhibit any relevant satellites over several orders of magnitude. By switching the photon counter with a CCD camera, a transient image is detected (Figure 1c) that shows the spatial resolution of our setup (~22.6 lp mm⁻¹ in x and 11.3 lp mm⁻¹ in y at 0.5 contrast), which is comparable to previous publications<sup>6,24,26,29</sup>. The final resolution is obtained by reducing the aperture of the FF to <1 mm, until the same (isotropic) resolution of ~11.3 lp mm⁻¹ is obtained in both directions.

With regards to the overall sensitivity of the upconversion process, previous measurements performed with the same setup<sup>20</sup> give an estimate of the signal attenuation at which the noise level is eventually reached. Starting from a typical pulse energy of ~1 nJ at a probe wavelength of 810 nm, an upconverted signal can still be detected after a damping of 8 decades, taking into account the fraction of the probe beam that is usually lost by diffuse reflection and the limited solid angle subtended by the collection optics. In addition, several other parameters can be adjusted to further enhance sensitivity, such as increasing the integration time of the detector, using a non-linear crystal with a higher upconversion efficiency and a larger angular acceptance (for example, bismuth triborate) or increasing the fluence on the crystal. To this purpose, it is particularly convenient to increase the gate beam intensity as much as possible to avoid any alteration of the sample.

## Data evaluation models

By accessing both the spatial and temporal evolution of light emerging from a complex environment, an incredibly rich picture emerges that needs to be analyzed with the proper evaluation tools. Exact modeling of energy transport in disordered samples is provided by the radiative transfer equation, which models light propagation as a random walk of energy packets undergoing scattering and absorption events. The characteristic length, $l_s$, over which a packet is scattered on average is called the scattering mean free path, whereas the anisotropic scattering

factor, *g*, expresses the average degree of directionality ranging from completely random (*g*=0) to perfectly forward (*g*=1). Both these parameters are directly linked to the microscopic structure of the sample and to the properties of its basic scattering elements. However, a similarity relation exists<sup>30,31</sup> equating a random walk with a certain $l_s$ and *g* to a completely isotropic random walk with a rescaled mean free path of $l_t=l_s/(1-g)$. Thanks to this relationship, light transport is often studied in the isotropic diffusion approximation (DA) regime, where a single diffusion coefficient, $D=l_tv/3$, is used as the only transport parameter, with a consequent loss of insight concerning the microscopic properties of the sample. Moreover, the convenient description offered by the DA is known to be valid only when the transport mean free path is much smaller than the sample size. Conversely, its applicability to optically thin media is often debated and a standard approach to measure transport parameters in this regime is still missing. Several attempts have been made to derive extensions of the diffusive model for optically thin media<sup>32–34</sup>, and alternative techniques to measure transport parameters have been proposed that avoid using diffusion theory in the first place<sup>35–43</sup> and are based on single-domain (for example, temporal, spatial or angular) investigations.

In this study, we compare analyses of our experimental data with the mainstream diffusion theory and a Monte Carlo model of light transport. While the former method retains its appeal due to its simplicity, the latter represents a direct implementation of the radiative transfer theory, providing a computationally expensive but asymptotically exact solution to the transport problem. This allows us to unveil some of the particularly deceptive ways that the DA can fail when dealing with certain types of media.

## RESULTS AND DISCUSSION

In the following section, we present a few relevant experimental cases where exploiting our newly developed ultrafast imaging (UFI) setup, we were able to reveal structural and optical features of the investigated samples that could not have been probed unambiguously with current state-of-the-art time-resolved techniques and data evaluation models. These claims are further validated and quantified by a direct comparison of the experimental data with Monte Carlo simulations. The main advantages offered by our experimental technique, as well as other possible applications, are discussed in the final subsection.

### Limited thickness artifacts

We first tested our setup with the simplest possible disordered slab, a homogeneous isotropic sample made of TiO₂ nanoparticles embedded in a polymer matrix (see Supplementary Information). The sample thickness was 203 µm, with an average refractive index at 810 nm of *n*=1.52.

Figure 2a shows a normalized time-resolved transmission profile, in addition to a fit using the DA relation, which shows excellent agreement with the experimental data. The fit was performed with the diffusion coefficient, *D*, and absorption coefficient, $\mu_a = l_a^{-1}$, as the only free parameters, assuming $D=l_tv/3$ and $v=c/n$. The fit gives $l_t=24.1 \pm 0.6$ µm and $l_a=12.0 \pm 1.1$ mm (albedo 0.998). The ratio between the thickness of the slab and its transport mean free path (referred to as the optical thickness) is commonly used to express the amount of multiple scattering in a sample. According to the fit, we obtained an optical thickness of ~8.4; a value of 8 is usually considered to be the rule-of-thumb lower limit above which diffusive transport is regarded as an acceptable approximation<sup>44</sup>. Nonetheless, the retrieved coefficients predicted a total absorption of $A \approx 8\%$, which
<figure><img src="image_5.png" /><figcaption>**Figure 2** Homogeneous sample data analysis. (a) Experimental data (black squares) acquired with the photon counter (PMT) compared to an exponential fit (solid green line), and a diffusion theory fit with *D* and *μ*<sub>a</sub> as free parameters (solid red line). The inset shows the asymptotic decay time expected according to the DA for a sample with *L*=203 µm, *n*=1.52 and *μ*<sub>a</sub>=0 µm<sup>-1</sup> (solid black line). The experimentally retrieved value (dashed green line) is incompatible with the DA unless some absorption is introduced. (b) A Monte Carlo (MC) solution (solid red line) of the radiative transfer equation showing perfect agreement with the experimental data for a non-absorbing sample. The integrated intensity from UFI frames and the PMT data (cyan circles) are shown. (c) A set of frames acquired at different delays with our UFI setup. Each frame is averaged over different disorder realizations and is displayed normalized to its own maximum intensity. (d) Comparison of the experimental and simulated mean square width.</figcaption></figure>

is much higher than expected because both the host polymer and the TiO₂ nanoparticles are known to have vanishing absorption at near infrared wavelengths. Even considering a non-absorbing sample with an equal thickness and refractive index, the lowest possible decay time that diffusion theory could predict (τ<sub>min</sub><sup>DA</sup> = 6.23 ps, see Figure 2a, inset) would still be appreciably longer than that determined experimentally (τ = 6.01 ± 0.08 ps). The DA accommodates for such discrepancies by introducing an absorption artifact. This can intui- tively be seen in its relationship for the asymptotic decay time, τ⁻¹ = π²D/L<sub>eff</sub>² + μ<sub>a</sub>ν (where L<sub>eff</sub> is the effective thickness depending on the refractive index boundary conditions), the value of which can be made arbitrarily small by assuming the presence of absorption.

In this case, the analysis of the time evolution of the mean square width of the transmitted profile is particularly interesting in that it is inherently free from this absorption-to-scattering crosstalk effect<sup>13,45,46</sup>. The mean square width is defined as the variance, $w^2(t) = \int r^2 I(\vec{r}, t) d^2 \vec{r} / \int I(\vec{r}, t) d^2 \vec{r}$, where $\vec{r}$ is a two-dimensional vector in the transmission plane, and can be determined if the full spatio-temporal evolution of the profile $I(\vec{r}, t)$ is known. Due to its normalization, any amplitude factor (such as absorption) applied to the profile will cancel out exactly at any time, leaving the mean square width unaffected. We therefore further investigated the sample using our imaging configuration (which could in principle also substitute the photon counter following frame integration, see Figure 2b). The mean square width, as extracted by the obtained frames, clearly shows a linear increase (Figure 2c and 2d) that can easily be interpreted within the DA, which predicts $w^2(t) = 4Dt$. The experimental slope, as obtained from a linear fit, corresponds to

$D_{\text{exp}}^{\text{fit}} = 1746 \pm 21 \, \mu\text{m}^2 \, \text{ps}^{-1} \rightarrow l_t = 26.6 \pm 0.3 \, \mu\text{m}$, which is appreciably larger than the value retrieved from the time-resolved curve. We therefore ran Monte Carlo simulations (see Supplementary Information) to resolve this discrepancy, and eventually succeeded in simultaneously reproducing both the experimental time-resolved and mean square width evolution datasets using a *single* transport parameter, $l_t$ (Figure 2b and 2d). The Monte Carlo inversion procedure gives a value of $l_t = 25.5 \pm 0.5 \, \mu\text{m}$ (again corresponding to an optical thickness of ~8), which is perfectly compatible with both datasets without introducing absorption.

It should be noted that the transport mean free path, as retrieved from the mean square width slope in the diffusive framework, results in just a 4% overestimation of the actual best value retrieved by the Monte Carlo simulations. Such a level of accuracy is acceptable for a wide number of applications, where UFI techniques could be exploited to extend the applicability of the DA in this intermediate thickness range. Moreover, as widely suggested in the literature<sup>45,46</sup>, estimating transport properties from the mean square width rather than from time-resolved data is much more accurate and straightforward because in the former case, *a priori* knowledge of the absorption coefficient, refractive index contrast and sample thickness is not required. Indeed, the simple 4Dt dependence does not involve these parameters or the extrapolated boundary conditions. This simple determination of the slope is also insensitive to the exact determination of the time axis origin and the exact time-window considered for the fitting, both of which represent long-standing issues in the evaluation of time-resolved data<sup>47,48</sup>. In addition, for the same reason that the mean square width is independent of the absorption, it is also effectively
<figure><img src="image_6.png" /></figure>

independent of the integration time, as well as of any laser power
fluctuations. Interestingly, the time-resolved DA can deceptively
compensate for its limited accuracy by introducing an absorption
artifact, but lacks a solution when the correct value of absorption is
considered. This point is particularly relevant for investigations of
optically thin specimens, such as biological tissues⁴⁹.

As previously mentioned, our UFI technique is capable of replacing
completely integrated measurements with comparable sensitivity
(Figure 2b). Exploiting the full spatio-temporal information, the
resulting set of frames can be decomposed into a mean square width
and a time-resolved curve, providing two almost independent datasets
in a single measurement session. Indeed, the former does not depend
on the integrated intensity of the individual frames, whereas the latter
is primarily linked to the transport along the z-axis. In this respect, our
perfectly matching Monte Carlo validation, which is conducted
independently on two experimental datasets that depend differently
on different sets of parameters, is a clear indication of the quantitative
accuracy and reliability of our setup over many orders of magnitude in
both the temporal and spatial domains.

Large-scale heterogeneity artifacts

We also used our UFI technique to investigate more diverse and
complex media. A new polymer slab (L=190 µm) with a roughly
doubled scatterer density was manufactured to obtain an optically
thicker sample. Again, a perfect fit could be obtained in the diffusive
framework (Figure 3a) if the presence of absorption was included at an
unlikely integrated value for the total absorption, A≈7%
(l<sub>t</sub>=12.0 ± 0.4 µm and l<sub>a</sub>=13.2 ± 1.2 mm). In sharp contrast to the
previous experiment, a non-absorbing Monte Carlo fit was unable to
reproduce the time-resolved data, even when the observed experi-
mental decay time was perfectly matched (Figure 3b, dotted black
line). Furthermore, the mean square width exhibited a steeper than
expected slope (D<sub>exp</sub> fit = 1050 ± 13µm<sup>2</sup> ps<sup>-1</sup> → l<sub>t</sub> = 16 ± 0.2µm) as if
the diffusion process was enhanced along the in-plane directions, as
shown in Figure 3c.

Driven by this discrepancy, we further investigated the sample with
a scanning electron microscope (SEM), revealing a heavily layered
modulation of the scatterer density compared to the more homo-
geneously dispersed one in our first sample (Figure 4a and 4b). Taking
advantage of the versatility of Monte Carlo simulations, we modeled
our sample after our SEM images as being composed of five layers
with three different thicknesses and densities arranged symmetrically
with respect to the central layer. A new Monte Carlo brute force fit
assuming a fixed geometry and three free parameters, l₁, l₂ and l₃,
was eventually able to perfectly reproduce both the time-resolved and
mean square width data simultaneously, as shown in Figure 3b and 3c.
We obtained transport mean free paths of l₁ = 3.5 ± 0.5 µm,
l₂ = 21.5 ± 0.5 µm and l₃ = 11.0 ± 0.5 µm, respectively, for the high,
low and intermediate-density layers. The absorption coefficients were
all set to µₐ = 0 µm⁻¹. In this case, the mean square width/
time-resolved decomposition clearly allowed an unexpected degree
of complexity to be revealed and characterized thanks to the rigid
constraints of the double boundary datasets. Resolving this composite
structure showed that a layered heterogeneity can, in principle, mimic
the effects of anisotropic transport and therefore revealed a previously
unaddressed physical effect. This is especially relevant given the
pervasiveness of layered media (for example, coatings, atmospheric
physics and biological tissues), which often exhibit counter-intuitive
features that standard time-resolved techniques are still unable to
explain, as in the case of light transport across the human forehead⁵¹.

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

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

<figure><img src="image_9.png" /><figcaption>**Figure 3** Heterogeneous sample data analysis. (a) Photon counter (PMT) data (black squares) and best diffusive (DA) fit (solid red line) giving *l*<sub>1</sub>=12 µm and *l*<sub>a</sub>=13.2 mm. (b, c) Comparison of the PMT and UFI data to a multi-layered simulation (solid red line) and two homogeneous simulations. The black-dotted line refers to the simulation that reproduced the experimental decay time (*l*<sub>1</sub>=14.2 µm), while the black-dashed line shows the simulation that reproduced the mean square width slope (*l*<sub>1</sub>=15.7 µm). The decay-time-reproducing simulation fails to reproduce the mean square width slope, while the simulation that shows perfect agreement with the mean square width slope poorly matches the time-resolved decay, therefore revealing the inconsistency of the homogeneity hypothesis. The Monte Carlo (MC) curves were rescaled to match the experimental peak count values.</figcaption></figure>

Ultra-thin samples

Despite the advantages offered by decomposition into a mean square
width and a time-resolved curve, it is clear that the raw, non-
integrated data provide an irreducible set of information. A dramatic
illustration of this point is obtained by studying a thin biological
sample, such as a small strip cut from the dried skin of a slice of
grapefruit. As the SEM image reveals (Figure 4c), its structure consists
<figure><img src="image_10.png" /></figure>

**Figure 4** SEM side-view images of the analyzed samples. (a) Heterogeneously disordered sample. The central inset shows a transverse section, which exhibits a layered symmetric structure that we modeled with a high-density central layer ($L_1=26$ µm), a low-density interstitial layer ($L_2=60$ µm) and an intermediate-density outer layer ($L_3=22$ µm). The heterogeneity represented by the core layer is also appreciable at the optical microscope (left inset). The right insets highlight the density difference between regions 1 and 2. (b) A transverse section of the homogeneously disordered reference sample. Despite showing signatures of scatterer clustering, small-scale clusters/heterogeneities are homogeneously dispersed across the sample thickness and do not exhibit any layered arrangement. As confirmed by the excellent agreement obtained with a single parameter Monte Carlo simulation and other work on similar samples<sup>50</sup>, light perceives an effective, homogeneous scatterer density in contrast to the previous tightly layered sample. (c) Lateral section of the dried skin of the grapefruit slice.

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

**Figure 5** Grapefruit skin data analysis. (a) The time shift induced by the sample in the probe path (black circles vs. black squares). Knowing its thickness from previous scanning electron microscopy images, the effective refractive index is found as $n=1+c\Delta t/L\approx1.31$. The inset shows the actual appearance of the investigated specimens. (b) The frame set, each of which is averaged over a small portion of the sample. The transmitted intensity distribution shows a quasi-ballistic elliptic wavefront traveling at slightly different speeds (see Supplementary Information). (c) Horizontal cross-cuts of the wavefront and a Monte Carlo fit with *g* and *l*<sub>s</sub> as free parameters. The fit returns values of *g*=0.7 and *l*<sub>s</sub>=150 µm. (d) Simulated radial intensity profiles in this extremely optically thin regime show a significant breakdown of the similarity relation.

of a conglomerate of small flakes forming a corrugated slab ~85 µm thick. The dried sample appears to be a brittle, almost transparent membrane (Figure 5a, inset).

In this range of optical (and absolute) thicknesses, standard experimental techniques fall short because of the extreme time scales

involved, causing common diffusive modeling to fail dramatically. Furthermore, most of the signal will be ballistically transmitted, carrying almost no information concerning the optical properties of the sample. As Figure 5a shows, the light that is scattered inside the grapefruit membrane is <1/50th of the original intensity. Our
ultrafast imaging technique allows us to selectively address the light
that was retained for a longer time inside the sample whereas at the
same time spatially inspecting it as it propagates along the main
(transverse) slab extension. Figure 5b shows a collection of frames
acquired over a time-window of ~5 ps after pulse injection, corre-
sponding to a total path length greater than 13 times the sample
thickness. We see that the light spreads through the membrane with a
well-defined wavefront traveling inside its main plane, resulting in a
dramatic departure from any prediction compatible with the diffusive
framework. However, standard time-resolved or steady state investiga-
tions would still just measure a single decay time and a bell-shaped
profile, respectively, which could deceptively support an inappropriate
interpretation in terms of the DA.

As we will show, our ultrafast time- and space-resolved experiment
provides a deep insight into the light transport properties of complex
systems, which requires the development of a novel analysis metho-
dology. We demonstrate a proof-of-concept procedure to assess in-
plane transport properties in this extremely optically thin regime. For
the sake of simplicity, we will illustrate our investigation on a cross-cut
of the profile assuming isotropic transport. The general treatment is
analogous and will be the subject of future work. Figure 5c shows a
comparison between experimental cross-cuts at different times and the
corresponding time evolution of an isotropic Monte Carlo simulation
that succeeds in reproducing the wavefront speed, intensity decay and
overall shape at all given times. To achieve this collective agreement,
the anisotropic scattering coefficient, g, needs to be considered as an
additional parameter to express the degree of directional correlation
between consecutive scattering events. As mentioned, diffusion theory
is degenerate in g as long as the transport mean free path $l_t$ is kept
constant, as expressed by the similarity relation, $l_t = l_s/(1-g)$. This is
why retrieving any information concerning the angular scattering
function is particularly challenging, and one must therefore resort to a
regime where the DA ceases to hold<sup>43</sup>. Indeed, as shown by Figure 5d,
the observed transverse intensity patterns exhibit a significant break-
down in this degeneracy condition in thin slabs, with g playing a
substantial role in determining the shape and evolution of the traveling
wavefront, even when $l_t$ is constant. In particular, while the outer
wavefront propagates almost ballistically (see Supplementary
Information), its instantaneous position and peak intensity vary
appreciably with different combinations of g and $l_s$, therefore allowing
both to be retrieved with an estimated precision of ~10%. The set of
simulated curves shown in Figure 5c is obtained for $g=0.7$ and
$l_s = 150$ µm (corresponding to an in-plane $l_t = 500$ µm). However,
quantitatively accurate figures for this particular specimen would
require an analysis involving fully anisotropic simulations.

The main purpose of this experiment is to showcase the vast
possibilities opened by rich time- and space-resolved outputs. Our g
estimation results from a collective fitting routine involving multiple
spatio-temporal datasets and stands in contrast to common and less
robust g determination techniques involving a single scalar, such as the
attenuation coefficient of a collimated beam<sup>52,53</sup>. Accurate determina-
tion of single scattering anisotropy is crucial to directly probe the
shape, orientation and optical properties of single microscopic
scatterers and might therefore be of interest to a broad range of
fields, from realistic computer rendering of participating media<sup>54,55</sup> to
photo-therapeutic diagnostics of biological tissues<sup>56,57</sup>. Moreover,
compared to other available techniques, the determination of g and
$l_s$ is achieved in a single measurement session, which also simplifies the
retrieval procedure. The peculiar doughnut-shaped intensity pattern
that enables such a rich characterization of the scattering properties is
of course not unique to this particular sample, and we observed it in a

variety of different biological samples of both plant and animal origin,
usually differing only by the actual contour shape of the propagating
wavefront. This supports the general validity and usefulness of our
proof-of-concept evaluation technique in studying extremely complex
media, included biological tissues.

Discussion

As illustrated, our optically gated imaging technique offers several
advantages, including higher quantitative fidelity over a large field of
view with respect to Fourier plane configurations and sub-ps time
resolution compared to other electronically gated imaging solutions.
Gaining direct access to such time scales not only enables the study of
optical properties in thin or inherently minute specimens (for
example, ocular fundus, vascular walls, skin dermis and dental
enamel), but is fundamental to investigate light transport at the
mesoscopic level. From a technical viewpoint, the main strength of our
setup is its simplicity and lack of complex electronics. The inherently
synchronous nature of the pulses renders the entire setup insensitive to
timing jitter. Power fluctuation/drifts also do not pose a problem
when measuring the mean square width evolution, which also allows
us to dynamically adjust the integration time to enhance the signal-to-
noise ratio. This flexibility should make it possible to replicate the
setup with cheaper and more easy to use hands-off fiber lasers,
therefore widening the effective applicability of the technique.

Accessing both transverse and axial transport simultaneously
through a full spatio-temporal investigation allows us to correctly
evaluate the applicability of the common diffusion framework, which
can be subject to appreciable absorption overestimation in samples
presenting a limited optical thickness or large-scale heterogeneities. An
accurate retrieval of the absorption coefficient is particularly relevant
in many biomedical applications where its value is directly linked to
quantities such as chromophores concentration or blood oxygenation
levels. The observed discrepancy, which would be not apparent with
state-of-the-art spatial or temporal measurements, is revealed by
comparing standard time-resolved transmitted intensity decay to the
evolution of the mean square width, both of which can be acquired
simultaneously. Nevertheless, the most relevant advantage offered by
direct access to the transient mean square width expansion is that it
enables a robust and simple interpretation. In sharp contrast to other
common observables, such as the total transmittance and its decay
time, the mean square width is independent of absorption and largely
unaffected by refractive index contrasts and sample thickness, there-
fore overcoming long-standing problems posed by their precise
assessment. One prominent case where this feature could be decisive
is that of the experimentally observed decrease in the diffusive
coefficient of a turbid slab with decreasing thickness<sup>58</sup>, which was
indirectly determined through the transmittance decay time. A spatio-
temporal measurement in terms of the mean square width evolution
would allow us to directly probe the diffusion coefficient irrespective of
the exact boundary conditions, the incorrect estimation of which has
been cited as a possible cause for this apparently anomalous behavior<sup>32</sup>.

A different scenario emerges when studying samples with extremely
low optical thickness, which revealed a peculiar transport regime
where $l_s$ and g, rather than $l_t$, are the meaningful parameters
determining how radiative energy spreads through the sample. This
finding illustrates a further practical situation where light transport
must be studied as a spatio-temporal phenomenon in its entirety. This
point has not yet been properly considered when studying transport in
thin samples. Moreover, we have demonstrated how this transport
regime enables accurate retrieval of otherwise barely accessible
transport parameters, especially for this class of samples where well-established experimental techniques are still missing.

In general, several applications will be facilitated by the three-dimensional characterization capabilities of the setup (see also Supplementary Information), opening new possibilities to properly investigate structurally anisotropic media that are of paramount relevance in many industrial and biomedical applications, and whose correct theoretical modeling in terms of light transport is still an object of lively debate in the literature⁵⁹,⁶⁰. Other relevant open questions that might be tackled directly with a spatio-temporal investigation include experimental verification of the polarization dependence predicted for the diffusion coefficient⁶¹, as well as the occurrence of non-monotonic mean square width evolution in samples exhibiting Anderson localization⁴⁵ or the lateral expansion of Lévy type transport, where the experimental limit posed by the finite thickness of the sample might be circumvented⁶². Similarly, an UFI measurement of intensity profiles reflected by macro-porous bulk materials could access light spreading within the first few hundred femtoseconds from which the time-dependent diffusion constant could be determined⁶³ leading to information concerning the structure factor of the material, which would not be obtainable at long times.

## CONCLUSIONS

The field of light transport in complex or structured media still presents a multitude of experimental challenges that could significantly benefit from the advent of novel multi-domain investigation techniques. In this study, we demonstrated an experimental technique to go beyond the ps resolution barrier for spatial investigations of transient intensity profiles emerging from arbitrary samples. Gaining access to this time scale with a broad field of view is key to the characterization of a large class of samples, ranging from thin biological membranes to complex optical microdevices. As we demonstrated, spatio-temporal investigations enable the study of heterogeneous or semi-transparent media that are difficult to characterize with present state-of-the-art techniques. Further, when dealing with more standard samples, combined access to both transverse and axial transport helps reveal the presence of flaws or artifacts in the common diffusion-based data evaluation framework. Remarkably, in many practical situations, we found that low optical thickness can be exploited as an advantage rather than a drawback, therefore offering a complementary framework to study the broad class of samples that do not meet the validity conditions required by the DA. Nonetheless, the setup that we propose is generally applicable and non-invasive, usable at different wavelengths, with a broad field of view and high-time resolution. We therefore envision that it could be applied to a wide range of photonic applications for both the sub-ps physics it allows to investigate and its convenient wide-field acquisition, which does not require scanning over the region of interest.

## ACKNOWLEDGEMENTS

Authors wish to thank G Mazzamuto, F Utel and J Bertolotti for technical assistance and T Svensson, E Alerstam, S Gigan and P Wasylczyk for fruitful discussions. This work is financially supported by the European Network of Excellence Nanophotonics for Energy Efficiency and the ERC through the Advanced Grant PhotBots (Proj. Ref. 291349).

## AUTHOR CONTRIBUTIONS

MB, LP and RS conceived and carried out the time-resolved experiments. LP and RS manufactured the samples and LP ran the simulations. All authors discussed the results and contributed to the writing of the paper.

1. Drexler W, Morgner U, Ghanta RK, Kärtner FX, Schuman JS *et al.* Ultra high-resolution ophthalmic optical coherence tomography. *Nat Med* 2001; **7**: 502–507.

2. Obrig H, Villringer A. Beyond the visible-imaging the human brain with light. *J Cereb Blood Flow Metab* 2003; **23**: 1–18.

3. Dunsby C, French PMW. Techniques for depth-resolved imaging through turbid media including coherence-gated imaging. *J Phys D: Appl Phys* 2003; **36**: R207.

4. Engelen RJP, Sugimoto Y, Gersen H, Ikeda N, Asakawa K *et al.* Ultrafast evolution of photonic eigenstates in k-space. *Nat Phys* 2007; **3**: 401–405.

5. Minardi S, Eilenberger F, Kartashov YV, Szameit A, Röpke U *et al.* Three-dimensional light bullets in arrays of waveguides. *Phys Rev Lett* 2010; **105**: 263901.

6. Bassi A, Brida D, D'Andrea C, Valentini G, Cubeddu R *et al.* Time-gated optical projection tomography. *Opt Lett* 2010; **35**: 2732–2734.

7. Velten A, Willwacher T, Gupta O, Veeraraghavan A, Bawendi MG *et al.* Recovering three-dimensional shape around a corner using ultrafast time-of-flight imaging. *Nat Commun* 2012; **3**: 745.

8. Eggebrecht AT, Ferradal SL, Robichaux-Viehoever A, Hassanpour MS, Dehghani H *et al.* Mapping distributed brain function and networks with diffuse optical tomography. *Nat Photon* 2014; **8**: 448–454.

9. Gao L, Liang JY, Li CY, Wang LV. Single-shot compressed ultrafast photography at one hundred billion frames per second. *Nature* 2014; **516**: 74–77.

10. Nakagawa K, Iwasaki A, Oishi Y, Horisaki R, Tsukamoto A *et al.* Sequentially timed all-optical mapping photography (STAMP). *Nat Photon* 2014; **8**: 695–700.

11. Li ZY, Zgadzaj R, Wang XM, Chang YY, Downer MC. Single-shot tomographic movies of evolving light-velocity objects. *Nat Commun* 2014; **5**: 3085.

12. Gundlach L, Piotrowiak P. Femtosecond Kerr-gated wide-field fluorescence microscopy. *Opt Lett* 2008; **33**: 992–994.

13. Sperling T, Bührer W, Aegerter CM, Maret G. Direct determination of the transition to localization of light in three dimensions. *Nat Photon* 2012; **7**: 48–52.

14. Farid N, Harilal S, Ding H, Hassanein A. Dynamics of ultrafast laser plasma expansion in the presence of an ambient. *Appl Phys Lett* 2013; **103**: 191112.

15. Grumstrup EM, Gabriel MM, Cating EEM, Van Goethem EM, Papanikolas JM. Pump-probe microscopy: Visualization and spectroscopy of ultrafast dynamics at the nanoscale. *Chem Phys* 2015; **458**: 30–40.

16. Shah J. Ultrafast luminescence spectroscopy using sum frequency generation. *IEEE J Quant Electron* 1988; **24**: 276–288.

17. Yodh AG, Kaplan PD, Pine DJ. Pulsed diffusing-wave spectroscopy: high resolution through nonlinear optical gating. *Phys Rev B* 1990; **42**: 4744–4747.

18. Trebino R. *Frequency-Resolved Optical Gating: The Measurement of Ultrashort Laser Pulses*. Berlin: Springer, 2000.

19. Savo R, Burresi M, Svensson T, Vynck K, Wiersma DS. Walk dimension for light in complex disordered media. *Phys Rev A* 2014; **90**: 023839.

20. Burresi M, Cortese L, Pattelli L, Kolle M, Vukusic P *et al.* Bright-white beetle scales optimise multiple scattering of light. *Sci Rep* 2014; **4**: 6075.

21. Potenza MAC, Minardi S, Trull J, Blasi G, Salerno D *et al.* Three dimensional imaging of short pulses. *Opt Commun* 2004; **229**: 381–390.

22. Firester AH. Image upconversion: part III. *J Appl Phys* 1970; **41**: 703–709.

23. Chiou W. Geometric optics theory of parametric image upconversion. *J Appl Phys* 1971; **42**: 1985–1993.

24. Faris GW, Banks M. Upconverting time gate for imaging through highly scattering media. *Opt Lett* 1994; **19**: 1813–1815.

25. Devaux F, Lantz E. Ultrahigh-speed imaging by parametric image amplification. *Opt Commun* 1995; **118**: 25–27.

26. Abraham E, Bordenave E, Tsurumachi N, Jonusauskas G, Oberlé *J et al.* Real-time two-dimensional imaging in scattering media by use of a femtosecond Cr⁴⁺: forsterite laser. *Opt Lett* 2000; **25**: 929–931.

27. Lantz E, Devaux F. Parametric amplification of images: from time gate to noiseless amplification. *IEEE J Select Top Quant Electron* 2008; **14**: 635–647.

28. Pedersen C, Karamehmedovi E, Dam JS, Tidemand-Lichtenberg P. Enhanced 2D-image upconversion using solid-state lasers. *Opt Expr* 2009; **17**: 20885–20890.

29. Dam JS, Pedersen C, Tidemand-Lichtenberg P. High-resolution two-dimensional image upconversion of incoherent light. *Opt Lett* 2010; **35**: 3796–3798.

30. Graaff R, Aarnoudse JG, de Mul FFM, Jentink HW. Similarity relations for anisotropic scattering in absorbing media. *Opt Eng* 1993; **32**: 244–252.

31. Van de Hulst HC. *Multiple Light Scattering: Tables, Formulas, and Applications*. Amsterdam: Elsevier; 2012.

32. Elaloufi R, Carminati R, Greffet JJ. Time-dependent transport through scattering media: from radiative transfer to diffusion. *J Opt A: Pure Appl Opt* 2002; **4**: S103.

33. Xu M, Cai W, Lax M, Alfano RR. Photon migration in turbid media using a cumulant approximation to radiative transfer. *Phys Rev E* 2002; **65**: 066609.

34. Garofalakis A, Zacharakis G, Filippidis G, Sanidas E, Tsiftsis D *et al.* Characterization of the reduced scattering coefficient for optically thin samples: theory and experiments. *J Opt A: Pure Appl Opt* 2004; **6**: 725–735.

35. Hammer M, Roggan A, Schweitzer D, Muller G. Optical properties of ocular fundus tissues-an *in vitro* study using the double-integrating-sphere technique and inverse Monte Carlo simulation. *Phys Med Biol* 1995; **40**: 963–978.

36. Kienle A, Patterson MS. Determination of the optical properties of turbid media from a single Monte Carlo simulation. *Phys Med Biol* 1996; **41**: 2221–2227.

37. Pifferi A, Taroni P, Valentini G, Andersson-Engels S. Real-time method for fitting time-resolved reflectance and transmittance measurements with a Monte Carlo model. *Appl Opt* 1998; **37**: 2774–2780.

38. Alerstam E, Andersson-Engels S, Svensson T. White Monte Carlo for time-resolved photon migration. *J Biom Opt* 2008; **13**: 041304.
<figure><img src="image_12.png" /></figure>

39 Leonetti M, López C. Measurement of transport mean-free path of light in thin systems. *Opt Lett* 2011; **36**: 2824-2826.

40 Nilsson AMK, Berg R, Andersson-Engels S. Measurements of the optical properties of tissue in conjunction with photodynamic therapy. *Appl Opt* 1995; **34**: 4609-4619.

41 Dam JS, Yavari N, Sørensen S, Andersson-Engels S. Real-time absorption and scattering characterization of slab-shaped turbid samples obtained by a combination of angular and spatially resolved measurements. *Appl Opt* 2005; **44**: 4281-4290.

42 Karlsson H, Fredriksson I, Larsson M, Strömberg T. Inverse Monte Carlo for estimation of scattering and absorption in liquid optical phantoms. *Opt Expr* 2012; **20**: 12233-12246.

43 Svensson T, Savo R, Alerstam E, Vynck K, Burresi M *et al*. Exploiting breakdown of the similarity relation for diffuse light transport: simultaneous retrieval of scattering anisotropy and diffusion constant. *Opt Lett* 2013; **38**: 437-439.

44 Elaloufi R, Carminati R, Greffet JJ. Diffusive-to-ballistic transition in dynamic light transmission through thin scattering slabs: a radiative transfer approach. *J Opt Soc Am A* 2004; **21**: 1430-1437.

45 Cherroret N, Skipetrov SE, Van Tiggelen BA. Transverse confinement of waves in three-dimensional random media. *Phys Rev E* 2010; **82**: 056603.

46 Hu HF, Strybulevych A, Page JH, Skipetrov SE, van Tiggelen BA. Localization of ultrasound in a three-dimensional elastic network. *Nat Phys* 2008; **4**: 945-948.

47 Alerstam E, Andersson-Engels S, Svensson T. Improved accuracy in time-resolved diffuse reflectance spectroscopy. *Opt Expr* 2008; **16**: 10440-10454.

48 Bouchard JP, Veilleux I, Jedidi R, Noiseux I, Fortin M *et al*. Reference optical phantoms for diffuse optical spectroscopy. Part 1-Error analysis of a time resolved transmittance characterization method. *Opt Expr* 2010; **18**: 11495-11507.

49 Tuchin V. *Tissue Optics: Light Scattering Methods and Instruments for Medical Diagnosis*. Bellingham: SPIE Press; 2007.

50 Vos WL, Tukker TW, Mosk AP, Lagendijk A, IJzerman WL. Broadband mean free path of diffuse light in polydisperse ensembles of scatterers for white light-emitting diode lighting. *Appl Opt* 2013; **52**: 2602-2609.

51 Comelli D, Bassi A, Pifferi A, Taroni P, Torricelli A *et al*. *In vivo* time-resolved reflectance spectroscopy of the human forehead. *Appl Opt* 2007; **46**: 1717-1725.

52 Marchesini R, Bertoni A, Andreola S, Melloni E, Sichirollo A. Extinction and absorption coefficients and scattering phase functions of human tissues *in vitro*. *Appl Opt* 1989; **28**: 2318-2324.

53 Pickering JW, Prahl SA, Van Wieringen N, Beek JF, Sterenborg HJCM *et al*. Double-integrating-sphere system for measuring the optical properties of tissue. *Appl Opt* 1993; **32**: 399-410.

54 Pauly M, Kollig T, Keller A. *Metropolis Light Transport for Participating Media*. Berlin: Springer, 2000.

55 Premoże S, Ashikhmin M, Tessendorf J, Ramamoorthi R, Nayar S. Practical rendering of multiple scattering effects in participating media. In: Jensen HW, Keller A editors. *Proceedings of the Fifteenth Eurographics Conference on Rendering Techniques*. Switzerland: The Eurographics Association, 2004, pp 363-374.

56 Nilsson AMK, Sturesson C, Liu DL, Andersson-Engels S. Changes in spectral shape of tissue optical properties in conjunction with laser-induced thermotherapy. *Appl Opt* 1998; **37**: 1256-1267.

57 Mourant JR, Freyer JP, Hielscher AH, Eick AA, Shen D *et al*. Mechanisms of light scattering from biological cells relevant to noninvasive optical-tissue diagnostics. *Appl Opt* 1998; **37**: 3586-3593.

58 Kop RHJ, de Vries P, Sprik R, Lagendijk A. Observation of anomalous transport of strongly multiple scattered light in thin disordered slabs. *Phys Rev Lett* 1997; **79**: 4369-4372.

59 Kienle A. Anisotropic light diffusion: an oxymoron? *Phys Rev Lett* 2007; **98**: 218104.

60 Alerstam E. Anisotropic diffusive transport: Connecting microscopic scattering and macroscopic transport properties. *Phys Rev E* 2014; **89**: 063202.

61 Vynck K, Pierrat R, Carminati R. Polarization and spatial coherence of electromagnetic waves in uncorrelated disordered media. *Phys Rev A* 2014; **89**: 013842.

62 Bertolotti J, Vynck K, Pattelli L, Barthelemy P, Lepri S *et al*. Engineering disorder in superdiffusive Lévy glasses. *Adv Funct Mater* 2010; **20**: 965-968.

63 Svensson T, Vynck K, Adolfsson E, Farina A, Pifferi A *et al*. Light diffusion in quenched disorder: Role of step correlations. *Phys Rev E* 2014; **89**: 022141.

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

This work is licensed under a Creative Commons Attribution 4.0
International License. The images or other third party material in this

article are included in the article's Creative Commons license, unless indicated otherwise
in the credit line; if the material is not included under the Creative Commons license, users
will need to obtain permission from the license holder to reproduce the material. To view a
copy of this license, visit http://creativecommons.org/licenses/by/4.0/

Supplementary Information for this article can be found on the *Light: Science & Applications*' website (http://www.nature.com/lsa).