1. Introduction

The radiative transfer equation (RTE) provides a full description of light propagation in diffusive media like biological tissue in the red and near infrared spectral range. Unfortunately the RTE is an integro-differential equation and the retrieval of solutions is an extremely computationally expensive process [1]. Solutions of the RTE are usually based on numerical methods like the finite element method [1, 2, 3] or other methods like the discrete ordinates method [1, 4], the spherical harmonics method [1], the integral transport method [1] and the path integral formalism [5]. Among the numerical procedures there are also statistical methods like the Monte Carlo (MC) that is largely used to reconstruct solutions of the RTE [1, 6]. To our knowledge, no analytical (closed-form) solutions of the RTE are available [1] and simpler approximated models are usually sought [2, 3]. The most widely and successfully used approach makes use of the diffusion equation (DE) to yield a variety of solutions capable of providing a modelling of light propagation in homogeneous, layered and perturbed diffusive media, both in the steady state and in the time-domain. Solutions of the DE are available for the slab and for the semi-infinite medium that represent typical geometrical schematics for most of the experiments carried out in biomedical optics applications. In the last decade it has been possible to develop noninvasive tissue spectroscopy and imaging for diagnostic and therapeutic applications (e.g. optical mammography, brain imaging, photodynamic therapy) due to the availability of these DE based easy-to-use tools for light dosimetry and for the assessment of tissue optical properties [7].

However, the investigations based on the DE are subjected to the intrinsic approximations of this theory. In particular for continuous wave instrumentation the basic configuration needs a source-detector separation ρ, typically in the range 20–40 mm, because photons received at shorter distances are not well described by the DE [6]. Moreover, in the time-domain early received photons must be disregarded because the validity of the DE is limited to photons that, in their migration, undergo many scattering events [8]. Conversely, by exploiting an RTE solution it can be foreseen to overcome the limits of the DE and to immediately obtain a twofold improvement: (i) the error determined by approximations of the theory would be drastically reduced, therefore modelling by simulations would be more accurate; (ii) measurements carried out in conditions where the use of the DE is strongly limited (like small volume samples, high absorption and/or low scattering media, short source-receiver distances and early times) could be more easily interpreted. For these reasons the availability of analytical solutions of the RTE would have generally a strong impact on many practical applications in the biomedical optics [9, 10]. A first example is the optical molecular imaging of small animals where the use of quantitative methods to estimate tumor growth or the effect of drugs is hampered by the rough approximation of the DE [11]. Further, the possibility of performing measurements at null or small source-receiver distance has been proposed as a new imaging modality offering improved contrast, sensitivity and spatial resolution and as a novel approach to local spectroscopy of tissue [8]. Finally, weakly diffusive media of interest can be found in the biomedical field (e.g. blood, due to the high absorption) and in the industrial field (e.g. vegetable oil, due to the weak scattering). In general the applications of RTE solutions are addressed to several research fields where a transport process is involved.

Ten years ago Paasschens [12] proposed an almost exact solution of the time dependent RTE in the infinite medium geometry and for isotropic scattering. Paasschens, using a path-integral method, gave explicit time domain solutions of the RTE for two and four dimensions, while in three dimensions the solution was given in terms of its Fourier transform. The time domain solution in three dimensions was obtained by an interpolation between the two exact solutions found for two and four dimensions. It is not surprising that Paasschens’ paper received scarce attention within the biomedical optics community. The solution proposed by Paasschens is too far from being a workable solution. The possibility of modelling a sample as an infinite medium is in fact rarely encountered in real applications.

In this paper, starting from the work of Paasschens [12] and following a heuristic approach, an analytical solution of the time dependent RTE for the semi-infinite medium is derived for the first time. The solution is obtained in the simple case of isotropic scattering and point-like isotropic source. The accuracy of the proposed solution is validated by a comparison with results of MC simulations. Moreover, for the more realistic case of anisotropic scattering and pencil beam source, the derived solution still shows an excellent behavior. The proposed solution can be directly used in applications where a reflectance approach is used to noninvasively characterize a sample. Also it represents a first step to derive analytical solutions of the RTE for other geometries (e.g. slab).

