Time-domain Raman analytical forward solvers

A set of time-domain analytical forward solvers for Raman signals detected from homogeneous diffusive media is presented. The time-domain solvers have been developed for two geometries: the parallelepiped and the finite cylinder. The potential presence of a background fluorescence emission, contaminating the Raman signal, has also been taken into account. All the solvers have been obtained as solutions of the time dependent diffusion equation. The validation of the solvers has been performed by means of comparisons with the results of “gold standard” Monte Carlo simulations. These forward solvers provide an accurate tool to explore the information content encoded in the time-resolved Raman measurements. c © 2016 Optical Society of America OCIS codes: (300.6450) Spectroscopy, Raman; (170.5660) Raman spectroscopy; (170.3660) Light propagation in tissues; (170.5280) Photon migration; (170.6280) Spectroscopy, fluorescence and luminescence. References and links 1. P. Matousek and M. D. Morris, Emerging Raman Applications and Techniques in Biomedical and Pharmaceutical Fields (Springer, 2010). 2. C. A. Froud, I. P. Hayward, and J. Laven, “Advances in the Raman Depth Profiling of Polymer Laminates,” Appl. Spectrosc. 57(12), 1468–1474 (2003). 3. N. J. Everall, “Confocal Raman microscopy: common errors and artefacts.,” Analyst 135(10), 2512–2522 (2010). 4. P. Matousek, and N. Stone, “Recent advances in the development of Raman spectroscopy for deep non-invasive medical diagnosis,” J. Biophotonics 1(5), 7–19 (2013). 5. K. Sowoidnich, J. H. Churchwell, K. Buckley, A. E. Goodship, A. W. Parker, and P. Matousek, “Photon migration of Raman signal in bone as measured with spatially offset Raman spectroscopy,” J. Raman Spectrosc. 47(2), 240–247 (2015). 6. P. Matousek, I. P. Clark, E. R. C. Draper, M. D. Morris, A. E. Goodship, N. Everall, M. Towrie, W. F. Finney, and A. W. Parker, “Subsurface Probing in Diffusely Scattering Media Using Spatially Offset Raman Spectroscopy,” Appl. Spectrosc. 59(4), 393–400 (2005). 7. N. Stone, R. Baker, K. Rogers, A. W. Parker, and P. Matousek, “Subsurface probing of calcifications with spatially offset Raman spectroscopy (SORS): future possibilities for the diagnosis of breast cancer,” Analyst 132(9), 899–905 (2007). 8. P. Matousek, E. R. C. Draper, A. E. Goodship, I. P. Clark, K. L. Ronayne, and A. W. Parker, “Noninvasive Raman Spectroscopy of Human Tissue In Vivo,” Appl. Spectrosc. 60(7), 758–763 (2006). 9. W. J. Olds, E. Jaatinen, P. Fredericks, B. Cletus, H. Panayiotou, and E. L. Izake, “Spatially offset Raman spectroscopy (SORS) for the analysis and detection of packaged pharmaceuticals and concealed drugs,” Forensic Sci. Int. 212(1-3), 69–77 (2011). 10. D. Yang and Y. Ying, “Applications of Raman Spectroscopy in Agricultural Products and Food Analysis: A Review,” Appl. Spectrosc. Rev. 46(7), 539–560 (2011). 11. C. Conti, C. Colombo, M. Realini, and P. Matousek, “Subsurface analysis of painted sculptures and plasters using micrometre-scale spatially offset Raman spectroscopy (micro-SORS),” J. Raman Spectrosc. 46(5), 476–482 (2015). 12. J. Johansson, S. Pettersson, and S. Folestad, “Characterization of different laser irradiation methods for quantitative Raman tablet assessment,” J. Pharm. Biomed. Anal. 39(3-4), 510–516 (2005). Vol. 24, No. 18 | 5 Sep 2016 | OPTICS EXPRESS 20382 #268152 Journal © 2016 http://dx.doi.org/10.1364/OE.24.020382 Received 10 Jun 2016; revised 13 Jul 2016; accepted 13 Jul 2016; published 25 Aug 2016 13. P. Matousek and A. W. Parker, “Bulk Raman Analysis of Pharmaceutical Tablets,” Appl. Spectrosc. 60(12), 1353– 1357 (2006). 14. M. V Schulmerich, J. H. Cole, K. A. Dooley, M. D. Morris, J. M. Kreider, S. A. Goldstein, S. Srinivasan, and B. W. Pogue, “Noninvasive Raman tomographic imaging of canine bone tissue," J. Biomed. Opt. 13(2), 020506 (2008). 15. J.-L. H. Demers, S. C. Davis, B. W. Pogue, and M. D. Morris, “Multichannel diffuse optical Raman tomography for bone characterization in vivo: a phantom study," Biomed. Opt. Express 3(9), 2299–2305 (2012). 16. J.-L. H. Demers, F. W. L. Esmonde-White, K. A. Esmonde-White, M. D. Morris, and B. W. Pogue, “Next-generation Raman tomography instrument for non-invasive in vivo bone imaging,” Biomed. Opt. Express 6(3), 793–806 (2015). 17. J. Steinbrink, H. Wabnitz, H. Obrig, A. Villringer, and H. Rinneberg, “Determining changes in NIR absorption using a layered model of the human head,” Phys. Med. Biol. 46(3), 879–896 (2001). 18. P. Matousek, M. Towrie, A. Stanley, and A. W. Parker, “Efficient Rejection of Fluorescence from Raman Spectra Using Picosecond Kerr Gating,” Appl. Spectrosc. 53(12), 1485–1489 (1999). 19. J. Wu, Y. Wang, L. Perelman, I. Itzkan, R. R. Dasari, and M. S. Feld, “Three-dimensional imaging of objects embedded in turbid media with fluorescence and Raman spectroscopy,” Appl. Opt. 34(18), 3425–3430 (1995). 20. P. Matousek, N. Everall, M. Towrie, and A. W. Parker, “Depth profiling in diffusely scattering media using Raman spectroscopy and picosecond Kerr gating,” Appl. Spectrosc. 59(2), 200–205 (2005). 21. F. Ariese, H. Meuzelaar, M. M. Kerssens, J. B. Buijs, and C. Gooijer, “Picosecond Raman spectroscopy with a fast intensified CCD camera for depth analysis of diffusely scattering media,” Analyst 134(6), 1192–1197 (2009). 22. S. Sundarajoo, E.L. Izake, W. Olds, B. Cletus, E. Jaatinen, and P. M. Fredericks, “Non-invasive depth profiling by space and time-resolved Raman spectroscopy, ” J. Raman Spectr. 44(7), 949–956 (2013). 23. E. L. Izake, B. Cletus, W. Olds, S. Sundarajoo, P. M. Fredericks, and E. Jaatinen, “Deep Raman spectroscopy for the non-invasive standoff detection of concealed chemical threat agents,” Talanta 94, 342–347 (2012). 24. A. Dalla Mora, E. Martinenghi, D. Contini, A. Tosi, G. Boso, T. Durduran, S. Arridge, F. Martelli, A. Farina, A. Torricelli, and A. Pifferi, “Fast silicon photomultiplier improves signal harvesting and reduces complexity in timedomain diffuse optics,” Opt. Express 23(11), 13937–13946 (2015). 25. A. Dalla Mora, D. Contini, S. Arridge, F. Martelli, A. Tosi, G. Boso, A. Farina, T. Durduran, E. Martinenghi, A. Torricelli, and A. Pifferi, “Towards next-generation time-domain diffuse optics for extreme depth penetration and sensitivity,” Biomed. Opt. Express 6(5), 1749–1760 (2015). 26. F. Martelli, T. Binzoni, A. Pifferi, L. Spinelli, A. Farina, and A. Torricelli, “There’s plenty of light at the bottom: statistics of photon penetration depth in random media,” Sci. Rep. 6, 27057 (2016). 27. D. Contini, F. Martelli, and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt. 36(19), 4588–4599 (1997). 28. F. Martelli, S. Del Bianco, A. Ismaelli, and G. Zaccanti, Light Propagation through Biological Tissue and Other Diffusive Media: Theory, Solutions and Software (SPIE, 2010). 29. N. Everall, T. Hahn, P. Matousek, A. W. Parker, and M. Towrie, “Photon Migration in Raman Spectroscopy,” Appl. Spect. 58(5), 591–597 (2004). 30. R. A. Desiderio, “Application of the Raman scattering coefficient of water to calculations in marine optics,” Appl Opt. 39(12), 1893–1894 (2000). 31. A. Bray, R. Chapman R, and T. Plakhotnik, “Accurate measurements of the Raman scattering coefficient and the depolarization ratio in liquid water,” Appl Opt. 52(11), 2503–2510 (2013). 32. N. Everall, T. Hahn, P. Matousek, A. W. Parker, and M. Towrie, “Picosecond Time-Resolved Raman Spectroscopy of Solids: Capabilities and Limitations for Fluorescence Rejection and the Influence of Diffuse Reflectance,” Appl. Spect. 55(12), 1701–1708 (2001). 33. H. S. Carslaw, and J. C. Jaeger, Conduction of Heat in Solids (Oxford University, 1959). 34. J. Wu, Y. Wang, L. Perelman, I. Itzkan, R. R. Dasari, and M. S. Feld, “Three-dimensional imaging of objects embedded in turbid media with fluorescence and Raman spectroscopy,” Appl. Spect. 34(18), 3425–3430 (1995). 35. A. Sassaroli, and F. Martelli, “Equivalence of four Monte Carlo methods for photon migration in turbid media,” J. Opt. Soc. Am. A 29(10), 2110–2117 (2012). 36. E. Alerstam, T. Svensson, and S. Andersson-Engels, “Parallel computing with graphics processing units for highspeed Monte Carlo simulation of photon migration,” J. Biomed. Opt. 13(6), 060504 (2008).


