1.

Introduction

The Monte Carlo (MC) method for radiative transport is known as the gold standard for the description of light scattering in turbid media, especially when their optical density does not satisfy the assumptions of the diffusive approximation. The presence of weakly scattering or even optically clear regions inside a scattering medium can alter significantly the overall propagation of light, which is relevant for a number of cases including for instance, in the biomedical field, the thin cerebrospinal fluid layer,¹ the vitreous humor,² and optically cleared tissues in general.³

In all these cases, the reduced scattering mean free path $l_{\mathrm{s}}^{\prime }$ of these materials is comparable to or even larger than the thickness of these layers, thus requiring an MC description of light propagation. Owing to the intrinsic rescalability of the MC trajectories and sample geometry,^4,5 this typically holds true irrespective of the actual (absolute) length scales involved. In practice, however, scattering samples that are both optically thin and physically thin pose additional challenges that have been rarely considered in the literature. A part of these difficulties arise from the fact that the energy velocity of light remains fixed upon geometrical rescaling operations, and therefore the typical time scales associated with the radiative transport process can become extremely fast for minute specimens, making time-domain experiments more difficult to perform as they require ultrafast laser sources and, more importantly, ultrafast detection apparatuses with sub-ps resolution.

Nonetheless, studying light transport at these extreme time scales can bring valuable insights. For instance, the details of early sub-ps transients are associated with the single scattering regime even in diffusive media,^6–9 reveal structural properties related to the heterogeneity of scattering samples^10,11 and carry information about the scattering process down to the micrometer scale, which is also relevant for the study of natural and biological materials.^12–14 Prominent examples include thin tissues or small bioptic samples,^15–18 cell suspensions¹⁹ or optical tissues which are often in the submillimeter range,^20–23 thus requiring femtosecond-scale pulsed excitation for their optical characterization.

As the micrometer-scale becomes relevant, analog considerations apply to the spatial domain too. While one can often consider a collimated laser beam as a reasonable representation of a pencil beam, this may not be valid anymore if the thickness of the scattering layer becomes comparable to the waist of the impinging beam. In a typical experiment on very small specimens, tightly focused beams are typically used, which however are associated with a strong divergence after a short Rayleigh length. Thus, the need for an accurate description of tightly focused beams is of particular relevance in these circumstances.

Finally, due to the fast transport dynamics in thin media, which is dominated by ballistic or quasiballistic propagation, large MC simulations are usually required to accumulate a sufficient statistics and reach the multiple scattering regime. Simulations comprising billions or even more trajectories are not infrequent in this case, which however require extra care to make sure that the generation of random variates follows the intended distributions and is not affected by finite precision truncation effects.

In this report, we discuss a few modifications that can be introduced into a standard MC implementation to mitigate the possible shortcomings related to ultrafast time-of-flight measurements. Simulation results are then compared with experimental data from an all-optical gating setup capable of recording the temporal evolution of the transmittance profile from a simple plane-parallel scattering slab at the transition between the ballistic and diffusive regime, showing good agreement with the numerical profiles.

2.

Materials and Methods

2.1.

Source Description

Light transport in optically thin samples is associated with several deviations from the simple diffusive behavior, which are particularly interesting to study in the time domain.^24–27 In this context, MC simulations are successfully used to reproduce and explain the observed transport properties, especially when the samples have a macroscopic extent.

For much smaller sample sizes, much shorter time scales become relevant which require the development of radically different experimental setups and detection schemes needed for sub-ps resolution,^8,28,29 which could explain why time-domain experiments in this regime are rarely reported.

Optical gating techniques offer a straightforward strategy to achieve the required resolution by all-optical means, typically exploiting the Kerr effect or a nonlinear frequency generation process. In either case, an ultrafast laser source is required for the fast gate pulse, which is typically obtained via broadband Ti:Sa laser sources. When operated in the femtosecond pulse regime, the onset of mode-locking via a combination of self-phase modulation and group velocity dispersion leads to the formation of stable pulses³⁰ with a characteristic temporal profile

Eq. (1)

$$I(t;\mathrm{\Delta }{\tau }_p)=\sqrt{\frac{{\gamma }_p}{2}}{\mathrm{sech}}^2({\gamma }_{p}t),$$

with ${\gamma }_p=2\ln (1+\sqrt{2})/\mathrm{\Delta }{\tau }_p$ and $\mathrm{\Delta }{\tau }_p$ being the full width at half-maximum of the pulse. This profile shape is characterized by exponentially decaying tails, which thus extend much further in time with respect to a Gaussian pulse profile. What can be typically measured is the autocorrelation of a ${\mathrm{sech}}^2$ pulse with itself, which still presents the same exponentially decaying tails and can still be expressed analytically as

Eq. (2)