2. Theory

Our derivation starts from the interpolating analytical formula for the time-dependent fluence rate [12], Φ_i(r, t), emitted by a pulsed point-like isotropic source (i.e. source=δ(r)δ(t)) in an infinite non-absorbing medium

where r is the distance from the source, Θ(x) is the step function (zero for x<0 and 1 for x>0), µ_s, µ′_s s and v are respectively the scattering coefficient, the reduced scattering coefficient and the speed of light in the medium, while the expression for the function G is provided as the following

At early times (small values of x) the first expression of Eq. (2) has to be preferred to the second one because it is more accurate, while for late times they tend to give the same values. The first term of Eq. (1) represents the ballistic peak, while the second term is the scattered component. The accuracy of Eq. (1) is better than 2% [12]. In the second term of Eq. (1), differently from the paper by Paasschens [12], we introduced µ′_s instead of µ_s. Since µ′_s=µ_s(1-g), where the symbol, g, denotes the anisotropy factor, for isotropic scattering (g=0) the two formalisms are identical. Furthermore, our approach would be useful to treat anisotropic scattering (g≠0) for which Eq. (2) can still be heuristically used. The effect of absorption can be introduced, according to the RTE, by multiplying the above time-domain solution obtained for the case of zero absorption by exp(-µ_avt).

When the semi-infinite geometry is considered the main problem to be addressed for obtaining the Green’s function is the modeling of the boundary condition with the external non scattering region. The boundary condition for the RTE and for an isotropic source embedded in the medium at a depth z ₀ would involve the radiance, I(r, ŝ, t), emerging from the medium along any direction ŝ. If the refractive indexes of both media are the same (as it is assumed in our derivation), the boundary condition results: I(r, ŝ, t)=0, for every point r≡(x, y, z) at the external boundary and for all directions ŝ inwardly directed. Although rigorous, this approach cannot be implemented to solve the RTE with an analytical method. In order to overcome this problem the same boundary condition used for the DE [13], usually named extrapolated boundary condition (EBC), is here heuristically used for the RTE. More in details, we derived the Green’s function by using the method of images [14]. The locations, r ⁺and r ^-, of the real source (positive) and of the image source (negative) are shown in Fig. 1. In this way the fluence rate will vanish at z=-z_e, which identifies the extrapolated boundary, and it is composed by two terms of fluence from an infinite medium, Φ_i: one positive from the source at r ⁺and the other negative from the source at r ^-. This heuristic approach leads to the following expression for the fluence rate in a semi-infinite medium, Φ_si(x, y, z, z ₀, t), with refractive index matched at the boundary

According to Ref. [13], we assume an extrapolated distance z_e=2/3µ′_s. Note that it differs in about 6%from the extrapolated end point, d_e, of the Milne problem [1] (d_e=0.7104/µ′_s) that is accounted as the distance from the interface at which the asymptotic component of the density of energy extrapolates to zero. Although this fact does not represent a ground for the heuristic boundary condition employed, it anyway suggests that for a source deeply sinked inside the medium the EBC tends within few per cents toward the boundary conditions of the classical Milne problem of photon transport. It can be also noticed that the differences between these two solutions become negligible for a source far from the boundary and for the cases analyzed in this work. For instance, for z ₀=1/µ′_s the differences observed are within 2% and the accuracy of these two boundary conditions is thus similar. The time-resolved reflectance emerging from the medium is derived by Fick’s law [15]

Fick’s law introduces an additional approximation at early times and the effects of this approximation decrease to zero for t≫1/(µ′_sv).

 

Fig. 1. Schematic of the medium, symbols used and positions of real and image sources, r ⁺and, r ^-, respectively.

Download Full Size | PDF