Introduction
Due to its excellent chemical specificity, Raman spectroscopy has been largely investigated in biomedical diagnosis and industrial applications [1].In most cases, a quasi-confocal approach was adopted, which limits the investigated volume to a maximum depth < 200 μm [2,3].For this reason, in the last decade, much interest was raised on the possibility to detect Raman signals well beneath the surface [4].The most popular approach was named Spacially Offset Raman Spectroscopy (SORS).SORS exploits the detection of Raman signals at increasing distances from the laser injection point to discriminate structures at different depth [4][5][6].Based on SORS, many applications have been proposed with huge industrial impact.Applications range from medical diagnostics (e.g.detection of cancer [7] or bone pathologies [8]), to security and forensic applications (e.g.identification of liquids or solids through plastic or glass containers [9]), and also to non-destructive food assessment [10] or analysis of paintings [11].Another approach was adopted for the analysis of pharmaceutical tablets, by exploiting a transmittance geometry [12,13].Similar concepts have also permitted to obtain tomographic Raman images through tissue simulating phantoms, biological specimens and, in vivo, on animals [14][15][16].
All approaches mentioned above are based on continuous wave (CW) illumination.A further different option is the adoption of a time-resolved (TR) injection and detection scheme.The TR scheme exploits the property that the more diffused photons are detected later on in time, the more they have penetrated deeply into the medium [17].A further advantage of this method is that the use of a time-gated scheme permits to reject a significant fraction of the contaminating fluorescence emission [18].The experimental feasibility of deep Raman spectroscopy using the TR approach was demonstrated using time-correlated single-photon counting [19], ultrafast (ps) Kerr-Gate spectroscopy [20] and fast gated camera [21][22][23].
It is important to note that the TR research direction, although possibly more informative than SORS, has to face with the complexity of the TR systems, which has hampered until nowadays a more widespread use.However, even if the translation of these methods and devices to Raman spectroscopy is not straightforward, a dramatic reduction of TR system complexity is expected by the use of novel compact photonics devices like Silicon Photomultipliers (SIPMs) [24].Thus, TR systems could definitevely increase their impact in the near future because of the great potentialities offered by this instrumentation to gain depth information [25,26].
In this framework, the development of rigorous mathematical algorithms, constituting one of the core elements of the TR systems, appears to be a fundamental issue.Moreover, theoretical models may also serve as a support for simulation studies in the development of novel measurement schemes, or in the understanding of the physics of Raman signal from its generation to its propagation and detection.For this reason, in the present work we address the rigorous modeling of TR diffused Raman based on the diffusion approximation (DA) to the Radiative Transfer Equation (RTE).
In summary, in the time domain exact analytical models describing propagation of Raman signals in diffusive media are still lacking.With the present work, we want to provide a general perspective of time domain diffuse Raman spectroscopy through its possible information content in terms of sensitiveness to the optical properties of the medium.We have developed a set of forward analytical solvers describing the TR Raman reflectance for a parallelepiped and a finite cylinder.The solvers have been obtained exploiting the DA to the RTE [27,28] and the Green's function theorem [28], and they have been systematically validated by "gold standard" Monte Carlo (MC) simulations for a wide set of physical conditions.Finally, the range of validity of an often used simplified heuristic model [6] for time-domain Raman signal is also discussed.