$$I_{\mathrm{ac}}={\gamma }_{\mathrm{ac}}{\mathrm{csch}}^2({\gamma }_{\mathrm{ac}}t)[{\gamma }_{\mathrm{ac}}t{\coth }^2({\gamma }_{\mathrm{ac}}t)-1],$$

with ${\gamma }_{\mathrm{ac}}=2\ln (1+\sqrt{2})/\mathrm{\Delta }{\tau }_{\mathrm{ac}}$. Notably, the expected spectrum for a ${\mathrm{sech}}^2$ pulse, as obtained by a Fourier transform operation, is still given by a ${\mathrm{sech}}^2$ curve as

Eq. (3)

$$\tilde{I}(\omega ;\mathrm{\Delta }{\tau }_p)=\frac{{\pi }^2}{2{\gamma }_p}\left[{\mathrm{sech}}^2\right(\frac{{\pi }^2}{2{\gamma }_p}(\omega -{\omega }_0)\left)\right].$$

Figure 1 shows an example of temporal and spectral characterization of a $\approx 100\mathrm{fs}$ pulse from a Ti:Sa laser exhibiting the typical ${\mathrm{sech}}^2$ pulse and spectrum profiles. The output beam from this laser is then used to pump a passive parametric oscillator to generate a synchronous replica at a wavelength of 1510 nm, to be used interchangeably as the probe or gate pulse in the optical gating experiment.

Fig. 1

(a) Autocorrelation measurement of the Ti:Sa pulse. The envelope profile formed by the autocorrelations fringes is fitted with the autoconvolution of a ${\mathrm{sech}}^2(t)$ pulse, returning a full width at half maximum for the convoluted signal equal to 145 fs, corresponding to a pulse width of $\mathrm{\Delta }{\tau }_{\mathrm{p}}=94\mathrm{fs}$. (b) Spectrum of the Ti:Sa pulse centered at ${\lambda }_0=810\mathrm{nm}$, fitted with a ${\mathrm{sech}}^2(\omega -{\omega }_0)$ function.

The hyperbolic secant square distribution can be expressed as the logistic distribution

Eq. (4)

$$f(x;\mu ,s)=\frac{\exp (-\frac{x-\mu }{s})}{s(1+\exp (-\frac{x-\mu }{s}){)}^2}=\frac{1}{4s}{\mathrm{sech}}^2\left(\frac{x-\mu }{2s}\right),$$

which has a closed-form cumulative distribution function

Eq. (5)

$$F(x;\mu ,s)=\frac{1}{1+\exp (-\frac{x-\mu }{s})}=\frac{1}{2}+\frac{1}{2}\tanh \left(\frac{x-\mu }{2s}\right),$$

whose inverse or quantile function is straightforward to implement to generate random ${\mathrm{sech}}^2$ variates in an MC routine

Eq. (6)

$$Q(\xi ;\mu ,s)=\mu +s\ln \left(\frac{\xi }{1-\xi }\right),$$

with $\xi$ uniformly distributed in the open interval (0, 1).

Experiments with small, submillimeter specimens typically require to focus the incoming beam tightly so that its beam waist is much smaller than the extent of the sample. Tightly focused beams are characterized by a strong divergence and a short Rayleigh length, which can easily become smaller or comparable to the thickness of the scattering layers. For instance, a beam focused down to a $10\text{-}\mu \mathrm{m}$ spot size can have a Rayleigh length of $\approx 50\mu \mathrm{m}$.

Due to its relevance, the accurate modeling of Gaussian beams has been addressed by several authors. Traditionally, trajectories are initialized to reproduce a certain focal distribution rather than an actual focal point.³¹ This approach has the advantage of being very easy to implement, even though it does not give the correct distribution of trajectories after or before the focal plane, which can be a critical factor when trying to evaluate accurately the deposition of energy in a tissue.

More advanced approaches have been proposed in recent years^32–34 which can be however more or less suitable for certain applications, lead to infinite trajectory densities or be significantly more computationally intensive.

A particularly elegant approach that is able to correctly reproduce the profile of a Gaussian beam at every position was proposed by Milsom³⁵ for the modeling of reverse saturable absorption. The approach is based on a stochastic ray-optics description which is perfectly compatible with an MC framework and based on the actual physical parameters used to describe the Gaussian beam.

Briefly, following this approach each ray is defined in Cartesian coordinates as the trajectory connecting two points separated along $z$ by a distance much larger than the Rayleigh length $z_R=\pi w_0^2/\lambda$, where $w_0$ is the spot size defined as the $1/e^2$ radius. A straightforward choice for the two $z$ coordinates are at the focusing lens and at the beam focal plane. Assuming that the beam starts with a certain waist $w_s$ at a distance $z=-d$ and arrives at the focal plane $z=0$ with a beam waist of $w_0$, then the pairs of points defining the rays in the bundle have coordinates $(x_s,y_s,-d)$ and $(x_f,y_f,0)$ given as

