## Abstract

The Green’s function of the time dependent radiative transfer equation for the semi-infinite medium is derived for the first time by a heuristic approach based on the extrapolated boundary condition and on an almost exact solution for the infinite medium. Monte Carlo simulations performed both in the simple case of isotropic scattering and of an isotropic point-like source, and in the more realistic case of anisotropic scattering and pencil beam source, are used to validate the heuristic Green’s function. Except for the very early times, the proposed solution has an excellent accuracy (>98% for the isotropic case, and >97% for the anisotropic case) significantly better than the diffusion equation. The use of this solution could be extremely useful in the biomedical optics field where it can be directly employed in conditions where the use of the diffusion equation is limited, e.g. small volume samples, high absorption and/or low scattering media, short source-receiver distances and early times. Also it represents a first step to derive tools for other geometries (e.g. slab and slab with inhomogeneities inside) of practical interest for noninvasive spectroscopy and diffuse optical imaging. Moreover the proposed solution can be useful to several research fields where the study of a transport process is fundamental.

©2007 Optical Society of America

## 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 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

_{s}*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

*µ*. Since

_{s}*µ′*=

_{s}*µ*(1-

_{s}*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(-

*µ*).

_{a}vtWhen 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*
_{0} 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, Φ

*: one positive from the source at*

_{i}**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*

_{0},

*t*), with refractive index matched at the boundary

According to Ref. [13], we assume an extrapolated distance *z _{e}*=2/3

*µ′*. Note that it differs in about 6%from the

_{s}*extrapolated end point*,

*d*, of the Milne problem [1] (

_{e}*d*=0.7104/

_{e}*µ′*) 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

_{s}*z*

_{0}=1/

*µ′*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]

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

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)}\u2044v$ and the other at ${t}_{2}=\sqrt{\left({\rho}^{2}+{\left({z}_{0}+2{z}_{e}\right)}^{2}\right)}\u2044v$. Moreover, at *t*
_{1} and *t*
_{2} 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*
_{1}=4.7 ps and *t*
_{2}=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*

_{1}and, with the exception of the singularities at

*t*

_{1}and

*t*

_{2}, 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

*µ*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

_{a}*z*

_{0}=1/

*µ′*and

_{s}*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}>0.2 mm. Both for RTE and for MC response the contribution of ballistic photons has not been included.

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*
_{1} (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*
_{1}, 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*
_{0}=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.

## 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*
_{1}<*t*≤*t*
_{1}+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.

## Acknowledgements

This work was supported by MIUR under the project PRIN2005 (prot. 2005025333).

## References and links

**1. **J. J. Duderstadt and W. R. Martin, *Transport Theory* (John Wiley&Sons, New York, 1979).

**2. **S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. **15**, R41–R93 (1999). [CrossRef]

**3. **A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. **50**, R1–R43 (2000). [CrossRef]

**4. **T. Feng, P. Edström, and M. Gulliksson, “Levenberg-Marquardt methods for parameter estimation problems in the radiative transfer equation,” Inverse Probl. **23**, 879–891 (2007). [CrossRef]

**5. **L. T. Perelman, J. Wu, Y. Wang, I. Itzkan, R. R. Dasari, and M. S. Feld, “Time-dependent photon migration using path integrals,” Phys. Rev. E **51**, 6134–6141 (1995). [CrossRef]

**6. **F. Martelli, M. Bassani, L. Alianelli, L. Zangheri, and G. Zaccanti, “Accuracy of the diffusion equation to describe photon migration through an infinite medium: Numerical and experimental investigation,” Phys. Med. Biol. **45**, 1359–1373 (2000). [CrossRef] [PubMed]

**7. **
For recent results: Special issue on recent development in biomedical optics, Phys. Med. Biol.49, N. 7 (2004). [PubMed]

**8. **A. Torricelli, A. Pifferi, L. Spinelli, R. Cubeddu, F. Martelli, S. Del Bianco, and G. Zaccanti, “Time-resolved reflectance at null source-detector separation: Improving contrast and resolution in diffuse optical imaging,” Phys. Rev. Lett. **95**, 078101 (2005). [CrossRef] [PubMed]

**9. **A. H. Hielscher, R. E. Alcouffe, and R. L. Barbour, “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys. Med. Biol. **43**, 1285–1302 (1998). [CrossRef] [PubMed]

**10. **H. Dehghani, S. R. Arridge, and M. Schweiger, “Optical tomography in the presence of void regions,” J. Opt. Soc. Am. A **17**,1659–1670 (2000). [CrossRef]

**11. **V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol. **23**, 313–320 (2005). [CrossRef] [PubMed]

**12. **J. C. J. Paasschens, “Solution of the time-dependent Boltzmann equation,” Phys. Rev. E **56**, 1135–1141 (1997). [CrossRef]

**13. **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**, 4587–4599 (1997). [CrossRef] [PubMed]

**14. **E. Zauderer, *Partial Differential Equations of Applied Mathematics*, (John Wiley&Sons, New York, 1989) Sec. 7.5, p. 484.

**15. **M. H. Lee, “Fick’s Law, Green-Kubo Formula, and Heisenberg’s Equation of Motion,” Phys. Rev. Lett. **85**, 2422–2425 (2000). [CrossRef] [PubMed]

**16. **F. Martelli, D. Contini, A. Taddeucci, and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation. II. Comparison with Monte Carlo results,” Appl. Opt. **30**, 4600–4612 (1997). [CrossRef]

**17. **G. Zaccanti, E. Battistelli, P. Bruscaglioni, and Q. N. Wei, “Analytic relationships for the statistical moments of scattering point coordinates for photon migration in a scattering medium,” Pure Appl. Opt. **3**, 897–905 (1994). [CrossRef]

**18. **F. Martelli, S. Del Bianco, and G. Zaccanti, “Perturbation model for light propagation through diffusive layered media,” Phys. Med. Biol. **50**, 2159–2166 (2005). [CrossRef] [PubMed]