Theory and methods
In this section, we review the theory of Raman scattering within the framework of the diffusion equation (DE) and we obtain analytical solutions for the Raman TR reflectance.The analytical theory is also developed for a background medium hosting distributed fluorescent molecules.The MC method used to generate "gold standard" data for the whole study is also presented.

Preamble
The origin of Raman scattering [1,4,29] is quite different from the classical Tyndall scattering.Tyndall scattering depends on the micro-structures of materials and on their heterogeneities, while Raman scattering depends on the intrinsic properties of molecules.Further, Tyndall scattering is an elastic interaction, only affecting the direction of photon migration, while Raman scattering is an inelastic phenomenon affecting direction and wavelength of the re-emitted pho-ton.This means that the two phenomena are characterized by independent scattering coefficients, but, more relevant, it also implies that a couple of independent RTEs need to be introduced for dealing with two different wavelengths λ and λ e .Raman scattering can be considered "equivalent" to an absorption interaction since the scattered photon is no longer available at the excitation wavelength, and thus can be accounted as a loss factor in the energy balance at this wavelength.We follow this approach in developing the analytical theory here described.
Figure 1 shows the excitation light at λ (isotropic point source at r s ), and the Raman light at λ e , which is emitted from a scatter at r and collected at point r.The global, absorption and scattering coefficients of the medium originate from: 1) the background medium and; 2) the Raman scattering molecules.The Raman signal is affected by Tyndall multiple scattering in the background at both excitation, λ, and emission, λ e , wavelengths.The global optical parameters: absorption coefficient, scattering coefficient, reduced scattering coefficient, and the refractive index of the investigated medium at λ are denoted, respectively, as μ a , μ s , μ s , and n i , and at λ e as μ ae , μ se , μ se , and n ie .For the background medium the absorption coefficient, the scattering coefficient, and the reduced scattering coefficient at λ are denoted, respectively, as μ ab , μ sb , μ sb , and at λ e as μ abe , μ sbe , μ sbe .We assume that Raman scattering molecules are distributed inside the whole volume of the medium with scattering coefficient μ sR at λ, and μ sRe at λ e .We also assume that μ sRe = 0, excluding the possibility of a second Raman scattering event on the same photon [4] (see next section).Being Raman interactions described by a scattering coefficient, their effects at λ, at least within the DE, can only be modeled with an equivalent effective "absorption" term.Thus, since Raman scattering at λ determines an instantaneous re-emission of light at λ e , for the purpose of the energy balance at λ, Raman interactions are equivalent to absorption events, and therefore, we can then write μ a = μ ab + μ sR , μ s = μ sb and μ s = μ sb .For a summary of the notation used see the column denoted Raman in Table 1.This approach may appear to be not orthodox from a physical point of view, but, what it really matters is to perform a correct energy balance at λ and λ e , without neglecting terms influencing the photon migration.In the next sections the validity of the proposed modeling of Raman scattering will be further discussed and demonstrated by means of MC simulations.