Eq. (7)

$$\begin{array}{cc}{x}_s=\frac{w_s}{2}{\mathrm{erf}}^{-1}(2{\xi }_1-1),& y_s==\frac{w_s}{2}{\mathrm{erf}}^{-1}(2{\xi }_2-1)\\ x_f=\frac{w_f}{2}{\mathrm{erf}}^{-1}(2{\xi }_3-1),& y_f=\frac{w_f}{2}{\mathrm{erf}}^{-1}(2{\xi }_4-1)\end{array},$$

where ${\xi }_{1\dots 4}$ is four independent random numbers uniformly distributed in the open interval (0, 1). An illustrative representation of a resulting ray bundle is shown in Fig. 2 for a beam focusing at the entrance of a sample ($z=0$) with a beam waist of $8\mu \mathrm{m}$.

Fig. 2

Views of a typical GRB focused in (0, 0, 0), with a beam waist of $8\mu \mathrm{m}$ in the focal plane. The envelope and intensity profiles formed by the rays follow the standard profile of a Gaussian beam at all $z$ planes.

Additionally, due to the ray-optics nature of the model, the so obtained Gaussian beam can be easily transformed using standard ray transfer matrix analysis when passing through interfaces, plates, lenses, and their combinations. The resulting bundle of rays will still follow the expected beam profiles as predicted by the theory.³⁵

3.

Results

The experimental measurements rely on an optical gating setup based on the collinear sum-frequency generation between a probe and a gate pulse impinging on a nonlinear crystal. This setup has been already used to characterize the time-resolved transmittance of different samples, including small specimens with a thickness of only a few micrometers, as described elsewhere,^9,43 and has been recently expanded to provide wide-field time-resolved imaging to study light propagation both in scattering media²⁹ and integrated photonic devices.⁴⁴

As a first measurement, the final thickness of the free-standing scattering slab can be estimated by performing a time-of-flight measurement through a transparent region of the sample, based on the time delay between the main transmitted peak and its first internally reflected replica. This measurement returned a thickness of $L=164\mu \mathrm{m}$, as shown in Fig. 3(a).

Fig. 3

Time-domain measurements. (a) Experimental time trace relative to a transparent region of the sample, together with an MC simulation for a slab with $n_{\mathrm{pol}}=1.509$ and thickness $L=164\mu \mathrm{m}$. The ${\mathrm{sech}}^2$ pulse shape shows the same qualitative features of the experimental curve, however the experimental data results from the convolution between two ${\mathrm{sech}}^2$ pulses with different widths. (b) Comparison between a transmittance measurement and an MC simulation with a transport mean free path of $41.2\mu \mathrm{m}$ and an absorption coefficient of ${\mu }_{\mathrm{a}}=4.5\times {10}^{-5}\mu {\mathrm{m}}^{-1}$, as determined by the mean square width growth data. The red and blue curves show simulation results obtained using a ${\mathrm{sech}}^2$ or a Gaussian pulse shape with equal full width at half-maximum. The difference between the two is only appreciable during the early transient (inset). (c) Mean square width data extracted from time-resolved images.

A standard time-domain transmittance measurement taken with a photo-multiplier tube is shown in Fig. 3(b), together with the output of two MC simulations, showing the difference between a ${\mathrm{sech}}^2$ and a Gaussian pulse shape. While changing the pulse shape has a negligible effect in the multiple scattering regime (as expected), their difference is clearly visible during the early transient, which bears useful information on single scattering parameters.⁹ Our experimental setup also allows us to study the evolution of the spatial profiles at different delays using a CCD camera (Andor iKon M912). A useful parameter that can be extracted from these spatiotemporal profiles, and which we use in conjunction with the transmittance decay rate, is their associated mean square width calculated as

Fig. 4

Time-resolved imaging measurements. (a) Set of acquired frames, each normalized to its maximum intensity. (b) Radial average of the spatial profiles for each frame, compared with MC data with $l_{\mathrm{s}}^{\prime }=41.2\mu \mathrm{m}$. Data points on the tails of early profiles are omitted for clarity. The amplitude of each curve is matched to that of the MC profile corresponding to the same delay.

The obtained spatiotemporal distributions are compared with those obtained by MC simulations in Fig. 4(b), showing good general agreement. The calibration of the spatial scale was characterized via transient imaging of a resolution target (not shown), returning a resolution of $\approx 13.5\mu \mathrm{m}/\mathrm{px}$. For the MC simulations, a ray-bundle beamwidth of $10\mu \mathrm{m}$ has been set based on a transient picture of the focal plane without the sample, which shows a single-pixel spot. The importance of using a Gaussian ray bundle (GRB) beam as opposed to a more common Gaussian pencil beam (GPB) spatial profile can be appreciated in Fig. 5. Therein, it can be seen that even for a sample with a thickness that is $4\times$ the transport mean free path, a residual difference between the two models is still appreciable in transmittance as a ringing effect due to the finite divergence of the GRB model. This difference is dominated by the early transient where the spatiotemporal profiles are most different. As can be seen, the GRB model shows a more gentle spatial profile as opposed to the sharp ballistic peak of the GPB case, even though the two models converge to very similar profile shapes after a short transient of $\approx 1\mathrm{ps}$.

