1. 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–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–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–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.

2. 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.

2.1. 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 photon. 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.

 

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.

Download Full Size | PDF

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 μ_afe = 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).

View Table

2.2. 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:

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

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⁻⁶) [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.

2.3. Simplified heuristic model in case μ′_sb = μ′_sbe, μ_ab = μ_abe and n_i = n_e

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

Given the typical low values of μ_sR (of the order of 10⁻⁶ mm⁻¹), Eq. (10) can be usually further simplified as

2.4. Solution for the Raman signal in a parallelepiped

2.4.1. 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]

 

Fig. 2 Schematic of the parallelepiped.

Download Full Size | PDF

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].

2.4.2. Time-resolved reflectance obtained with the EBPC

According to the EBPC the TR Raman reflectance is [28]

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]

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.

2.5. Solution for the Raman signal with a background fluorescence

2.5.1. 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 μ_af at λ, and μ_afe = 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 + μ_af, μ′_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

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

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

2.5.2. Time-resolved reflectance with EBPC

The TR reflectance can be calculated using the EBPC, similarly to Eq. (17), as

2.5.3. Time-resolved reflectance with Fick’s law

The time resolved reflectance can be also calculated by using Fick’s law as

2.5.4. 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 re-emitted (i.e. Eq. (24) for τ → 0). The series of Φ_e0Fluo converges much faster than Φ_eFluo. The fluence rate Φ_e0Fluo is simply given by

Similarly, the improved numerical calculation can be obtained for the fluorescence reflectance term as

2.6. 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:

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₃ ∈ [0, 1] is extracted and, if w₃ < η, 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₄ ∈ [0, 1] sampling the exponential decay time (Eq. (21)):

3. Results

3.1. 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] for 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⁷ photons reached the detector. The analytical models have been computed using the formulae presented in Sec. 2.5 applied to a 100×100×100 mm³ 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.

 

Fig. 3 Comparison for the TR Raman reflectance between analytical formulae and MC in absence (a, e, g) and in presence (c) of background fluorescence. For all the figure panels μ′_sb = μ′_sbe = 1, μ_abe = 7 · 10⁻³ mm⁻¹, μ_sR = 10⁻³ mm⁻¹, n_i = n_ie = 1.4, n_o = n_oe = 1, g = 0. Specifically: (a) μ_ab = 5 · 10⁻³ mm⁻¹, μ_af = 0, ρ = 20 mm; (c) μ_ab = 5 · 10⁻³ mm⁻¹, μ_af = 7 · 10⁻³ mm⁻¹, η_e=1 and τ = 2 ns, ρ = 20 mm; (e) μ_ab = 5 · 10⁻³ mm⁻¹, μ_af = 0, ρ = 10 mm; (g) μ_ab = 10⁻² mm⁻¹, μ_af = 0, ρ = 20 mm.

Download Full Size | PDF

Although the probability of generating a Raman photon is low, with typical values of about 10⁻⁶, for the validation of the analytical models we have assumed unrealistic values of this probability with values around 10⁻³. This was obtained by selecting μ_sR = 10⁻³ mm⁻¹. 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: μ_af (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.

3.2. 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. 4. The ranges of the background optical properties in Fig 4 include those of biological media in the wavelength range of [700, 900] nm. The general trend is similar to what observed in 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.

 

Fig. 4 Dependence of the temporal profile of Raman photons for different values of μ_abe and μ′_sb. The common parameters were ρ = 10 mm, n_i = n_ie = 1.4, n_o = n_oe = 1 and μ_sR = 10⁻⁶ mm⁻¹, μ_ab = 0.01 mm⁻¹. Specifically: (a) μ_abe ∈ [0, 0.035] mm⁻¹ and μ′_sb = μ′_sbe = 1mm⁻¹ and; (b) μ′_sb = μ′_sbe ∈ [0.5, 2] mm⁻¹ and μ_abe = 0.01 mm⁻¹.

Download Full Size | PDF

3.3. 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.

 

Fig. 5 Comparison of the temporal profile of the DE Raman signal (Eq. (20)), and of the heuristic model (Eq. (12)). The common parameters were: ρ = 10 mm, n_i = n_ie = 1.4, n_o = n_oe = 1, μ_sR = 10⁻⁶ mm⁻¹, μ′_sb = μ′_sbe = 1mm⁻¹, and μ_ab = 0.01 mm⁻¹.

Download Full Size | PDF

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)):

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⁻¹ 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.

 

Fig. 6 Relative per cent error of the heuristic model (Eq. (12)) as compared to the DE solver (Eq. (20)). The common parameters were ρ = 10 mm, n_i = n_ie = 1.4, n_o = n_oe = 1, μ_sR = 10⁻⁶ mm⁻¹, μ_ab = 0.01 mm⁻¹ and μ′_sb = 1mm⁻¹. Specifically: (a) μ_abe ∈ [0, 0.02] mm⁻¹ and μ′_sbe = 1mm⁻¹ and; (b) μ′_sbe ∈ [0.5, 1] mm⁻¹ and μ_abe = 0.01 mm⁻¹.