General Raman model based on the time-dependent diffusion equation
Light propagation through biological tissues and other highly scattering media can be described with the time-dependent DE [27,28] due to the diffusive regime of propagation that is often established in these media.Therefore, the Raman signal can be described by two coupled diffusion equations for the time-dependent photon fluence rate, Φ(r, t) and Φ e (r, t), at λ and λ e : where D = 1/(3μ s ) and D e = 1/(3μ se ) are the diffusion coefficients at λ and λ e and v and v e are the speed of light in the medium at λ and λ e .The term q(r, t) is the real source term at λ, while q e (r, t) represents the virtual source term, at λ e , generated by Raman scatterers.To facilitate the reading, the dependence of Φ(.) and Φ e (.) on variables other than r and t will be explicitly written along the manuscript only when necessary.Equations ( 1) and ( 2) are coupled by the physical relationship existing between the excitation fluence rate, Φ(r, t), and the source term q e (r, t).As we have highlighted above, at λ the effect of Raman scattering is "equivalent" to an "absorption" interaction.Then, μ sR can be accounted in Eq. ( 1) as an absorption term so that μ a = μ ab + μ sR .The source term q e (r, t), at λ e , accounts for the number of photons generated at this wavelength per unit volume and time due to Raman scattering events at λ. Thus, using the fluence rate at λ, Φ(r s , r , μ a , μ s , n i , t , λ), the source term can be written as q e (r , t , λ e ) = μ sR Φ(r s , r , μ a , μ s , n i , t , λ). ( Assuming an isotropic source term of unitary strength q(r, t) = δ 3 (r − r s )δ(t), then Φ in Eq. ( 3) is the Green's function, G, of the problem and q e can be written as By substituting Eq. (4) in Eq. ( 2), the following equation at λ e is obtained: According to the Green's function theorem [28], the fluence rate at λ e , Φ e (r, t), due to the Raman scattering at λ is given by with G e the Green's function at λ e .We stress that the Green's functions G and G e are related to the same medium, so that the different notation is maintained only to keep the distinction between λ and λ e .Putting all explicitly, Equation (7) is the convolution in time of G with the same Green's function calculated at the site of the Raman scatterers.We remind that Eq. ( 7) is obtained under the assumption that only single scattering Raman events contribute to the solution.In principle, also multiple scattering Raman events can be possible.But, given the low values for the probability of Raman scattering (usually less than 10 −6 ) [30,31], their contribution to the calculation can be considered negligible.Moreover, in Eq. (7) it is implicitly assumed that the angular re-emission of Raman scattering is the same of Tyndall scattering since no distinctions on this point are made in the theory.Starting from Eq. ( 7) it is possible to derive the solution of the Raman signal for any desired geometry.In particular, Raman scatters can be homogeneously distributed inside the medium, but they can also occupy a small portion of the whole medium resulting as an inclusion.In this work we focus our analysis on geometries where the Raman scatters are homogeneously distributed inside the medium.In the previous literature we have an example of a heuristic model describing the Raman signal [29,32] obtained under the hypothesis that the optical properties of the background medium at excitation (λ) and emission (λ e ) wavelengths are the same.This model can be true or approximately true in some particular cases of photon migration.Here we want to provide an analytical justification of the heuristic model obtained by Everall et al. [29,32].Let's first re-write the general solution Φ(r, t) by invoking the very well known scaling relationship for the absorption coefficient [27,28] (Beer-Lambert's law) as The evaluation of the fluence rate at λ e can be obtained under the hypothesis /n e and also assuming the same scattering function of the medium at λ and λ e and for Tyndall and Raman scattering.This means that the Raman photons at λ e follow the same identical trajectories, and with the same probability, as they were Tyndall photons at λ and thus subjected to the optical properties of the medium at λ. Thus, the fluence rate at λ e due to Raman photons, Φ e (r, μ ae , μ se , t, λ e ), can be synthetically written as By using the property of Eq. ( 8) and the described constraints on the optical parameters, the right term of Eq. ( 9) can be rigorously expressed as i.e. the fluence rate of the whole photons reaching the detector (photons at λ or λ e are subjected to the same optical parameters) minus the fluence rate of the photons that did not undergo a Raman scattering.
Given the typical low values of μ sR (of the order of 10 −6 mm −1 ), Eq. ( 10) can be usually further simplified as The same relation can be obtained for the TR reflectance.On this ground the heuristic model for the TR reflectance at λ e , R eHeur , can be written as with ρ source detector distance and R the standard solution of the DE at λ.The solution of the DE for the TR reflectance for the slab obtained with the extrapolated boundary condition and Fick's law can be used for R (see for instance Eq. (4.27) in [28]).Equations ( 11) and ( 12) are both the desired Everall et al. model [29,32] obtained under the simplified conditions.Thus, it is possible to directly derive the Green's function of the Raman signal from the Green's function of the background medium and the Raman signals can be obtained by multiplying the diffuse reflectance signal by μ sR vt.The physical interpretation of this simplified model is that the generated Raman photons at the emission wavelength continue their migration in the medium as they were at the excitation wavelength.We can also state that the longer is the photon time-offlight, the higher is the probability of Raman emission.This heuristic model can only be applied for homogeneous media.

Time-resolved fluence rate
The solution for the parallelepiped (Fig. 2) for the fluence rate is obtained using the Green's function calculated with the eigenfunction method.The TR Green function G(r, r , t, t ) for the fluence rate, making use of the extrapolated boundary conditions (EBC) at the boundaries of the parallelepiped, can be written as [28,33] and The coefficients of Eqs. ( 14) and ( 15) will be used in all the next series expansions implemented for the parallelepiped.The coefficient A, that accounts for the effects of Fresnel reflections at the parallelepiped boundaries, is a function of the refractive index of the internal (n i ) and external (n o ) medium (A = 1 when n i = n o see [28]).The notation A e will be used when the coefficient A pertains to n ie .Making use of Eq. ( 13) and of the ortho-normality property of the eigenfunctions, the solution of Eq. ( 7) for the TR fluence rate, assuming an isotropic delta source of unitary strength placed at (0, 0, z s = 1/μ s ) used to model an external pencil beam of unitary strength impinging onto the parallelepiped at (0, 0, 0), i.e. q(r, t) = δ(x)δ(y)δ(z − z s )δ(t), becomes For obtaining this equation we have assumed, as it is usual, that the extrapolated lengths depend only slightly on λ, and thus, L xe L x , L ye L y , and The same approximation is used throughout this paper.Once more we remind that μ a = μ ab + μ sR and μ ae = μ abe .A similar solution to Eq. ( 16) for the fluence rate is obtained in appendix A also for the finite cylinder.
The TR Raman reflectance from the external surface of the parallelepiped at any point r of the boundary can be obtained using Fick's law or a hybrid approach denoted EBPC [28].For both the approaches the Raman fluence rate, Φ eRaman , (Eq. ( 16)) is first calculated using the EBC, then Fick's law or the partial current boundary condition (PCBC) is applied to Eq. ( 16) for calculating the emerging reflectance from the medium [28].

Time-resolved reflectance obtained with the EBPC
According to the EBPC the TR Raman reflectance is [28] R eRamanEBPC (x, y, t) = Φ e Raman (x,y,z=0,t ) By substituting Eq. ( 16) in Eq. ( 17) we obtain the TR Raman reflectance as 2.4.3.Time-resolved reflectance with Fick's law The TR Raman reflectance from the medium can be also calculated by using the classical Fick's law as [28] R eRamanFick (x, y, t) = D e ∂Φ e Raman (x,y,z=0,t ) ∂z .
By inserting Eq. ( 16) in Eq. ( 19) the TR Raman reflectance becomes The two expressions of R eRamanEBPC and R eRamanFick show only slight differences related to the way the flux is calculated [28].In the Results section comparisons with the results of MC simulations show that within the limits of the DA the two formulas are practically equivalent.

Time-resolved fluence rate
It may happen that a fluorescence signal survives at the detector's site when Raman spectroscopic measurements are carried out or, alternatively, fluorescence and Raman signals can be indeed used together for imaging and/or spectroscopic purposes [34].For this reason, we generalize the previous solution in the event that, together with Raman scatterers, fluorescent molecules are also uniformly distributed inside the medium generating a background fluorescence at λ e (observed wavelength).The fluorescence signal, as the Raman scattered light, is affected by multiple Tyndall scattering at both λ and λ e .We assume that a fluorophore is uniformly distributed inside the background medium with absorption coefficient μ a f at λ, and μ a fe = 0 at λ e .As in the previous section, we also assume that Raman scatterers are distributed inside the whole volume of the medium with scattering coefficient μ sR at λ, and μ sRe = 0 at λ e .The global optical parameters: absorption coefficient, reduced scattering coefficient, and the refractive index of the medium are denoted as μ a , μ s , and n at λ, and as μ ae , μ se , and n e at λ e , respectively.Therefore, we implement the DE solution for a medium having the global optical properties, μ a = μ ab + μ sR + μ a f , μ s = μ sb at λ and, μ a = μ abe and μ se = μ sbe at λ e , respectively.For this purpose we can still use Eqs.( 1) and (2) provided the optical properties are those shown in Table 1 for the case of Raman+Fluorescence.
Again, the DEs at λ and λ e are coupled by the excitation fluence rate, Φ, which determines the source term q e (r, t) of the DE at λ e .However, the fluorescent component of the source term q e (r, t) is more complicated than the Raman component.In fact, not all the fluorescence interactions at λ produce a photon at λ e .Moreover, fluorescent photons at λ e , generated at time t, may be derived from fluorescence interactions occurring at a time t different from t (t < t).For these reasons, it is necessary to define: 1) a quantum efficiency η e and; 2) a probability density function p Fluo (t | t ).The parameter η e is the quantum efficiency of the fluorophore in the wavelength range Δλ e around λ e , i.e. the ratio of the number of photons emitted in Δλ e to the number of absorbed photons.The function p Fluo (t | t ) is expressed as ( p Fluo (t | t )dt = 1) and describes the probability density that an interaction at t will generate a photon at t.The function Θ(.) is the Heaviside function.The parameter τ is the fluorescence lifetime of the fluorophore (i.e. the average time that the fluorophore spends in the excited state prior to returning to the ground state).In practice, the production of fluorescent photons must be suitably "weighted" by η e and p Fluo (t | t ).Therefore, in analogy with Sec.
2.2 the coupling of the DE at λ and λ e is given by the following source term that combines fluorescent and Raman emitted photons where G is the Green's function for the fluence rate at λ.By substituting Eq. (22) in the DE at λ e , and making use of the Green's function theorem, the fluence rate at the emission wavelength, Φ e (r, t), due to the fluorescent photons generated into the interval dλ e , is with G e as the Green's function of the medium for the fluence rate at λ e .The solution of Eq. ( 23) is given by two terms, one, Φ eFluo (r, t), given by fluorescent photons and the other, Φ eRaman (r, t), given by Raman photons.The term Φ eRaman is the one calculated in the previous section with Eq. ( 16).For the term Φ eFluo we need to perform the integral of Eq. ( 23) for the fluorescence source term.For the parallelepiped the TR solution is obtained by making use of the TR Green's function given by Eq. ( 13) and applying the orthonormality property of the eigenfunctions.The solution given by Eq. ( 23) for the TR fluence rate, Φ eFluo , assuming an isotropic delta source of unitary strength placed at (0, 0, z s = 1/μ s ), i.e.
Table 1.Summary of the notations used for the optical properties of the medium at λ and λ e in the case of a pure Raman model (column Raman) or a hybrid Raman and fluorescence model (column Raman+Fluorescence).The subscript b denotes the background, the subscript R denotes the Raman scattering and the subscript f denotes the fluorescence.In this work, we have assumed μ sRe = 0 and μ a fe = 0.The optical parameters appearing in column Raman have to be used with Eqs. ( 18), ( 20) and (39), while optical parameters appearing in column Raman+Fluorescence have to be used with Eqs. ( 24), ( 28) and (32).