We remark a feature of Eq.(4) that needs to be known in its use. This feature is the tail of the ballistic peak due to radiation which has undergone a single forward scattering event, already mentioned by Paasschens [12] for Eq.(1). This tail determines two logarithmic singularities in Eq.(4), one at $t_1=\sqrt{\left({\rho }^2+z_0^2\right)}⁄v$ and the other at $t_2=\sqrt{\left({\rho }^2+{\left(z_0+2z_e\right)}^2\right)}⁄v$. Moreover, at t ₁ and t ₂ we have the singularities due to the ballistic peak of the two sources. These effects arise from the delta source considered in the solution and are confined to very early times (e.g. t ₁=4.7 ps and t ₂=8.7 ps for ρ=1 mm, µ′_s=1 mm^-1 and v=0.3 mm/ps). We want to stress that Eq. (4) is zero for t <t ₁ and, with the exception of the singularities at t ₁ and t ₂, it is a continuous function.

3. Results

The proposed solution has been obtained with a heuristic hypothesis that was validated by comparing the results of Eq. (4) with the results of MC simulations, here used as a gold standard. Details on the MC code used can be found in Ref. [16]. The correctness of the routines used to determine the scattering points was checked comparing some statistical parameters with analytical formulae [17]. The scattering functions used for the MC simulations were based on the Henyey-Greenstein (HG) model.

The comparison between the results of MC and RTE has been done for the time-resolved reflectance in a semi-infinite medium with µ′_s=1 mm^-1. Since both Eq. (4) (denoted RTE in the figures) and the MC use the same dependence on µ_a the comparison is done for a non absorbing medium without any loss of generality. The results obtained with the DE theory (Eq. (36) of Ref. [13] where only the first term of the series is retained) are also shown. The simulations have been firstly carried out for an isotropic source and for isotropic scattering as it is assumed in the derivation of Eq. (1). The position of the isotropic point-like source is considered at z ₀=1/µ′_s and v is assumed to be 0.3 mm/ps in the medium. In Fig. 2 (a) and (b) the reflectance from a semi-infinite medium is shown at ρ=1 mm and 5 mm, respectively. In Fig. 2 (c) and (d) the plots show the ratio of the MC results to those of RTE and DE. In Fig. 2 (a) and (b) the error bars on the MC data are omitted since they have similar size as the symbols. In Fig. 2 (c) and (d) the error bars are obtained from the statistical standard deviations on the MC data. Similar results have been obtained for ρ in the range 0.1–20 mm and for z ₀>0.2 mm. Both for RTE and for MC response the contribution of ballistic photons has not been included.

 

Fig. 2. Comparison between the results of MC, RTE and DE for the time-resolved reflectance from a non absorbing semi-infinite medium with µ′_s=1 mm^-1, ρ=1 and 5 mm, where an isotropic point-like source at depth z ₀=1/µ′_s and isotropic scattering is assumed. Figures (a) and (b) show the reflectance and figures (c) and (d) show the ratio of the MC results to those of RTE and DE.

Download Full Size | PDF

Figures 1(c) and 1(d) show that, except for the very first temporal windows (4.7≤t≤12 ps for ρ=1 mm or 16.9≤t≤26 ps for ρ=5 mm), the RTE solution is within 2% in agreement with the results of the MC simulations. Approximately, for ρ=0.1, 1, 5, 10 and 20 mm the error is less than 2% just 8 ps after t ₁ (that corresponds to 2.4 mean free paths). This accuracy is similar to that of Eq. (1). Thus even at very short distances the accuracy of Eq. (4) is similar to that of Eq. (1), assuring that the introduction of the heuristic boundary condition does not diminish the accuracy of the Paasschens’ solution. These results also show for the DE an accuracy significantly worse than the proposed heuristic solution. Moreover the solution from the DE is nonzero also for t<t ₁, due to the non-causality hypothesis assumed in this approach. The good agreement obtained at early times suggests that the CW reflectance obtained from an integration of Eq. (4) works well even at short distances and for high values of µ_a.