Fig. 5

Numerical comparison between (a) spatial and (b) early spatiotemporal profiles obtained with a GPB or GRB source term focused on at the entrance of a $164\text{-}\mu \mathrm{m}$-thick sample with $l_{\mathrm{s}}^{\prime }=41.2\mu \mathrm{m}$. Especially at early times (spatiotemporal profiles cover a time span between 0 and 1.2 ps), the GRB source results in a more natural spatial spread of transmitted intensity due to its natural divergence.

4.

Discussion

In this paper, we presented an MC adaptation of experimental data measured in the spatial and temporal domains on thin scattering films with a transmittance decay rate of the order of few picoseconds. This requires the use of tightly focused, ultrashort (i.e., sub-ps) laser pulses, which we have also modeled in a simple MC implementation inspired by MCML. Due to the fast decay rate, and the need to gather enough statistics to represent spatial profiles at different delays, MC simulations must be run with a large number of trajectories. When generating random numbers required to gather such a large statistics, care should be taken to avoid hitting the limitations of simple PRNGs and the direct implementation of quantile functions involving logarithms, as these could be associated with undesired truncation effects due to the finite precision of floating-point representations.

We should note that, in recent years, several alternative strategies to the accurate but computationally expensive MC approach have also been proposed to improve the description of transient light transport in small specimens. For instance, some authors suggested introducing an effective frequency-dependent diffusion coefficient,^26,45 heuristic definitions for the boundary conditions⁴⁶ or crossover lengths between the ballistic and diffusive regimes.⁴⁷ Additionally, ongoing progress is being made with a number of analytical and semi-analytical methods to solve the radiative transfer equation (RTE) in the time domain, some of which can also be applied to model the transmittance through bounded geometries. Efficient $P_3$ solutions exist for the time-domain transmittance through a slab geometry,⁴⁸ while other solutions are being developed following the rotated reference frame method, which however are used to study either the spatial domain or semi-infinite geometries^49,50 and may still exhibit small deviations for samples with $L\approx l_{\mathrm{s}}^{\prime }$. Computationally efficient solutions have been also reported using a colocated radial basis function method,⁵¹ the pseudospectral method,⁵² or within the framework of the broad class discrete ordinate methods,⁵³ e.g., by combining it with spectral methods and Laplace transforms,⁵⁴ convection schemes,⁵⁵ finite-volume,⁵⁶ or Fourier continuation methods.⁵⁷ Finally, a promising approach has been recently proposed based on a dynamic radiative transfer system,⁵⁸ which should be in principle capable of modeling light transport down to the femtosecond scale. However, to the best of our knowledge, none of these methods has been tested against experimental transmittance data on such short time scales.

Indeed, experimental studies on transient transport optically thin layers are still sporadic in the literature, despite the existing theoretical and numerical evidence that light transport in this geometry exhibits intriguing and counterintuitive properties.^29,59 Given the increased availability of compact femtosecond-pulsed laser sources offering turn-key operation simplicity, optical gating detection schemes may become more widespread and bridge the gap toward the time and spatial scales that are relevant for submillimeter samples. Notably, all-optical gating techniques can typically provide a temporal resolution that is at least one order of magnitude higher than streak cameras, while providing wide-field (2D) imaging capabilities. Moreover, the probe-gate illumination scheme typically allows to enhance the upconverted signal by just increasing the intensity on the gate beam, and is therefore compatible with the study of delicate tissues and samples.

Robust and efficient methods to generate and handle large MC simulations will be increasingly needed to support the experimental investigations in this range of samples and geometries, as well as to validate the modeling approaches being developed for the description of early transients. At the sub-ps scale, the density of trajectories is potentially spanning over huge dynamic ranges,⁵⁸ requiring the accurate simulation of extremely rare events such as exceedingly long steps tens of times larger than the average mean free path which may eventually dominate the transport properties at very large delays. The contribution of such extremely rare events remains likely negligible for most practical purposes: however, it could be interesting to check whether these extreme outliers arising naturally in semi-transparent media are compatible and/or necessary for the verification of the invariance property for the average pathlength under Lambertian illumination, which should in principle remain valid to an arbitrary degree of precision.^60,61 Not least, on a more fundamental level, studying this extreme transport regime could provide insights regarding how far we can push the analogy between the scattering of light waves and the random walk model, on which all MC implementations implicitly rely.