Raman
We note that, for the evaluation of the Eq. ( 23) we have previously calculated the first term of q e (r , t ), q eFluo (r , t ), (Eq.( 22)) as with G given by Eq. (13).

Time-resolved reflectance with Fick's law
The time resolved reflectance can be also calculated by using Fick's law as with R eRamanFick given by Eq. ( 20) and R eFluoFick given by (28)

Improved numerical calculation
The formula obtained above for Φ eFluo suffers from slow computational convergence so that usually many terms of the series are needed for a correct evaluation of the fluence rate.A simple way to speed-up the calculation is achieved by dividing the calculation of Φ eFluo into two steps.At first is calculated the fluence rate, Φ e0Fluo , describing fluorescence photons promptly emitted (i.e.Eq. ( 24) for τ → 0).The series of Φ e0Fluo converges much faster than Φ eFluo .The fluence rate Φ e0Fluo is simply given by ) Then, the calculation of the whole term Φ eFluo is obtained by taking into account that fluorescent photons are not promptly re-emitted.Therefore, Φ eFluo may be expressed as the convolution of Φ e0Fluo with the exponential decay of the fluorophore, i.e.
This calculation in two steps is computationally mush faster (even ten times) than the direct calculation of Eq. ( 24).Moreover, the implementation of Eq. ( 30) can be easily modified in order to account for more complicated multi-exponential decays that sometimes are required to describe the re-mission of some specific fluorophores.Similarly, the improved numerical calculation can be obtained for the fluorescence reflectance term as When the reflectance is calculated with the EBPC, R e0Fluo is obtained by dividing Φ e0Fluo by a factor 2A e according to Eq. ( 26).When the reflectance is calculated with Fick's law, R e0Fluo is obtained as

