Deducing effective light transport parameters in optically thin systems

We present an extensive Monte Carlo study on light transport in optically thin slabs, addressing both axial and transverse propagation. We completely characterize the so-called ballistic-to-diffusive transition, notably in terms of the spatial variance of the transmitted/reflected profile. We test the validity of the prediction cast by diffusion theory, that the spatial variance should grow independently of absorption and, to a first approximation, of the sample thickness and refractive index contrast. Based on a large set of simulated data, we build a freely available look-up table routine allowing reliable and precise determination of the microscopic transport parameters starting from robust observables which are independent of absolute intensity measurements. We also present the Monte Carlo software package that was developed for the purpose of this study.


Introduction
Light represents a powerful and versatile tool for investigating the optical properties of complex media. The radiative transport features of a disordered sample are indeed deeply related to its microscopic structural [1][2][3] and chemical properties [4,5], which makes light propagation studies a desirable non-invasive and nondestructive diagnostic tool in the most diverse applications. Light transport in disordered media can be described in terms of the microscopic optical parameters l s , l a and g which represent, respectively, the scattering and absorption lengths and the single scattering anisotropy. In a random walk framework for light transport, the two lengths l s and l a (whose reciprocals are the scattering and absorption coefficients m s and m a ) represent the mean distance traveled between two consecutive scattering and absorption events, respectively. With regards to g, it is defined as the average cosine of the scattering angle θ, and it represents the first moment of the angular distribution of the scattered light intensity, also known as the phase function. All these parameters provide different insights on the sample composition, such as the density and scattering strength of its constituents, the concentration of pigments or absorbing chromophores, and the shape and size of the single scatterers. Accurate retrieval procedures for these parameters are therefore of paramount importance in a number of applications. The scattering anisotropy g is particularly relevant in the biomedical field, where tissues usually feature strongly anisotropic scattering whose characterization yields insights on the presence and properties of scattering elements such as blood cells, mitochondria, inner membranes and other cellular organelles [6]. In practice, however, after multiple scattering events the information on scattering anisotropy is typically lost and a new isotropic scattering length called transport mean free path = l l g 1 In this regime, the Diffusion Approximation (DA) represents an incredibly robust theoretical framework. As such it provides a complete set of straightforward analytical expressions to describe light transport both in the spatial and temporal domain, thus allowing to solve the inverse transport problem and eventually retrieve l t and l a . Nonetheless it significantly fails in those circumstances where its fundamental assumption of almost isotropic radiance is not verified, a condition that is truly fulfilled only asymptotically in space and time. As a consequence, a vast effort is devoted to characterize its precision, accuracy and validity range, which are known to become defective in low albedo scattering systems and, most importantly, media whose extension is not large enough to allow the onset of a multiple scattering regime [8,9]. Among these, thin slabs and membranes represent a very common and relevant geometry, especially in the biological field [10], where many tissues such as the ocular fundus [11], vascular walls [12], cellular cultures [2], skin dermis [13], skull bones [14] or dental enamel [15]-to name a few-are inherently available only within narrow ranges of minute sizes and thicknesses. For these systems, ballistic and quasi-ballistic light can give a significant contribution that must be taken into account. On one hand, low-order scattering leaves transmitted light mostly unaltered, thus enabling notable applications such as optical tomography or imaging through turbid media [16][17][18]: in this class of applications, the diffuse component of transmitted light represents unwanted background that must be properly filtered. On the other hand, only the diffuse component bears information as for the microscopic optical properties of the sample, whose retrieval is the subject of this work. In both cases, the ballistic and diffuse components can be straightforwardly separated in the time domain. However, while ballistic light by its own definition always preserves the same properties, the interpretation of the diffuse component must be carefully considered when the optical thickness decreases, i.e. when the diffusion approximation is no longer expected to hold. In this scenario, experimental data evaluation must rely on more accurate methods such as refined approximations of the Radiative Transfer Equation [19][20][21][22] or Monte Carlo simulations [23,24], the latter providing an exact solution to the Radiative Transfer Equation for any sample geometry, given that large enough statistics are collected.
Nonetheless, due to its simplicity, the diffusion approximation still retains a large appeal, and great efforts are constantly made to extend its validity range introducing all sorts of minor modifications [25][26][27]. Even at its standard formulation, diffusion theory casts an incredibly simple prediction on transverse transport which could be profitably applied in a thin slab geometry, since boundary and confinement effects play a less relevant role along the slab's main extension. A light beam impinging on a disordered slab is usually transmitted with an enlarged spatial profile which, according to the diffusion approximation, is Gaussian with a standard deviation growing linearly in time as = w t Dt 4 2 ( ) with a slope determined by the diffusion coefficient D. This quantity is more generally referred to as the Mean Square Width (MSW), and can be defined for an arbitrary intensity distribution r I t , ( )through the expression ò ò r r r r r r r where ρ is the distance from the optical axis on the exit plane. Experimentally, the temporal evolution of the mean square width can be measured from a collection of discrete, spatio-temporally resolved profiles r I t , i ( ) [28,29]. Its linear increase can be considered as a signature for transverse diffusion, with the obtained slope representing a valuable observable since the MSW is remarkably unaffected by the presence of absorption (which cancels out exactly at any t i ). Even more surprisingly, the diffusion approximation predicts the mean square width slope to be independent from both the slab thickness and its refractive index contrast with the surrounding environment. These properties are unique to the MSW, whereas typically other observables depend critically on absorption, slab thickness and boundary conditions [30,31]. Nevertheless, transverse transport has been so far largely disregarded, supposedly because it was hardly accessible experimentally.
The aim of this study is to characterize diffusive light transport in the infinite slab geometry over a wide range of optical parameters, testing for the first time to our knowledge how the linear prediction = w t Dt 4 2 ( ) from the diffusion approximation performs as the optical thickness of the sample decreases. In this respect, our results fill the gap with the already extensive characterization available in the literature regarding the breakdown of the diffusion model along the propagation axis [8,9,29,[31][32][33][34][35][36][37][38], finally providing a more complete picture of diffusive transport in this intermediate scattering regime for which only early-time characterizations are generally available. Moreover we tackle the problem of transverse transport in a confined geometry, which is one of fundamental interest especially when considering that this configuration allows the onset of multiple scattering also in semi-transparent media which are not usually associated with this transport regime [39,40]. Furthermore we envision that also from an experimental point of view our study will allow to exploit the large potential of novel optical investigation techniques based on the mean square width. To this purpose, we present a freely available look-up table (LUT) interface as a straightforward tool to fully take advantage of our simulated dataset and retrieve effective light transport parameters for optically thin systems, starting from robust observables which are independent of absolute intensity measurements.

Results and discussion
In the first part of our work we perform a systematic Monte Carlo study over a range of different optical properties aimed at testing the validity of the diffusion approximation (DA). In particular we test the linear prediction = w t Dt , against w t 2 ( ) values obtained using (1) on a single photon basis. Assuming a slab thickness of = L 1 mm we simulated different optical thicknesses (OT) ranging from 1 to 10 by varying the transport mean free path l t between 0.1 and 1mm. Each OT simulation has been run for 11 different values of the anisotropy factor g, defined as the average cosine of the scattering angle θ, between 0 and 0.99 and 16 values of the refractive index contrast n ranging from 0.6 to 2.2, for a grand total of 2816 simulations of 10 9 photons each (see Methods). Photons are injected normally to the slab plane from a d d r t ( ) ( ) pencil beam source and collected at the opposite surface (i.e. we only consider transmitted intensity). For the sake of convenience, since Fresnel reflection coefficients depend solely on the relative refractive index contrast = n n n in out , we kept = n 1 in constant while varying n out in order to have a consistent time scale over the whole set of simulations. The real time scale of any single simulation can be recovered by simple multiplication by the actual value of n in . Figure 1(a) shows a subset of simulated mean square width data obtained for typical optical properties of relevance in the bio-optical field (n = 1.4,g = 0.9), exhibiting a perfectly linear increase also in the optically thinnest case. As previously specified, mean square width values are exactly independent of absorption, which thus has been excluded from the simulations. The retrieved values of D evaluated as 1/4 of the variance slope have been divided by the expected value = D l c 3 DA t and arranged in a hyper-surface of relative deviations ( figure 1(b)). The obtained volume is sampled in a discrete set of points in the (n, g, -OT 1 ) space where the single simulations have been performed. Noise coming either from limited statistics or fitting uncertainties is largely reduced by applying a local regression algorithm using weighted linear least squares and a 2nd degree polynomial model as provided by the Loess MATLAB model ( figure 1(c)). This allows us to eventually obtain an accurate, smoothly sampled volume suitable for finer interpolation (figure 1(d)).

Mean square width expansion
A few comments are in order. Firstly, the present investigation is intended to focus on long-time/asymptotic transport. To this purpose, the diffusion coefficient D has been evaluated by the linear slope of the mean square width (MSW) in a time window ranging from 4 to 9 decay times, as determined from time-resolved curves. Depending mainly on the optical thickness of the sample, there is an early-time range where the MSW exhibits a super-linear increase. We carefully checked that the aforementioned fitting time range was always largely excluding such non-linear time range in order to address safely the asymptotic slope, as confirmed for example in figure 1(a).
Secondly, it is well known that most biological soft tissues share a refractive index equal or close to = n 1.4 in [41]. This is supposedly the reason why refractive index variations have so far been disregarded in similar multiparameter investigations [42][43][44][45]. Nonetheless, we included the refractive index contrast as a simulation parameter because, especially in the case of thin slabs, the range of interest for n is undoubtedly wider, spanning from well below 1 to as high as 2. The first case for small n is of interest for cases where specimens are enclosed in glass slides, or laid or immersed in different substrates/solutions, whereas the high values for n have been included envisioning possible applications of our study to metal oxides and similar highly scattering materials, which are extremely relevant, for instance, for coatings and in photovoltaics. Looking at the obtained data, two features are immediately noticeable. Firstly, the diffusion approximation appears to always underestimate the actual spreading rate, of course recovering agreement for higher optical thicknesses as expected. A second, finer feature occurs in the close proximity of n = 1, particularly evident at low g and OT values. Both these features arise from the interplay between geometric and boundary conditions. In particular, we found that the presence of internal reflections in a thin layer geometry helps to selectively hold inside the slab those photons that happen to draw statistically longer steps, as we discuss in more detail in a related work [39]. For the purpose of this study, it is worth noting that the mean square width slope exhibits a distinct pattern of characteristic deviations from the diffusion approximation, which can therefore be exploited as a guide to unambiguously retrieve the intrinsic microscopic transport properties of a given sample.

Decay time
It is also interesting to consider how the decay time is affected by different configurations. In the diffusion approximation, for a non-absorbing medium we have: e is the effective thickness of the medium and z e represents the extrapolated boundary conditions for a given refractive index contrast [46,47].
Figures 2(a) and (b) respectively show a typical set of time-resolved transmittance decays and the hypersurface of relative deviations from the diffusion approximation prediction. Two main features are worth commenting when comparing to the previous results for the mean square width slope. First of all the obtained lifetime deviations are more significant, reaching down to 20 % of the expected value for the highest values of g and n. It is indeed known that > n 1 refractive index contrasts are harder to be taken into account even when applying appropriate boundary corrections and even at higher optical thicknesses. Secondly, deviations in both directions are possible, with the t t DA ratio taking values both above and below 1. This might help explaining some experimental evidences obtained in thin disordered samples that are still debated at a fundamental level [19,[48][49][50][51][52], as further discussed in a related paper [39]. Our findings stress the importance of an accurate and precise modeling of the index contrast, which we think has been often overlooked, for example when a symmetric averaged contrast is used to model asymmetric experimental configurations [19,50,52].
Despite the vast literature regarding the validity range of the diffusion approximation in the time domain [8,9,[31][32][33][34][35]38], a comprehensive understanding of the interplay between optical thickness, refractive index contrast and absorption is still a debated topic. It is a common assumption that the diffusion approximation fails gradually with decreasing optical thickness, with = OT 8 being customarily considered as the lower threshold under which the introduced error starts to be significant [9]. Nevertheless, it was recently shown how a non- absorbing, > OT 8 slab sample withñ 1.5 exhibits a transmittance decay time such that the diffusion approximation is unable to provide any real solution at all [29], thus suggesting that the breakdown of the diffusion approximation might step in abruptly depending on the interplay between different parameters other than the optical thickness. In the same work, by resorting to Monte Carlo simulations it was possible to correctly retrieve the transport mean free path and absorption coefficient beyond the evaluation capabilities of the diffusion approximation. Notably, the observed deviations from both the experimental decay time and the mean square width expansion are in perfect agreement with our simulations.

Absorption
By combining our mean square width and decay time analysis some unique insight on the effects of absorption is gained. As already mentioned, absorption has not been considered as a simulation parameter since its effects can be accounted for exactly both in w t 2 ( ) and τ. Concerning the mean square width, its value is completely unaffected at any time by the presence of any amount of absorption, while the asymptotic decay time is expected to shift exactly to τ being the decay time in the non-absorbing case. This is equivalent to saying that whenever we know the specific  The problem with absorption is that both scattering and absorption can deplete specific intensity from a given position, time and direction (an effect sometimes referred to as absorption-to-scattering cross-talk). Hence, retrieving its unknown value from experimental data has been to date a challenging task. Besides that, it is often reported that the diffusion approximation is expected to hold only for weakly absorbing media since the onset of the properly diffusive regime requires long trajectories to contribute dominantly to transport properties, whereas these are selectively suppressed by absorption [53]. This explains why absorption is often considered as a major hindrance in the correct assessment of transport properties [54][55][56], if not even an invalidating condition for certain optical parameter measurements [57][58][59]. For this reason, techniques capable of directly accessing the mean square width recently aroused a great deal of interest given the absorptionindependent nature of the variance expansion [28,60,61]. From our perspective, the full potential of MSW measuring techniques is still to be fully unraveled, and will eventually play a key role among the most accurate characterization techniques of both scattering and absorption properties, as we will demonstrate in the following section.

Look-up table approach
In recent years, the increasing availability of massively parallel computation capabilities is fostering Monte Carlo simulations as a viable method to solve the inverse transport problem, i.e.retrieving the microscopic properties that give rise to certain temporal or spatial profiles. Two different Monte Carlo-based approaches are represented by fitting and look-up table (LUT) routines respectively. In the former method, one is usually required to solve iteratively the forward problem until some convergence condition between experimental and simulated data is satisfied. Due to the impractical computational load, fitting routines to date have exploited rescaling properties of the radiative transfer equation to adapt a limited set of pre-simulated Monte Carlo data to experimental measurements [11,34,42,[62][63][64][65][66]. To limit the occurrence of 'scaling artifacts', fitting approaches need to perform rescaling on a single-photon basis and therefore need to store each exit time and position separately [64,65]. Bin-positioning strategies are also known to represent a possible source of artifacts, requiring complex correction strategies to be deployed [66]. Finally, while for the semi-infinite geometry a single dimensionless Monte Carlo simulation can be rescaled both in terms of absorption (exploiting (4)) and scattering mean free path, the computational cost increases in the case of finite thickness geometries, where different scattering mean free paths values must be simulated separately. This is why only few examples can be found in the literature that deal with this configuration [42].
Taking advantage of the large set of simulations performed, we developed a look-up table (LUT) routine to demonstrate the capabilities of a combined mean square width/decay time approach. In the literature, several LUT approaches have been proposed, based on both experimental [44] and simulated data [43,45,[67][68][69][70][71][72]. In contrast with a standard fit, look-up table routines rely on single scalar parameters directly linked to transport properties, such as, typically, the total amount of transmitted/ballistic/reflected light from a slab. This triplet of observables, often referred to as T tot , T coll and R tot , has been extensively exploited to retrieve optical parameters through Monte Carlo-LUT routines [43,45,[67][68][69][70][71][72], despite the fact that such absolute intensity measurements are extremely challenging [43,70] and prone to unpredictable systematic errors [72]. Moreover, while the natural propensity of Monte Carlo-based data evaluation techniques is towards the study of optically thin media, integrated transmittance/reflectance quantities do actually loose their usefulness as the sample thickness decreases [11], since the acquired signal will be eventually dominated by light that has been either specularly reflected or ballistically transmitted through the sample, thus carrying almost no information about its properties. Even more fundamentally, one of the main assumptions in the integrating sphere theory is that of a lambertian diffuse profile [67], which is clearly not holding for thinner systems.
In this respect, for the first time to our knowledge, we designed a look-up table relying on observable quantities that are both free from any absolute intensity assessment and well into the multiple scattering regime, such as the asymptotic decay time and the mean square width expansion slope. This offers several advantages over existing solutions: both the asymptotic decay time and the mean square width slope can be measured without any reference to the excitation intensity, therefore there is no need to calibrate the source or the detector; decay time determination is also not strictly connected to any particular detection geometry, which stands in contrast with other typical techniques often requiring a particular configuration of collection fibers, integrating spheres or angular measurements.
because of the asymptotic nature of both τ and w t 2 ( ), the actual temporal response function or spatial excitation spot size are eventually irrelevant to their accurate determination.
precise determination of the origin of the time axis (i.e. the exact time of pulse injection), while being dramatically relevant in many analogue situations, is here made completely irrelevant since both the decay time and the linear increase of the mean square width do not exhibit any critical dependence on the exact delay at which they are determined, provided that it is sufficiently large.
with respect to Monte Carlo-based fitting routines, a look-up table routine is more suitable for real time solving of the inverse problem since it does not involve any iterative procedure. While this guarantees better performances, on the downside we must note that it is less clear how to define the uncertainty of retrieved values. Typical approaches involve mapping the relative error of retrieved parameters over a broad range of independent simulations, in order to give a numerical estimate.
several issues typically associated to fitting routines are also obviated. Once that the two scalar parameters are calculated with the proper, original binning, they can be virtually rescaled arbitrarily without the risk of introducing any binning-related artifact.
the problem of correct bin positioning is also removed. Midpoint positioning adopted in our case represents the exact solution for the linear increase of the mean square width. While this is not the case for a monoexponential decay, it can be again shown trivially that midpoint positioning does leave the decay constant exactly unmodified.
as we will demonstrate in the following, a possible implementation of our LUT routine allows to retrieve m a and l t . It must be nonetheless stressed again that none of the simulations that make up our look-up table include absorption, which stands in contrast with the customary practice of most LUT approaches demonstrated to date. This allows to deal with a simulation phase space of reduced dimensionality, which represents a huge saving on the computational burden of Monte Carlo simulations.
-finally, since there is no need to directly simulate absorption nor add it after the simulation, it is also not necessary to store exit times and positions on a single photon basis. This allows to hugely reduce the output size for each simulation and streamline its handling, thus allowing for better statistics to be collected.
For the sake of simplicity, our present Monte Carlo-LUT demonstration is limited to the retrieval of pairs of transport parameters, e.g.l t and g (and therefore l s ) assuming that absorption is known, or l t and m a with known g, which is a realistic and common assumption in similar works, especially those involving biological samples )and n L L in 0 respectively. Let us consider first the general case of a medium with unknown scattering mean free path and absorption coefficient, but known g. Given the sample thickness and the internal/external refractive indices, it is possible to rescale the continuously interpolated version of the mean square width hyper-surface ( figure 1(d)) in order to slice it at the given refractive index contrast n. The obtained 2D surface will feature an iso-level curve corresponding to the measured mean square width slope, which will basically give the expected -OT 1 , i.e.l t , in a completely absorption-independent way. For the test sample described above, fitting the simulated mean square width values yields a slope of m -337750 m ps 2 1 . By intersecting this iso-level curve with the simulated value of g = 0.95, we eventually retrieve the best -OT 1 estimate ( figure 3(a)). Nonetheless, even if the scattering anisotropy is not known a priori, plugging reasonably bounded values into the routine helps getting an estimate of how an uncertainty on g spreads over l t and eventually m a . Once that also -OT 1 is determined, it is sufficient to read the expected absorption-free decay time value stored in t . The second implementation of our routine allows to retrieve l t and g assuming that m a is known. A common case is that of vanishing absorption, which is often encountered when studying for example metal oxide powders with NIR radiation. We therefore reconsider the same iso-level curve of figure 3(a) which we replot in figure 3(c).  [29] recently highlighted the robustness of the mean square width as a valuable experimental observable for the retrieval of the microscopic transport parameters. In that work, the authors consider a homogeneous isotropic sample made of TiO 2 nanoparticles (g = 0.6, with vanishing absorption at the working wavelength of 810 nm) embedded in a polymer matrix; sample thickness was measured to be 203 μm and the average refractive index at the working wavelength is 1.52, close to that of many biological tissues. The observed experimental decay time and mean square width slope were t = 6.01 ps and m -6984 m ps 2 1 and a value of = m l 25.5 m t was found with a bruteforce Monte Carlo inversion procedure against the experimental data, involving the simulation of many combination of optical parameters. Here, by feeding the same experimental parameters into our look-up table routine, we instantly find a value of = m l 25.7 m t , in good agreement with the value reported in the paper. Thorough evaluation of errors should be performed on a wide range of parameters, both from simulations and experimental data, which is beyond the scope of this work. Nonetheless it is clear that, especially at lower thicknesses where the diffusion approximation is more defective, our routine offers very accurate retrieving capabilities as compared to other slab geometry fitting and/or LUT approaches [42]. Uncertainties as low as a few percent with respect to simulated data have been demonstrated in other works for the semi-infinite geometry using integrated intensities as the input parameters. It might be questioned whether this kind of uncertainty is realistic, since integrated sphere measurements themselves suffer of both random and systematic errors of similar magnitude in the first place [43]. On the contrary, the slope of the mean square width and the transmittance decay time can be typically determined with better precision, accuracy and robustness, since their scalar value depends on the collective arrangement of multiple points of a curve.
As a last point, it is interesting to discuss on possible extensions of our routine applicability. At least a third input observable in addition to the decay time and the mean square width slope needs to be known in order to retrieve simultaneously all three transport parameters at once from an unknown medium. A possible candidate could be represented by the asymptotic tail slope of a steady state profile, which should exhibit an appreciable dependence on g at lower optical thicknesses. For all practical purposes, this asymptotic tail slope would feature all the previously listed advantages, with the possible exception of the last one, because of the need to add absorption ex post. Other relative parameters could be exploited, taking advantage of their g dependence, such as the transmittance rising time [38], even though such parameters would clearly not benefit from all the aforementioned points. To sum up, look-up table methods are very general in their nature and consequently can be profitably applied in a number of practical use cases. Of course, in order to tackle more complex samples (e.g. multilayered or anisotropic slabs) more observables are needed. Nonetheless, we believe that, whenever possible, mean square width and decay time measurements should always be preferred and included in every LUT-based retrieval routine because of their intrinsic robustness.

Software
For the purpose of this work we developed a new multilayered Monte Carlo C++ software library called MCPlusPlus, whose source code can be freely consulted, used, and contributed at http://www.lens.unifi.it/ quantum-nanophotonics/mcplusplus. This software was developed from scratch aiming at enriching existing multi-layer Monte Carlo software such as MCML [23] or CUDAMCML [24]. Being developed entirely in C++, the program takes extensive advantage of object-oriented programming paradigm (OOP), which is particularly suited to model a random walk problem [73]. Since pieces of code can be encapsulated in reusable objects, OOP offers several advantages including scalability, modularity, ease of maintenance and abstraction. Of equal importance is the fact that OOP naturally lends itself as a tool to describe a high-level interface to the software itself. Indeed, MCPlusPlus comes as a shared library rather than an executable package. A Python interface to the library is also provided so that simulations are extremely easy to set up and run through very simple scripts. Scriptability proves indeed to be very useful as it easily allows to loop over the parameter space when building a look-up table or performing a brute-force fit. We believe this to be a key strength of our package which improves considerably on existing multi-layer Monte Carlo software and will hopefully encourage its adoption.
Random walk implementations of light transport fall into the category of so-called 'embarrassingly parallel' problems, whose solution can largely benefit from the increasing availability of parallel computing architectures such as GPUs (graphics processing units) and multi-core CPUs. Yet, despite delivering fastest performance, GPUs still present some limitations. We therefore developed our software to be run on CPUs in order to better meet our goal. In particular, as our approach is based on a look-up table, our main targets are numerical accuracy, reliability and reproducibility rather than execution speed. Indeed, current GPU implementations of the light transport problem in scattering media are more often targeted at realistic biomedical applications involving complex meshes which would otherwise present an overwhelming computational burden, rather than fundamental and statistical studies such those underlying this study. However we must note that, despite running on CPUs, MCPlusPlus performance is not sacrificed as we can still exploit the ubiquitous multi-core architecture of modern computers; performance close to GPU is soon matched on a small computing cluster or even on a single multi-core workstation. CPU code also ensures maximum hardware compatibility, while GPUbased implementations are hardware or even vendor specific. This should make the adoption of MCPlusPlus as broad as possible so to also encourage further expansion of the dataset on which our look-up table is based. Additionally, developing software for a pure CPU architecture has generally less complications than writing GPU-compatible code; plenty of software libraries and high-quality pseudo-random number generators (PRNG) are widely available for the CPU, providing us with the flexibility and freedom that we needed for the purpose of this work.
Finally, the magnitude of simulations performed for this and a related study based on the same software [39] is particularly significant, which alone poses several challenges. In particular, performing simulations with a large number of walkers (10 10 ) requires the use of 64-bit PRNGs in place of the more common 32-bit implementations, which would introduce a statistically significant truncation in the step length distribution. Accordingly, the correct representation of 64-bit generated random variates requires the use of long double floating point notation. Both these requirements are straightforward on a CPU architecture, as opposed to GPUs, supporting our preference for the former.
Regarding the simulation output, MCPlusPlus allows to collect statistics on photon exit times, exit points and exit angles. The software provides a powerful histogramming interface that allows to specify any number of simple or bivariate histograms, so that both steady state or time-resolved statistics can be extracted very easily. Nevertheless, if desired it is also possible to obtain raw, non-histogrammed data on a photon-by-photon basis. This can be useful for some particular applications, for example it allows one to introduce photon absorption after a single simulation has been run. All output is provided in the widespread HDF5 binary file.

Simulations
The look-up table (LUT) dataset which we described is based directly on the hyper-surfaces shown in figures 1 and 2, which result from a fine sampling of the (n, g, respectively, while varying l s , g, and n out . For each configuration we simulated 10 9 photons emitted from a pencil beam source impinging normally on the slab. Photons are propagated inside the scattering material through a standard random walk algorithm, where scattering lengths follow an exponential distribution and scattering angles are generated using the well known Henyey-Greenstein function. To build the look-up table, we first fit a single exponential decay to the time-resolved transmitted intensity curve to obtain a first estimate of the decay time t ¢. We then repeat the same fit, this time on the range t ¢ 4 -t ¢ 9 , to obtain the final estimate for τ. This ensures that the fitting is done at times long enough to extract the asymptotic value of τ and adds consistency to the fitting method between different simulations. Once τ is obtained, we find D by performing a linear fit on the mean square width w t 2 ( ) as a function of time (see figure 1), limiting the fit to the same t 4 -t 9 range; the lower limit is chosen so as to always exclude the early-time photons before the onset of the diffusive regime, while the upper limit is to avoid the noise found at very long times due to insufficient statistics. It might be called into question whether it is appropriate to use the decay time as a time unit for the mean square width evolution, since the former is mainly determined by transport properties along the thickness direction, while the latter occurs along the plane. A range based on the decay time provides indeed a convenient way of defining a consistent, self-tuning fitting window across the whole dataset. This simple choice is also advocated under practical reasons, since the decay time is undoubtedly the actual temporal unit that eventually dictates-both in real and numerical experiments-the signal-to-noise ratio. In this respect, every diffusion coefficient within our simulated phase space has been determined under equal noise conditions. No less important, limiting our investigation to a long-time window is also relevant under a more technical point of view: i.e. it renders irrelevant for all practical purposes the specific choice of both the spatial source distribution and the phase function. For the same reason, polarization can be safely ignored and does not play any role in our analysis.
Values obtained for τ and D for each simulation are eventually arranged in the form of a hyper-surface as shown in figures 1(b) and 2(b). In order to neutralize the noise originating from statistic fluctuations and fit uncertainty, we consider each simulated n-slice separately and smooth the data through a Loess fitting routine (range parameter set to 0.25) as shown in figures 1(c) and 2(c). Smoothed slices are eventually put back together to perform a cubic interpolation along the refractive index contrast axis to obtain a hyper-surface for D and τ that can be evaluated for any triplet in the (n, g, -OT 1 ) parameter space (figures 1(d) and 2(d)). Interpolation has been performed separately on the  n 1 and  n 1 regions of the parameter space due to the sharp firstderivative discontinuity occurring in n = 1.

Conclusions
In this paper we studied the radiative transfer problem in the infinitely extended slab geometry and compared its exact numeric solution (obtained by means of Monte Carlo simulations) to the predictions obtained by diffusion theory. For the first time, to our knowledge, both transverse and axial transport are addressed using, respectively, the mean square width growth rate and the transmittance decay time as figures of merit. These two observables are to be measured at late times, thus providing valuable insight on light transport properties well into the diffusive regime. Our investigation provides a complete characterization of how the diffusive approximation gradually fails over a range of optical thicknesses from 10 to 1-a range where resorting to the diffusion approximation is a questionable yet common practice-considering both the refractive index contrast and the scattering anisotropy. The mean square width growth rate is of particular interest since its standard prediction according to the diffusion approximation depends solely on the diffusion coefficient, as opposed to other timeresolved observables where also the thickness and refractive index contrast usually play a critical role.
With regards to the mean square width expansion rate, we found that deviations from the simple diffusion approximation prediction exist especially at low optical thicknesses and scattering anisotropy, and that they are always in the form of an under-estimation of the actual rate. Notably, when considering the case of high scattering anisotropy which is most relevant in many applications, the magnitude of the observed deviation remains limited even at low optical thicknesses as compared to the lifetimes, which instead deviate significantly even including extrapolated boundary conditions. In cases where quantitative accuracy is critical or the sample is optically thin, non-approximated methods should be adopted. In this context, readily available look-up tables (LUT) based on Monte Carlo simulations offer a convenient strategy to obviate the computational cost of numerical simulations in retrieving the microscopic transport properties, while retaining comparable accuracy. Based on the large simulated dataset, we presented a novel look-up table based on the combination of the decay time and mean square width slope as input parameters, which offer a series of relevant advantages over existing LUT solutions. Prominently, they can be measured without the need of any absolute intensity measurement nor time axis origin determination. On the simulation side, the look-up table is easily built since both the lifetime and the mean square width slope are unaffected by binning strategies, plus there is no need to simulate absorption nor add it subsequently in contrast with previously available solutions. Finally we apply our LUT strategy against a real experimental test sample with controlled properties. The retrieved values are comparable to those obtained by a brute-force Monte Carlo fit.
The creation of a custom look-up table is streamlined by taking advantage of a fully scriptable Monte Carlo software, which we developed. A documented version of the software is freely available at http://www.lens. unifi.it/quantum-nanophotonics/mcplusplus, which we hope will encourage feedback, adoption and contribution. The presented look-up table can be queried online at http://www.lens.unifi.it/quantumnanophotonics/mcplusplus/lut as an additional evaluation tool for samples with an intermediate optical thickness.