The validation of the presented heuristic solution has been done for the case of an isotropic source and of isotropic scattering. However, when light propagates through biological tissues or real media, this exact situation is almost never found. In order to test the applicability of Eq. (4) in more realistic conditions, with the MC code we simulated the case of a pencil light beam impiging a semi-infinite medium where a scattering function with anisotropy factor g=0.9, typical of biological tissue, was assumed. In Fig. 3 (a) and (b) the reflectance is shown at ρ=1 mm and 5 mm respectively. The ratio of the MC results to those of RTE and DE is plotted in figures (c) and (d). Equation (4) and the DE are evaluated with the source at z ₀=1/µ′_s. The distance ρ=1 mm could be an actual value for measurements at null source-detector distance. These results show an agreement between the heuristic solution and the MC results similar to that of Fig. 2 and we still can see clear advantages in using this solution instead of the DE. These results support the use of Eq. (4) for analyzing measurements on tissue at short source-receiver distances and at early times.

 

Fig. 3. Comparison between the results of MC, RTE and DE for the time-resolved reflectance from a non absorbing semi-infinite medium with µ′_s=1 mm^-1, ρ=1 and 5 mm. The MC simulations assume a pencil beam impinging the medium at z=0 and anisotropic scattering with g=0.9. Equation (4) for the RTE and Eq. (36) of Ref. [13] for the DE have been used with z ₀=1/µ′_s. Figures (a) and (b) show the reflectance and figures (c) and (d) show the ratio of the MC results to those of RTE and DE.

Download Full Size | PDF

4. Discussion and conclusions

The theory we have presented has been validated for refractive index matched at the boundary with the external medium. This situation has been chosen because it significantly simplifies the boundary condition and thus a better accuracy of the solution is expected. All the excellent agreement observed in the comparisons we showed refer to matched refractive index at the boundary. The agreement becomes worse for mismatched refractive index. However, we point out that for many applications the use of a matching liquid makes possible to carry out measurements with the refractive index matched at the boundary.

One of the main drawbacks of using RTE solutions is the long computation time that is usually required when numerical methods are used to solve integro-differential equations. The RTE solution proposed in this paper has the advantage to be an analytical closed-form solution with a negligible computation time. This fact is particularly important for inverse problems where Eq. (4) can be used to retrieve the optical properties of the medium in short time with a procedure that does not require to handle complex mathematics.

In conclusion, we have presented for the first time a heuristic analytical solution of the time dependent RTE for a semi-infinite medium based on an almost exact solution for the infinite medium where isotropic scattering is assumed. The comparisons with MC simulations for a point-like isotropic source and isotropic scattering (Fig. 2) have shown excellent performances: Except for the first temporal windows (i.e. approximately t ₁<t≤t ₁+8 ps for v=0.3 mm/ps and µ′_s=1 mm^-1) the proposed solution shows an error less than 2% even in regions where the DE shows errors larger than 10%. Also the comparisons with MC results for a medium with anisotropic scattering (g=0.9) illuminated by a pencil light beam (Fig. 3) showed an excellent agreement. This is surprising, since the proposed solution is based on the response for the infinite medium that is exact only for isotropic scatterers. These surprising good results are probably due to a compensation effect between the different approximations introduced to obtain the heuristic solution (boundary conditions, pencil beam modelled with a point-like isotropic source, isotropic scattering). Anyway, the excellent agreement shown by the comparisons in Fig. 3 for the temporal response even at short distances and at early times, suggests that also the CW reflectance obtained integrating Eq. (4) gives accurate results even at short distances and for high values of absorption for diffusive media with anisotropic scattering.

Thus, the use of Eq. (4) can significantly improve the analysis of measurements from tissue at short source-receiver distances and at early times with respect to solutions based on the DE. The use of this new formula could improve the retrieval of the optical properties from CW measurements at short distances since the DE has a poor accuracy in modelling light close to the source. These kind of measurements can be particularly important for small volume samples of biological tissues where the use of short interfiber distances is mandatory.We also stress that, following a similar procedure, it is possible to derive a similar expression to Eq. (4) for the slab geometry. Finally, we point out that the perturbation theory can also be implemented for the solution here described following the lines of Ref. [18].

The importance of the proposed solution is not only restricted to the biomedical optics field but it is in general addressed to a wide variety of physical phenomena where transport processes are involved.