Monte Carlo simulations
It is useful, for the comprehension and validation of the analytical approach proposed in the previous sections, to address and simulate Raman scattering, also in the presence of background fluorescence, using a Monte Carlo method without introducing approximations in the modelling of photon transport [28,29].Following the notations of Sec.2.5 we define an interaction coefficient, μ eλ , at the excitation wavelength as: where the right hand side is the sum of the four coefficients of: absorption in the background, Tyndall scattering in the background, Raman scattering and fluorescence excitation.A random number w 1 uniformly distributed in the interval [0, 1] is generated.This number is related to photon path, , by the following relation [35]: Once Eq. ( 34) is inverted to obtain we have to decide which interaction the photon will undergo.To sample the probability of the four different events we extract a second random number w 2 uniformly distributed in [0, 1], so that: the photon is Raman scattered, otherwise if w 2 ≤ 1 the photon possibly excites fluorescence.
When the photon is absorbed its propagation is stopped; when it is Tyndall scattered, the photon is tagged as "λ" (the wavelength remains λ), the scattering angle is sampled [35] and it moves of a step-length equal to ; when it is Raman scattered the photon is tagged as "λ e " (from now on the photon can only have Tyndall scattering events at λ e ), an isotropic scattering angle is sampled and it moves of a step-length equal to ; when it possibly excites fluorescence a further random number w 3 ∈ [0, 1] is extracted and, if w 3 < η, the photon is tagged again as "λ e " and it propagates with a step-length (from now on the photon can only have Tyndall scattering events at λ e ), otherwise the photon is killed.Furthermore, if the photon excites fluorescence, a temporal delay is extracted using a new random number w 4 ∈ [0, 1] sampling the exponential decay time (Eq.( 21)): The photons tagged as "λ e " represent photons at the emission wavelength then, assuming that no further Raman scattering events are undergone by photons, and they propagate as described above using a new interaction coefficient μ eλ e given by: At the beginning of the simulation it is possible to select the kind of interaction (Raman, fluorescence, Raman+fluorescence or Tyndall) and the number of photons that, after surviving at the detector site, need to be saved and stored by the simulation.In particular, the path-length of the tagged photons arriving at the detector site is saved and then a time-of-flight histogram is created so the TR reflectance can be re-constructed from the detected photons packets.