Download Full Size | PDF

4. 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⁻⁶ mm⁻¹) 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 absorption 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₀. 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]

(38)$$\begin{array}{ll}G(\mathbf{r},{\mathbf{r}}^{\prime },t)\hfill & =2\frac{v}{\pi a^{\prime 2}{s}^{\prime }}\sum _{l=1}^{\infty }\sum _{n=1}^{\infty }\frac{J_0(\rho {\lambda }_l)J_0({\rho }^{\prime }{\lambda }_l)}{J_1^2(a^{\prime }{\lambda }_l)}\text{sin}\left[K_n(z+2AD)\right]\text{sin}\left[K_n(z^{\prime }+2AD)\right]\hfill \\ \hfill & \text{exp}\left[-{\lambda }_l^{2}Dv(t-t^{\prime })\right]\text{exp}\left[-K_n^{2}Dv(t-t^{\prime })\right]\text{exp}\left[-{\mu }_{a}v(t-t^{\prime })\right],\hfill \end{array}$$

with a′ = a + 2AD, s′ = s + 4AD, K_n = nπ/s′, λ_l roots of J₀(λ_la′) = 0, J₀ Bessel function of the first kind of order zero and J₁ 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, 0, z_s = 1/μ′_s) used to model an external pencil beam of unitary strength impinging onto the cylinder at (0, 0, 0), i.e. q(r, t) = δ(x) δ(y) δ(z − z_s), δ(t), becomes

(39)$$\begin{array}{l}{\mathrm{\Phi }}_{\mathit{eRamanCyl}}(\mathbf{r},t)=2\frac{{\mu }_{sR}v{v}_e}{\pi a^{\prime 2}{s}^{\prime }}\sum _{l=1}^{\infty }\sum _{n=1}^{\infty }\frac{J_0(\rho {\lambda }_l)}{J_1^2(a^{\prime }{\lambda }_l)}\text{sin}\left[K_n(z+2A_e{D}_e)\right]\text{sin}\left[K_n(z_s+2AD)\right]\\ \times {\left[\left(D_e{v}_e-Dv\right)\left({\lambda }_l^2+K_n^2\right)+\left({\mu }_{ae}{v}_e-{\mu }_{a}v\right)\right]}^{-1}\\ \times \left\{\text{exp}\left[-\left({\lambda }_l^2+K_n^2\right)Dvt-{\mu }_{a}vt\right]-\text{exp}\left[-\left({\lambda }_l^2+K_n^2\right)D_e{v}_{e}t-{\mu }_{ae}{v}_{e}t\right]\right\},\end{array}$$

with a′ = a + 2AD, s′ = s + 4AD, K_n = nπ/s′, λ_l roots of J₀(λ_la′) = 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).

 

Fig. 7 Schematic of a diffusive cylinder.

Download Full Size | PDF

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).

Acknowledgments

We acknowledge partial support from OILTEBIA (Optical Imaging and Laser TEchniques for BIomedical Applications, ITN) under Grant 317526. The authors wish to thank Prof. Giovanni Zaccanti for his useful suggestions.

References and links

1. P. Matousek and M. D. Morris, Emerging Raman Applications and Techniques in Biomedical and Pharmaceutical Fields (Springer, 2010). [CrossRef]  

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). [CrossRef]   [PubMed]  

3. N. J. Everall, “Confocal Raman microscopy: common errors and artefacts.,” Analyst 135(10), 2512–2522 (2010). [CrossRef]   [PubMed]  

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). [CrossRef]  

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. RamanSpectrosc. 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). [CrossRef]   [PubMed]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]  

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). [CrossRef]  

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). [CrossRef]   [PubMed]  

13. P. Matousek and A. W. Parker, “Bulk Raman Analysis of Pharmaceutical Tablets,” Appl. Spectrosc. 60(12), 1353–1357 (2006). [CrossRef]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]  

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). [CrossRef]   [PubMed]  

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 time-domain diffuse optics,” Opt. Express 23(11), 13937–13946 (2015). [CrossRef]  

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). [CrossRef]  

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). [CrossRef]  

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). [CrossRef]  

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). [CrossRef]  

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). [CrossRef]  

30. R. A. Desiderio, “Application of the Raman scattering coefficient of water to calculations in marine optics,” Appl Opt. 39(12), 1893–1894 (2000). [CrossRef]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]  

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). [CrossRef]  

36. E. Alerstam, T. Svensson, and S. Andersson-Engels, “Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration,” J. Biomed. Opt. 13(6), 060504 (2008). [CrossRef]