Validation of the DE solutions by comparisons with MC simulations
The analytical solutions derived in Sec.2.5 have been compared with the results of MC simulations based on the sampling rules described in Sec.2.6.The code was accelerated by means of a Graphical Processing Unit (GPU) and has been written adapting the CUDA-based code published by Alerstam et al. [36] a laterally infinite slab.A large number of comparisons between analytical models and MC simulations have been performed to validate the formulae of Sec.2.5 and those shown in this section are only a representative sample of the whole work done.Simulated reflectance data have been generated with the optical parameters appearing in the caption of Fig. 3.All the simulations have been run until 10 7 photons reached the detector.The analytical models have been computed using the formulae presented in Sec.2.5 applied to a 100x100x100 mm 3 diffusive cube.The boundary effects due to photons escaping the lateral faces of the cube are negligible for the values of the optical properties used in this section.
Although the probability of generating a Raman photon is low, with typical values of about 10 −6 , for the validation of the analytical models we have assumed unrealistic values of this probability with values around 10 −3 .This was obtained by selecting μ sR = 10 −3 mm −1 .This choice, needed for reducing the computation time of the MC simulations, does not affect the validation of the analytical models.In fact, the accuracy of the analytical models, provided that μ sR μ sb , is not affected by the value of μ sR considered.In Fig. 3 four examples of comparison of the analytical models (EBPC and Fick solutions) with the results MC simulations for the TR Raman-Fluorescence reflectance are shown.Starting from the first simulation (Fig. 3(a)-3(b)), the others have been obtained by varying one parameter at a time: μ a f (Fig. 3(c)-3(d)), another with a different ρ (Fig. 3(e)-3(f)) and the remaining one with a different μ ab (Fig. 3(g)-3(h)).We note that the optical properties used for the background medium are typical of several biological tissues.It is also worth to observe the excellent agreement between analytical models and MC results.
The slight differences between Fick and EBPC model are typical of the DA and, as for many other solutions based on the DE, they are mainly affecting the early times of the TR reflectance.The Fick and EBPC solutions slightly differ on the peak of the signal.In particular, EBPC overestimates the photon fluence rate.At late times and in absence of fluorescence the two formulae converge (Fig. 3(b)).
In the presence of a background fluorescence (Fig. 3(c)), the agreement between Fick's law formula and MC is still very good and slightly worse for the EBPC formula because the fluence rate overestimation on the peak is spread out at late times.This is due to the convolution operation appearing in Eq. ( 31).This effect is clearly noticeable in Fig. 3(d) by plotting the ratio of analytical solutions and MC.But we can note that this overestimation of the signal has a mild effect on the shape of the TR reflectance generated with the EBPC.

Dependence of the Raman signal on the optical properties
The factor μ sR appearing in Eqs.(18) and (20) clearly shows that the TR Raman reflectance is proportional to the Raman scattering coefficient.In practice, μ sR is just a simple amplitude factor and does not influence the shape of the temporal profile.However, Eqs. ( 18) and ( 20) also show that the TR Raman reflectance depends on absorption and reduced scattering coefficients in a more complex manner.In Fig. 4 we have reported this feature by displaying the TR Raman reflectance with: 1) varying μ abe , fixed μ ab and fixed μ sb = μ sbe and; 2) varying μ sb = μ sbe and fixed μ ab = μ abe .The values of the optical parameters appear in the caption of Fig. standard TR diffuse optical spectroscopy, with a decrease of signal more prominent at late times for increased absorption (Fig. 4(a)), and at early times for increased scattering (Fig. 4(b)).In conclusion, Fig. 4 show that the TR Raman signal carries information not only of the Raman phenomenon, but also on the background properties of the medium.

Validity range of the heuristic model
As already demonstrated in Sec.2.3 the heuristic model proposed in [29,32] by Everall et al. is actually a simplified model derived under the constraints of μ sb = μ sbe , μ ab = μ abe and n i = n e .Here we want to study the range of validity of the heuristic model by means of the forward solvers discussed and validated in the present work.More specifically, we study the effect given by the release of the above constraints.Figure 5 shows an example of the TR Raman signal (Eq.( 20)) and of the heuristic model (Eq.( 12)) for a case with μ sb μ sbe (see figure caption for the values of the optical parameters).The figure shows that the heuristic model fails in describing the correct behavior.
In Fig. 6 the differences observed in Fig. 5 are studied in a larger number of cases in order to provide a more complete view of the validity range of Eq. (20).To this aim we define a relative error factor, (t), comparing the heuristic model (Eq.( 12)) with the rigorous solution (Eq. ( 20)): (t) = 1 − R e Heur (ρ,t ) R e RamanFick (ρ,t ) .
The error (t) is depicted in Fig. 6 in per cent for different values of μ abe (Fig. 6(a)) or of μ sbe (Fig. 6(b)) while keeping both μ ab and μ sb fixed.The relative error is null when μ ab and μ abe have the same values, while it rapidly deviates with huge errors whenever μ abe is significantly altered from μ ab .For example, in biological tissues, for Raman shifts in the range of [30,50] nm a change of μ abe in the range of [0.005, 0.01] mm −1 can be expected for an excitation wavelength in the range of [700, 800] nm.In case we have a change of scattering between λ and λ e , the heuristic model is substantially valid at late times, but yields quite inaccurate predictions at early times.Thus, Eq. ( 12) can be used as a first-order approximation to understand the physics of diffuse Raman photons, but a more rigorous approach is needed when the absorption properties of the medium are not constant over the spectral range of interest.

Conclusions
We have presented, within the framework of the DE, a set of analytical forward solvers for the calculation of the Raman TR reflectance from a homogeneous parallelepiped and a homogeneous cylinder, provided the Raman scatterers are uniformly distributed in the whole medium.The obtained formulas have been validated by means of comparisons with the results of "gold standard" MC simulations.In the proposed analytical solutions we have also included the option to account for the influence of a background fluorescence, on the Raman signal, generated by a fluorophore uniformly distributed inside the medium.Further, we have provided an exact theoretical justification of a well known heuristic model; widely used for the study of the Raman signal when the optical properties of the medium show weak variations between the excitation and the emission wavelength.Thanks to the present results, it has been possible to carefully study the range of validity of the heuristic model.The comparisons between analytical models and MC simulations of Sec. 3 have shown that in the diffusive regime the analytical solutions (Eqs.( 18) and ( 20)) of the DE for the TR Raman reflectance show indistinguishable results compared to the MC simulations.The differences between Fick and EBPC model are compatible with the diffusion approximation and both can be used to describe the Raman signal with an error of few per cent.Similar conclusions can be extended to the effect of a background fluorescence that can be accurately calculated by the DE analytical solutions (Eqs.( 28) and ( 26)) with an error of few per cent.The computation time of such models is of the order of seconds, while the computation time of the corresponding MC simulation, although the simulation was accelerated by means of a GPU, was of the order of hours.Moreover, when the MC simulations are carried for realistic values of the Raman scattering coefficient (of the order of 10 −6 mm −1 ) the computation time can be even thousand time larger, becoming prohibitively high.The comparisons with the results of MC simulations have shown that the heuristic model (Eqs.(11) and ( 12)) has significant deficiencies when the absorp- tion coefficient varies between the excitation and the emission wavelength, while the changes of scattering affect the accuracy of the model only at early times.In practice, the proposed forward solvers based on the DE are suitable tools for studying the TR Raman signal in diffusion conditions, while the heuristic model can only provide an approximate solution that is consistent with the DE when the optical properties at excitation and emission wavelengths are coincident.
In conclusion, the proposed forward solvers can be a useful tool for exploring the information content encoded in the TR Raman measurements.The use of these solvers will allow one to depict, in a wide perspective, possible new applications derived from time-domain Raman measurements.

A. Appendix: DE solution for a finite cylinder
We provide here the soluton of the DE for the Raman signal for the finite cylinder.This solution for the fluence rate is obtained with the same procedure used for the parallelepiped in Sec.2.4, provided the radial eigenfunctions are taken equal to the Bessel functions of the first kind of order zero J 0 .Given a finite cylinder (see Fig. 7) of thickness s, radius a, with the z axis along the main axis of the cylinder, and with ρ the distance between the position r and the z axis, the Green's function with the EBC of the DE for the finite cylinder of Fig. 7 can be written as [28,33] G (r, r', t) = 2 v with a = a + 2AD, s = s + 4AD, K n = nπ/s , λ l roots of J 0 (λ l a ) = 0, J 0 Bessel function of the first kind of order zero and J 1 Bessel function of the first kind of order one.Making use of Eq. ( 38) and of the ortho-normality property of the eigenfunctions, the solution of Eq. ( 7) for the TR fluence rate, assuming an isotropic delta source of unitary strength placed at (0, with a = a+2AD, s = s+4AD, K n = nπ/s , λ l roots of J 0 (λ l a ) = 0.The above equation has been obtained with the approximations a a e and s s e .From Eq. (39) the TR reflectance can be obtained with Eq. ( 17) or with Eq. (19).
We stress that the presented solution can be useful to speed-up the calculations since by using Eq.(39) the computation time can be reduced of a factor ten compared to Eq. ( 16).

Fig. 1 .
Fig. 1.Diagram of a photon path trajectory in a x, z plane projection from the source (red circle) to the Raman scatter (dashed line) and from the Raman scatter to the detector (continuous line).Green circles are Tyndall interactions.The yellow circle is a Raman interaction.