Light diffusion in a turbid cylinder. I. Homogeneous case

doi:10.1364/OE.18.009456

« Show journal navigation

Light diffusion in a turbid cylinder. I. Homogeneous case

André Liemert and Alwin Kienle »View Author Affiliations

Optics Express, Vol. 18, Issue 9, pp. 9456-9473 (2010)
http://dx.doi.org/10.1364/OE.18.009456

View Full Text Article

Acrobat PDF (1101 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Article Navigation

▸ Article Outline
– Abstract
1. Introduction
2. Theory
3. Results
4. Discussion
– References and links

Article Tools
Full Text
Article Info
References
Cited By
Figures
Metrics
Download PDF

Viewing Equations in the
HTML Article

Optimal Viewing Recommendations

Full Text
Article Info
References (18)
Cited By
Figures (8)
Metrics

Abstract

This paper is the first of two dealing with light diffusion in a turbid cylinder. The diffusion equation was solved for a homogeneous finite cylinder that is illuminated at an arbitrary location. Three solutions were derived for an incident δ-light source in the steady-state, frequency, and time domains, respectively, applying different integral transformations. The performance of these solutions was compared with respect to accuracy and speed. Excellent agreement between the solutions, of which some are very fast (< 10ms), was found. Six of the nine solutions were extended to a circular flat beam which is incident onto the top side. Furthermore, the validity of the solutions was tested against Monte Carlo simulations.

1. Introduction

Light propagation in turbid media, e.g. biological tissue, is described with different equations depending on the considered applications and the required accuracy. On the microscopic scale Maxwell equations are usually used [1

F. Voit, J. Schäfer, and A. Kienle, “Light Scattering by Multiple Spheres: Comparison between Maxwell Theory and Radiative-Transfer-Theory Calculations,” Opt. Lett. 34, 2593–2595 (2009). [CrossRef] [PubMed]

], whereas on the mesoscopic scale transport equation and on the macroscopic scale diffusion equation are normally applied [2

A. Ishimaru, Wave Propagation and Scattering in Random Media , (Academic Press, New York, 1978).

The advantage of using the diffusion equation is that analytical solutions can be obtained for a variety of geometries. The disadvantage is that due to its approximations its applicability is restricted and, therefore, has to be checked against more precise theories, e.g. transport theory. Analytical solutions of the diffusion equation were found, e.g. for infinite, semi-infinite, and slab geometries [3–5

M. S. Patterson, B. Chance, and B. C. Wilson, “Time Resolved Reflectance and Transmittance for the Noninvasive Measurement of Tissue Optical Properties,” Appl. Opt. 28, 2331–2336 (1989). [CrossRef] [PubMed]

], further for a sphere [6

S. R. Arridge, M. Cope, and D. T. Delpy, “The Theoretical Basis for the Determination of Optical Pathlength in Tissue: Temporal and Frequency Analysis,” Phys. Med. Biol. 37, 1531–1560 (1992). [CrossRef] [PubMed]

, 7

B.W. Pogue and M.S. Patterson, “Frequency Domain Optical Absorption Spectroscopy of Finite Tissue Volumes Using Diffusion Theory,” Phys. Med. Biol. 39, 1157–1180 (1994). [CrossRef] [PubMed]

], for a parallel-epiped [8

A. Kienle, “Light Diffusion Through a Turbid Parallelepiped,” J. Opt. Soc. Am. A 22, 1883–1888 (2005). [CrossRef]

], and for a cylinder [6

, 7

B.W. Pogue and M.S. Patterson, “Frequency Domain Optical Absorption Spectroscopy of Finite Tissue Volumes Using Diffusion Theory,” Phys. Med. Biol. 39, 1157–1180 (1994). [CrossRef] [PubMed]

, 9

A. Sassaroli, F. Martelli, D. Imai, and Y. Yamada, “Study on the Propagation of Ultra-Short Pulse Light in Cylindrical Optical Phantoms,” Phys. Med. Biol. 44, 2747–2763 (1999). [CrossRef] [PubMed]

Besides these homogeneous geometries solutions of the diffusion equation were also derived for inhomogeneous geometries. Especially, layered geometries aroused interest because many different parts of the body exhibit layered structures. In this study we present solutions of the diffusion equation for a homogeneous cylinder, whereas in the accompanying paper solutions for a cylinder having an arbitrary number of layers are given [10

A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. II. Layered case,” accepted (2010).

In the field of biomedical optics the applicability of the derived solutions is twofold. First, in phantom experiments often cylindrical geometries are used, so that it is obviously required to apply the corresponding solutions. Second, for in-vivo experiments the cylindrical geometry can be used to describe approximately parts of the human body or of animals, e.g. the arm, the leg, or the finger. Of course, these examples are not exactly described by a cylinder but, definitely, the assumption of a cylindrical geometry is more exact than the mostly used semi-infinite geometry. We note that the derived equations can be used to obtain the corresponding solutions of the heat conduction equation, as the diffusion equation includes the heat equation as a special case. In this article we present three different solutions of the diffusion equation for a homogeneous finite cylinder in the steady-state, frequency, and time domains, respectively. For these nine solutions the assumed pencil beam is incident onto the cylinder barrel or onto the cylinder top (or bottom). The solutions were compared among each other for the different domains. We discuss how the solutions given in the literature so far are related to ours. In addition, we derived the solutions of the diffusion equations for a circular flat beam with a finite diameter incident onto the cylinder top (or bottom). Finally, the solutions were compared to results obtained by the Monte Carlo method, which is a numerical solution of the transport equation. In general, the derived solutions can be easily implemented using standard programming languages. In addition, we provide a freely available executable code on the internet [11

http://www.uni-ulm.de/ilm/index.php?id=10020200.

2. Theory

2.1. Diffusion Theory

In this section we present different solutions of the diffusion equation for a finite homogeneous cylinder. In the first part the finite Hankel and cosine transforms are used to obtain an ordinary differential equation for the z-coordinate, see Fig. 1. The boundary conditions in z-direction are then implemented applying two different approaches in the frequency domain: use of an infinite series of point sources (version A) and use of hyperbolic sine functions (version B), and applying two approaches in the time domain: use of the finite sine transform (version A) and again use of an infinite series of point sources (version B). In addition, for these six cases we give also the solutions for a flat beam incident onto the cylinder top. In the second part we use the finite sine transform in z-direction and then solve the resulting modified Bessel differential equation both in the frequency and time domains. We note that the results in the steady-state domain can be easily obtained from those derived in the frequency domain by setting the modulation frequency ω = 0.

In literature Arridge et al. presented a two-dimensional solution of the diffusion equation for a cylinder in the frequency domain solving the modified Bessel differential equation [6

]. In addition, they gave the solution of the diffusion equation for a finite cylinder in the frequency and time domains, but only for the incident plane. Sassaroli et al. extended this solution to locations outside the incident plane and for the extrapolated boundary condition [9

]. Pogue and Patterson indicated the solution of the diffusion equation for a finite cylinder in the time domain, which corresponds to version B in this study, but they did not present a derivation of the solution [7

B.W. Pogue and M.S. Patterson, “Frequency Domain Optical Absorption Spectroscopy of Finite Tissue Volumes Using Diffusion Theory,” Phys. Med. Biol. 39, 1157–1180 (1994). [CrossRef] [PubMed]

Figure 1 shows the cylindrical geometry used in the simulations. The center of the coordinate system is the middle of the cylinder top. The radius and height of the cylinder are denoted as a and l_z , respectively. The incident light beam position can be chosen arbitrarily, e.g. it can be incident onto the top (or bottom) and onto the barrel. It is assumed that the incident light beam can be represented by a point source in the turbid medium in direction of the beam at a distance of l/(μ′_s + μ_a ) from the location of (perpendicular) incidence at the cylinder barrel, see Fig. 1, where μ_a and μ′_s are the absorption and reduced scattering coefficients of the homogeneous cylinder. The position of the point source is given in cylindrical coordinates r⃗₀ = (ρ ₀, ϕ ₀,z ₀). In case of an oblique incident beam the Snel’s law has to be applied to determine the position of the point source. Extrapolated boundary conditions are used to describe the influence of the border between the turbid medium and the surrounding, characterized by the refractive indices, n ₁ and n ₀, respectively. As in earlier papers we use the fluence rate and the flux term for calculation of the reflectance in the steady-state domain [12

A. Kienle and M. S. Patterson, “Improved Solutions of the Steady-State and the Time-Resolved Diffusion Equations for Reflectance from a Semi-Infinite Turbid Medium,” J. Opt. Soc. Am. A 14, 246–254 (1997). [CrossRef]

], whereas for the other solutions only the flux term (Fick’s law) is applied. The rationale for this is the better agreement to Monte Carlo simulations when using the fluence and flux terms for calculation of the spatially resolved reflectance [12

2.1.1. Solutions derived via the finite Hankel transform

We start with the diffusion equation in the frequency domain

ΔΦ (\vec{r}, ω) - (\frac{μ_{a}}{D} + i \frac{ω}{Dc}) Φ (\vec{r}, ω) = - \frac{1}{D} δ (\vec{r} - {\vec{r}}_{0}),

(1)

Fig. 1. Scheme of the cylinder used in the calculations. Exemplarily, the beam is incident perpendicular onto the cylinder top at r⃗ = (ρ ₀, ϕ ₀,0). The isotropic point source is located at r⃗ = (ρ ₀, ϕ, z ₀)

where Φ, D = 1/(3μ′_s), c, and ω denote the fluence rate, the diffusion coefficient, the speed of light in the turbid medium and the angular frequency of the intensity modulated light, respectively. Equation (1) can be rewritten by using cylindrical coordinates as follows

\frac{\partial^{2}}{\partial ρ^{2}} Φ + \frac{1}{ρ} \frac{\partial}{\partial ρ} Φ + \frac{1}{ρ^{2}} \frac{\partial^{2}}{\partial ϕ^{2}} Φ + \frac{\partial^{2}}{\partial z^{2}} Φ - (\frac{μ_{a}}{D} + i \frac{ω}{Dc}) Φ = - \frac{1}{Dρ} δ (ρ - ρ_{0}) δ (ϕ - ϕ_{0}) δ (z - z_{0}),

(2)

where we use the abbreviation Φ = Φ(r⃗, ω). First, a cosine transform of the form

Φ (ϕ, m) = \int_{0}^{2 π} Φ (ϕ') \cos (m (ϕ - ϕ')) dϕ', m = 0,1,2, \dots

(3)

is applied to Eq. (2) in order to eliminate the ϕ derivative. Thus, we obtain

\frac{\partial^{2}}{\partial ρ^{2}} Φ + \frac{1}{ρ} \frac{\partial}{\partial ρ} Φ - \frac{m^{2}}{ρ^{2}} Φ + \frac{\partial^{2}}{\partial z^{2}} Φ - (\frac{μ_{a}}{D} + i \frac{ω}{Dc}) Φ = - \frac{\cos (m (ϕ - ϕ_{0}))}{Dρ} δ (ρ - ρ_{0}) δ (z - z_{0}) .

(4)

Equation (4) can be reduced to an ordinary differential equation by using the finite Hankel transform of m^th order

Φ (s_{n}) = \int_{0}^{a'} ρ Φ (ρ) J_{m} (s_{n} ρ) d ρ,

(5)

where s_n are the positive roots of the Bessel functions of first kind and m^th order divided by a′, i.e. J_m (a′s_n ) = 0. The roots of the Bessel functions are determined using Halley’s method and applying the start values described by Abramowitz and Stegun [17

M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions , (Dover Publications, New York, 1971).

]. The upper limit of the integral a′ results from the extrapolated boundary condition Φ(ρ = a′) = 0, where we use a′ = a + z_b . The extrapolation length z_b is calculated with

z_{b} = \frac{1 + R_{eff}}{1 - R_{eff}} 2 D,

(6)

where R_eff represents the fraction of photons that is internally diffusely reflected at the boundary to the non-scattering surrounding. R_eff was calculated by solving the integrals presented by Haskell et al [13

R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. McAdams, and B. J. Tromberg, “Boundary Conditions for the Diffusion Equation in Radiative Transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). [CrossRef]

]. For the finite Hankel transform the relation

\int_{0}^{a'} ρ (\frac{\partial^{2}}{\partial ρ^{2}} Φ + \frac{1}{ρ} \frac{\partial}{\partial ρ} Φ - \frac{m^{2}}{ρ^{2}} Φ) J_{m} (s_{n} ρ) d ρ = - s_{n}^{2} Φ (s_{n}) + a' s_{n} J_{m + 1} (a' s_{n}) Φ (ρ = a')

(7)

is valid [14

E. Meissel and G.B. Mathews, A Treatise on Bessel Functions and Their Applications to Physics , (Bibliobazaar, 2008).

]. By applying the finite Hankel transform to Eq. (4) and by considering the boundary condition at ρ = a′, we obtain the following ordinary differential equation where Φ = Φ(s_n , ϕ, m, z, ω)

\frac{\partial^{2}}{\partial z^{2}} Φ - (s_{n}^{2} + \frac{μ_{a}}{D} + i \frac{ω}{Dc}) Φ = - \frac{1}{D} J_{m} (s_{n} ρ_{0}) \cos (m (ϕ - ϕ_{0})) δ (z - z_{0}) .

(8)

We define G(s_n , z, ω as the solution of

\frac{\partial^{2}}{\partial z^{2}} G (s_{n}, z, ω) - α^{2} G (s_{n}, z, ω) = - \frac{1}{D} δ (z - z_{0}),

(9)

where α ² = s ² _n + μ_a /D + iω/(Dc). This equation is similar to the equation we derived solving the diffusion equation with lateral infinitely large extensions [15

A. Kienle, M. S. Patterson, N. Dögnitz, R. Bays, G. Wagnières, and H. van den Bergh, “Nonin-vasive Determination of the Optical Properties of Two-Layered Turbid Media,” Appl. Opt. 37, 779–791 (1998). [CrossRef]

]. In addition, we note that for the case of an N-layered cylinder the corresponding N-layered equation has to be solved, see accompanying paper [10

A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. II. Layered case,” accepted (2010).

]. Equation (8) is solved by using the solution of Eq. (9) given below in version A and B. Thus, Φ can be obtained with

Φ (s_{n}, ϕ, m, z, ω) = G (s_{n}, z, ω) J_{m} (s_{n} ρ_{0}) \cos (m (ϕ - ϕ_{0})) .

(10)

We note that Eqs. 3 and 5 together give the following transform

Φ (s_{n}, ϕ, m) = \int_{0}^{2 π} \int_{0}^{a'} ρ Φ (ρ, ϕ') J_{m} (s_{n} ρ) \cos (m (ϕ - ϕ')) d ρ d ϕ',

(11)

with which Eq. (8) was obtained from Eq. (2). Next, an inversion formula of Eq. (11) is obtained by, first, deriving the inversion formula of Eq. (3) and then of Eq. (5). By summing both sides of Eq. (3) over all m ≥ 0, the following equation is obtained

\sum_{m = 0}^{\infty} Φ (ϕ, m) = \int_{0}^{2 π} Φ (ϕ') \sum_{m = 0}^{\infty} \cos (m (ϕ - ϕ')) d ϕ' .

(12)

By applying the relation

\sum_{m = 0}^{\infty} \cos (m (ϕ - ϕ')) = \frac{1}{2} + πδ (ϕ - ϕ')

(13)

we get

\int_{0}^{2 π} Φ (ϕ') [\frac{1}{2} + πδ (ϕ - ϕ')] dϕ' = Φ (ϕ, 0) + π Φ (ϕ)

(14)

and, finally, we obtain

Φ (ϕ) = \frac{1}{2 π} \sum_{m = 0}^{\infty} (2 - δ_{m, 0}) Φ (ϕ, m) .

(15)

An inverse transform for the second transform, the Hankel transform (Eq. (5)), is found by using the completeness relation for the Bessel functions [14

E. Meissel and G.B. Mathews, A Treatise on Bessel Functions and Their Applications to Physics , (Bibliobazaar, 2008).

]

\sum_{n = 1}^{\infty} \frac{ρ J_{m} (s_{n} ρ) J_{m} (s_{n} ρ_{0})}{J_{m + 1}^{2} (a' s_{n})} = \frac{{a'}^{2}}{2} δ (ρ - ρ_{0}) .

(16)

By applying Eq. (15) and Eq. (16) to Eq. (11), we find the inversion formula as follows

Φ (ρ, ϕ) = \frac{1}{π {a'}^{2}} \sum_{m = - \infty}^{\infty} (2 - δ_{m, 0}) \sum_{n = 1}^{\infty} Φ (s_{n}, ϕ, m) J_{m} (s_{n} ρ) J_{m + 1}^{- 2} (a' s_{n}) .

(17)

By using Eq. (10), the solution of Eq. (1) in the frequency domain is given by

Φ (\vec{r}, ω) = \frac{1}{π {a'}^{2}} \sum_{m = 0}^{\infty} (2 - δ_{m, 0}) \cos (mφ) \sum_{n = 1}^{\infty} G (s_{n}, z, ω) J_{m} (s_{n} ρ_{0}) J_{m} (s_{n} ρ) J_{m + 1}^{- 2} (a' s_{n}),

(18)

where φ = ϕ - ϕ ₀. In order to simplify Eq. (18) the following recurrence relation is applied

J_{m + 1} (a' s_{n}) = \frac{2 m}{a' s_{n}} J_{m} (a' s_{n}) - J_{m - 1} (a' s_{n}) = - J_{m - 1} (a' s_{n}) .

(19)

Thus, we obtain for the general case of a pencil beam incident onto an arbitrary location of the cylinder

Φ (\vec{r}, ω) = \frac{1}{π {a'}^{2}} \sum_{m = - \infty}^{\infty} \cos (mφ) \sum_{n = 1}^{\infty} G (s_{n}, z, ω) J_{m} (s_{n} ρ_{0}) J_{m} (s_{n} ρ) J_{m + 1}^{- 2} (a' s_{n}) .

(20)

For the special case where the pencil beam is incident onto the center of the cylinder top, r⃗ = (0,0,0), Eq. (20) becomes

Φ (\vec{r}, ω) = \frac{1}{π {a'}^{2}} \sum_{n = 1}^{\infty} G (s_{n}, z, ω) J_{0} (s_{n} ρ) J_{1}^{- 2} (a' s_{n}) .

(21)

As stated above, the light transmitted from the cylinder is calculated using Fick’s law. The results for the reflectance and transmittance from the top and bottom of the cylinder are easily obtained and, thus, are not given. For the case of the transmittance from the cylinder barrel the following equation

T (ϕ, z, ω) = - D \partial Φ (ρ, ϕ, z, ω) / \partial ρ |_{ρ = a}

(22)

is used to obtain

T (ϕ, z, ω) = \frac{D}{π {a'}^{2}} \sum_{m = - \infty}^{\infty} \cos (mφ)

\times \sum_{n = 1}^{\infty} G (s_{n}, z, ω) J_{m} (s_{n} ρ_{0}) (s_{n} J_{m + 1} (a s_{n}) - \frac{m}{a} J_{m} (a s_{n})) J_{m + 2}^{- 2} (a' s_{n}) .

(23)

For the special case of central incidence we get

T (z, ω) = \frac{D}{π {a'}^{2}} \sum_{n = 1}^{\infty} s_{n} G (s_{n}, z, ω) J_{1} (a s_{n}) J_{1}^{- 2} (a' s_{n}) .

(24)

The results presented above were obtained by assuming an incident δ-source. For the derivation of the solution of a finite flat beam incident onto the center of the cylinder top we use the diffusion equation for the rotationally symmetric case (no ϕ-dependence)

\frac{\partial^{2}}{\partial ρ^{2}} Φ + \frac{1}{ρ} \frac{\partial}{\partial ρ} Φ + \frac{\partial^{2}}{\partial z^{2}} Φ - (\frac{μ_{a}}{D} + i \frac{ω}{Dc}) Φ = - \frac{1}{D} S (\vec{r}, ω),

(25)

where S(r⃗, ω) is a rotationally symmetric source. For a flat beam we have

S (\vec{r}, ω) = \frac{1}{π ρ_{w}^{2}} circ (\frac{ρ}{ρ_{w}}) δ (z - z_{0}),

(26)

where ρ_w is the radius of the flat beam and circ is defined as

circ (\frac{ρ}{ρ_{w}}) = {\begin{matrix} 1, ρ \leq ρ_{w} \\ 0, ρ > ρ_{w} \end{matrix} .

(27)

By employing the finite Hankel transform of zero order, Eq. (25) becomes

\frac{\partial^{2}}{\partial z^{2}} Φ - (\frac{μ_{a}}{D} + s_{n}^{2} + i \frac{ω}{Dc}) Φ = - \frac{1}{D} S (s_{n}, z, ω) = - \frac{1}{D} \frac{J_{1} (s_{n} ρ_{w})}{π s_{n} ρ_{w}} δ (z - z_{0}) .

(28)

By using the solution of Eq. (9), G(s_n ,z,ω), Eq. (28) is solved with

Φ (s_{n}, z, ω) = \frac{J_{1} (s_{n} ρ_{w})}{π s_{n} ρ_{w}} G (s_{n}, z, ω) .

(29)

The inverse relation for the finite Hankel transform of zero order is

Φ (\vec{r}, ω) = \frac{2}{{a'}^{2}} \sum_{n = 1}^{\infty} Φ (s_{n}, z, ω) J_{0} (s_{n} ρ) J_{1}^{- 2} (s_{n} a') .

(30)

Thus, the solution of the diffusion equation for an incident flat beam is given by

Φ (\vec{r}, ω) = \frac{2}{{π a'}^{2} ρ_{w}} \sum_{n = 1}^{\infty} G (s_{n}, z, ω) \frac{J_{0} (s_{n} ρ) J_{1} (s_{n} ρ_{w})}{s_{n} J_{1}^{2} (s_{n} a')} .

(31)

We note that Eq. (31) is also used for solving the diffusion equation for N-layered turbid media, see the accompanying paper [10

A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. II. Layered case,” accepted (2010).

]. The final step for obtaining the complete solutions is the derivation of the Green’s function, G, see Eq. (9). The boundary conditions in z-direction for the homogeneous cylinder are given by

Φ (ρ, ϕ, z = - z_{b}, ω) = 0

Φ (ρ, ϕ, z = l_{z} + z_{b}, ω) = 0 .

(32)

Below we present two solutions in the frequency and two solutions in the time domains.

Solutions in the frequency domain

Version A

We start with the ordinary differential equation given in Eq. (9). It can be solved with the method of image sources where the exponential functions are summed to satisfy the boundary conditions for the z-direction given in Eq. (32). The solution is

G (s_{n}, z, ω) = \frac{1}{2 Dα} \sum_{k = - \infty}^{\infty} \exp (- α ∣ z - z_{1 k} ∣) - \exp (- α ∣ z - z_{2 k} ∣),

(33)

where

z_{1 k} = 2 k l_{z} + 4 k z_{b} + z_{0}

z_{2 k} = 2 k l_{z} + (4 k - 2) z_{b} - z_{0}

α = \sqrt{\frac{μ_{a}}{D} + s_{n}^{2} + i \frac{ω}{Dc}} .

(34)

Thus, by inserting Eq. (33) in Eq. (20) we obtain the fluence rate for the pencil beam incident at an arbitrary location onto the cylinder as

Φ (\vec{r}, ω) = \frac{1}{2 πD {a'}^{2}} \sum_{m = - \infty}^{\infty} \cos (mφ) \sum_{n = 1}^{\infty} \frac{J_{m} (s_{n} ρ_{0}) J_{m} (s_{n} ρ)}{α J_{m + 1}^{2} (a' s_{n})}

\times \sum_{k = - \infty}^{\infty} \exp (- α ∣ z - z_{1 k} ∣) - \exp (- α ∣ z - z_{2 k} ∣) .

(35)

For the case when the pencil beam is incident at the center of the cylinder Eq. (35) is simplified to

Φ (\vec{r}, ω) = \frac{1}{2 πD {a'}^{2}} \sum_{n = 1}^{\infty} \frac{J_{0} (s_{n} ρ)}{α J_{1}^{2} (a' s_{n})} \sum_{k = - \infty}^{\infty} \exp (- α ∣ z - z_{1 k} ∣) - \exp (- α ∣ z - z_{2 k} ∣) .

(36)

Version B

There is another possibility to solve Eq. (9) for the boundaries given in Eq. (32). By applying the method described in reference [15

] the Green’s function for a one-dimensional layer is given in the interval 0 ≤ z < l_z as

G (s_{n}, z, ω) = \frac{\exp (- α ∣ z - z_{0} ∣) - \exp (- α (z + z_{0} + 2 z_{b}))}{2 Dα} - \frac{\sinh [z_{0} + z_{b}]}{Dα}

\times \frac{\sinh [α (z + z_{b})] \exp (- α (l_{z} + 2 z_{b}))}{\sinh [α (l_{z} + 2 z_{b})]} .

(37)

We note that for the evaluation of Eq. (37) no approximation is needed, in contrast to the solutions presented earlier [15

]. In general, the solution calculated with Eq. (37) is faster than the solution obtained with Eq. (33), if l_z is small.

Finally, Eq. (37) is inserted in Eq. (20) to obtain the fluence rate for a pencil beam incident at an arbitrary position, in Eq. (21) for a pencil beam incident in the middle of the cylinder top, or in equation Eq. (31) for an incident flat beam.

Solutions in the time domain

Version A

For the derivation of the solution of the fluence rate in the time domain we start with the ordinary differential equation [Eq. (8)]. The boundary conditions given in Eq. (32) are fulfilled by using a finite sine transform:

Φ (λ_{k}) = \int_{z_{1}}^{z_{2}} Φ (z) \sin (λ_{k} (z - z_{1})) dz, λ_{k} = kπ / (z_{2} - z_{1}), k = 1,2, \dots,

(38)

where z ₁ = -z_b and z ₂ = l_z +z_b . By applying Eq. (38) and the following relation

\int_{z_{1}}^{z_{2}} \frac{\partial^{2}}{\partial z^{2}} Φ (z) \sin (λ_{k} (z - z_{1})) d z = λ_{k} [Φ (z_{1}) - {(- 1)}^{k} Φ (z_{2})] - λ_{k}^{2} Φ (λ_{k})

(39)

we obtain from Eq. (8)

(s_{n}^{2} + λ_{k}^{2} + \frac{μ_{a}}{D} + i \frac{ω}{Dc}) Φ = \frac{1}{D} J_{m} (s_{n} ρ_{0}) \sin (λ_{k} (z_{0} + z_{b})) \cos (m (ϕ - ϕ_{0})) .

(40)

Thus, the fluence rate in the frequency domain is given by

Φ (s_{n}, ϕ, m, λ_{k}, ω) = \frac{c J_{m} (s_{n} ρ_{0}) \sin (λ_{k} (z_{0} + z_{b})) \cos (m (ϕ - ϕ_{0}))}{Dc (s_{n}^{2} + λ_{k}^{2}) + μ_{a} c - iω} .

(41)

The fluence rate in the time domain is obtained from the solutions in the frequency domain using the following Fourier transform

\int_{- \infty}^{\infty} \exp (- at) ε (t) \exp (- iωt) dt = \frac{1}{a + iω}, a > 0,

(42)

where ε(t) describes the Heaviside step function. By employing Eq. (42), Eq. (41) becomes

Φ (s_{n}, ϕ, m, λ_{k}, t) = c J_{m} (s_{n} ρ_{0}) \sin (λ_{k} (z_{0} + z_{b})) \cos (m (ϕ - ϕ_{0}))

\times \exp (- μ_{a} ct) \exp (- Dc λ_{k}^{2} t) \exp (- Dc s_{n}^{2} t), t > 0 .

(43)

We use the completeness relation for the sine function to obtain an inversion formula of Eq. (38)

Φ (z) = \frac{2}{z_{2} - z_{1}} \sum_{k = 1}^{\infty} Φ (λ_{k}) \sin (λ_{k} (z - z_{1})) .

(44)

Finally, the inversion formula given in Eq. (20) is used to obtain the time resolved fluence rate in real space

Φ (\vec{r}, t) = \frac{2}{π {a'}^{2} (l_{z} + 2 z_{b})} \exp (- μ_{a} ct) \sum_{k = 1}^{\infty} \exp (- Dc λ_{k}^{2} t) \sin (λ_{k} (z_{0} + z_{b})) \sin (λ_{k} (z + z_{b}))

\times \sum_{m = - \infty}^{\infty} \cos (mφ) \sum_{n = 1}^{\infty} \exp (- Dc s_{n}^{2} t) J_{m} (s_{n} ρ) J_{m} (s_{n} ρ) J_{m + 1}^{- 2} (a' s_{n}) .

(45)

For a pencil beam incident onto the center of the cylinder top Eq. (45) is simplified to

Φ (\vec{r}, t) = \frac{2 c}{π {a'}^{2} (l_{z} + 2 z_{b})} \exp (- μ_{a} ct) \sum_{k = 1}^{\infty} \exp (- Dc λ_{k}^{2} t) \sin (λ_{k} (z_{0} + z_{b})) \sin (λ_{k} (z + z_{b}))

\times \sum_{n = 1}^{\infty} \exp (- Dc s_{n}^{2} t) J_{0} (s_{n} ρ) J_{1}^{- 2} (a' s_{n}) .

(46)

Version B

We also present a second possibility to solve Eq. (8). Considering again the boundary conditions given in Eq. (32) the series of point sources already applied for version A of the frequency solutions is used, see Eq. (33), where α is now

α = \sqrt{μ_{a} / D + s_{n}^{2} + iω / Dc}

. By employing the Fourier transform

\int_{- \infty}^{\infty} \frac{1}{\sqrt{πt}} \exp (- at) \exp (\frac{b^{2}}{4 t}) ε (t) \exp (- iωt) d t = \frac{\exp (- ∣ b ∣ \sqrt{a + iω})}{\sqrt{a + iω}}, a > 0

(47)

the solution for the homogeneous cylinder in the time domain is for an arbitrary incident pencil beam

Φ (\vec{r}, t) = \frac{1}{2 π a^{' 2}} \sqrt{\frac{c}{Dπt}} \exp (- μ_{a} ct) \sum_{k = - \infty}^{\infty} \exp (- \frac{{(z - z_{1 k})}^{2}}{4 Dct}) - \exp (- \frac{{(z - z_{2 k})}^{2}}{4 Dct})

\times \sum_{m = - \infty}^{\infty} \cos (mφ) \sum_{n = 1}^{\infty} \exp (- Dc s_{n}^{2} t) J_{m} (s_{n} ρ_{0}) J_{m} (s_{n} ρ) J_{m + 1}^{- 2} (a' s_{n}) .

(48)

For the pencil beam incident at the middle of the cylinder top we get

Φ (\vec{r}, t) = \frac{1}{2 π a^{' 2}} \sqrt{\frac{c}{Dπt}} \exp (- μ_{a} ct) \sum_{k = - \infty}^{\infty} \exp (- \frac{{(z - z_{1 k})}^{2}}{4 Dct}) - \exp (- \frac{{(z - z_{2 k})}^{2}}{4 Dct})

\sum_{n = 1}^{\infty} \exp (- Dc s_{n}^{2} t) J_{0} (s_{n} ρ_{0}) J_{1}^{- 2} (a' s_{n}) .

(49)

2.1.2. Solutions derived via modified Bessel differential equation

In this section we derive the solutions of the homogeneous diffusion equation for an incident pencil beam by solving the modified Bessel differential equation. We start in the frequency domain and, then, give the time domain solution.

The modified Bessel differential equation is obtained by neglecting the z-dependency in Eq. (2)

\frac{\partial^{2}}{\partial ρ^{2}} Φ + \frac{1}{ρ} \frac{\partial}{\partial ρ} Φ + \frac{1}{ρ^{2}} \frac{\partial^{2}}{\partial ϕ^{2}} - α^{2} Φ = - \frac{1}{D} δ (ρ - ρ_{0}) δ (ϕ - ϕ_{0}),

(50)

where

α = \sqrt{μ_{a} / D + iω / Dc}

. This differential equation can be solved considering the radial boundary condition Φ(ρ = a′, ω) = 0. The solution can be derived by using a formula given in reference [6

] as

Φ (ρ, φ, ω) = \frac{1}{2 πD} [K_{0} (α \sqrt{ρ^{2} + ρ_{0}^{2} - 2 ρ ρ_{0} \cos φ}) .

- \sum_{m = - \infty}^{\infty} \cos (mφ) \frac{I_{m} (α ρ_{0}) K_{m} (αa') I_{m} (αρ)}{I_{m} (αa')}] .

(51)

I_m and K_m denote the modified Bessel functions of first and second kind of m^th order. By setting

α = \sqrt{μ_{a} / D + λ_{k}^{2} + iω / Dc}

and multiplying the right side of Eq. (50) with sin(λ_k (z ₀ + z_b )), we obtain the same result as if we had performed a finite sine transform related to the z-coordinate, compare Eq. (38). By applying the inverse transform, see Eq. (44), we get the fluence rate for an arbitrary incident beam position

Φ (\vec{r}, ω) = \frac{1}{4 πD} \sum_{k = - \infty}^{\infty} \frac{\exp (- μ_{eff} r_{1})}{r_{1}} - \frac{\exp (- μ_{eff} r_{2})}{r_{2}} - \frac{1}{πD (l_{z} + 2 z_{b})}

\times \sum_{k = 1}^{\infty} \sin (λ_{k} (z_{0} + z_{b})) \sin (λ_{k} (z + z_{b})) \sum_{m = - \infty}^{\infty} \cos (mφ) \frac{I_{m} (α ρ_{0}) K_{m} (αa') I_{m} (αρ)}{I_{m} (αa')},

(52)

where

r_{1} = \sqrt{ρ^{2} + ρ_{0}^{2} - 2 ρ ρ_{0} \cos φ + {(z - z_{1 k})}^{2}}

r_{2} = \sqrt{ρ^{2} + ρ_{0}^{2} - 2 ρ ρ_{0} \cos φ + {(z - z_{2 k})}^{2}}

μ_{eff} = \sqrt{μ_{a} / D + iω / Dc} .

(53)

The transmitted light from the cylinder barrel is found by using the ρ derivative for the modified Bessel functions, see Eq. (22). The result is

T (φ, z, ω) = \frac{a - ρ_{0} \cos φ}{4 π} \sum_{k = - \infty}^{\infty} (μ_{eff} + 1 / r_{1}) \frac{\exp (- μ_{eff} r_{1})}{r_{1}^{2}} - (μ_{eff} + 1 / r_{2}) \frac{\exp (- μ_{eff} r_{2})}{r_{2}^{2}}

+ \frac{1}{π (l_{z} + 2 z_{b})} \sum_{k = 1}^{\infty} \sin (λ_{k} (z_{0} + z_{b})) \sin (λ_{k} (z + z_{b}))

\times \sum_{m = - \infty}^{\infty} \cos (mφ) \frac{I_{m} (α ρ_{0}) K_{m} (αa') [α I_{m + 1} (αρ) + m / a I_{m} (αa)]}{I_{m} (αa')} .

(54)

For a pencil beam incident onto the middle of the cylinder top the fluence rate is

Φ (\vec{r}, ω) = \frac{1}{4 πD} \sum_{k = - \infty}^{\infty} \frac{\exp (- μ_{eff} r_{1})}{r_{1}} - \frac{\exp (- μ_{eff} r_{2})}{r_{2}}

- \frac{1}{πD (l_{z} + 2 z_{b})} \sum_{k = 1}^{\infty} \sin (λ_{k} (z_{0} + z_{b})) \sin (λ_{k} (z + z_{b})) \frac{K_{0} (αa') I_{0} (αρ)}{I_{0} (αa')},

(55)

and the transmittance is obtained with

T (z, ω) = \frac{a}{4 π} \sum_{k = - \infty}^{\infty} (μ_{eff} + 1 / r_{1}) \frac{\exp (- μ_{eff} r_{1})}{r_{1}^{2}} - (- μ_{eff} + 1 / r_{2}) \frac{\exp (- μ_{eff} r_{2})}{r_{2}^{2}}

- \frac{1}{π (l_{z} + 2 z_{b})} \sum_{k = 1}^{\infty} α \sin (λ_{k} (z_{0} + z_{b})) \sin (λ_{k} (z + z_{b})) \frac{K_{0} (αa') I_{1} (αa)}{I_{0} (αa')} .

(56)

For the solution in the time domain, the frequency domain solution, Φ(r⃗, ω), is used and split into real and imaginary parts [17

M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions , (Dover Publications, New York, 1971).

]. Then, the inverse Fourier transform is calculated by using the FFT algorithm.

Alternatively, we used the ‘damping’ rule of the Laplace transform for separation of that part in Eq. (52) which contains the z-coordinate [18

A. Liemert and A. Kienle, “Light Diffusion in N-layered Turbid Media: Frequency and Time Domains,” J. Biomed. Opt. 15, 025002 (2010). [CrossRef] [PubMed]

]. The result is

Φ (\vec{r}, t) = \frac{1}{4 πD \sqrt{πDct}} \sum_{k = - \infty}^{\infty} \exp (- \frac{{(z - z_{1 k})}^{2}}{4 Dct}) - \exp (- \frac{{(z - z_{2 k})}^{2}}{4 Dct})

\times [\frac{1}{2 t} \exp (- μ_{a} ct) \exp (- \frac{ρ^{2} - ρ_{0}^{2} - 2 ρ ρ_{0} \cos ρ}{4 Dct})

- \frac{1}{2 π} \sum_{m = - \infty}^{\infty} \cos (mφ) \int_{- \infty}^{\infty} \frac{I_{m} (μ_{eff} ρ_{0}) K_{m} (μ_{eff} a') I_{m} (μ_{eff} ρ)}{I_{m} (μ_{eff} a')} \exp (iωt) d ω],

(57)

and the transmittance is obtained with

T (φ, z, t) = \frac{1}{4 π \sqrt{πDct}} \sum_{k = - \infty}^{\infty} \exp (- \frac{{(z - z_{1 k})}^{2}}{4 Dct}) - \exp (- \frac{{(z - z_{2 k})}^{2}}{4 Dct})

\times [\frac{a - ρ_{0} \cos φ}{4 Dc t^{2}} \exp (- μ_{a} ct) \exp (- \frac{a^{2} - ρ_{0}^{2} - 2 a ρ_{0} \cos ρ}{4 Dct})

+ \frac{1}{2 π} \sum_{m = - \infty}^{\infty} \cos (mφ) \int_{- \infty}^{\infty} \frac{I_{m} (μ_{eff} ρ_{0}) K_{m} (μ_{eff} a') I_{m} (μ_{eff} a)}{I_{m} (μ_{eff} a')} \exp (iωt) d ω],

(58)

where I′_m(μ_effa) = ∂I_m (μ_effρ)/∂ρ∣_ρ=a = μ_eff I _m+1(μ_effa)+m/aI_m (μ_effa).

The integral is solved by applying the inverse discrete Fourier transform. The calculations using Eq. (58) are approximately 2 orders of magnitude faster than those calculated without using the ‘damping’ rule.

2.2. Monte Carlo Simulations

The solutions of the diffusion equation are compared to results obtained by Monte Carlo simulations. The Monte Carlo program we developed is able to handle different geometries. For the simulations presented in this study a cylindrical geometry is employed. The light is incident perpendicular to the cylinder top both as a pencil light beam and as a circular flat beam having a finite radius. The refractive index of the cylinder and the surrounding medium is set to n ₀ = n ₁ = 1.0. The Henyey-Greenstein phase function is used for calculating the scattering angles using an anisotropy factor of g = 0.8. 10⁷ photons were used in Monte Carlo simulations.

3. Results

In this section we, first, compare the different solutions of the diffusion equation derived in the steady-state, frequency, and time domains among themselves. Then, we compare the derived solutions of the diffusion theory with Monte Carlo simulations.

3.1. Comparison of the different solutions of the diffusion equation

For the first comparison the cylinder is illuminated by a pencil beam. It is incident onto the cylinder barrel perpendicular to the barrel at a depth of z = 8mm, see inset of Fig. 2. The steady-state transmittance is calculated around the cylinder barrel at a depth of z = 10mm. The height and radius of the cylinder are l_z = 20mm and a = 5mm, respectively. The reduced scattering coefficient is μ′_s = 0.9mm^-1, whereas the absorption coefficient is varied (μ_a = 0.005,0.01,0.015mm^-1). For the calculations shown in Fig. 2 we used Eq. (54) for ω = 0. As expected the transmitted light decreases up to the opposite point relative to the incident beam (at ϕ = π) and the transmittance is smaller for increasing absorption coefficients.

Next, we compare the different solutions of the steady-state diffusion equation for the cylinder shown in Fig. 2. The relative difference of two of the three solutions [(Eqs. (35)–(54)]/[Eq. (35)] is shown in Fig. 3. The relative differences are in the range of 10^-10, thus, the agreement is excellent. The relative differences between the two steady-state solutions calculated via the Hankel transform (version A and version B) are even smaller (data not shown).

Subsequently, the solutions in the time domain are compared. Figure 4 shows the time resolved transmittance from a homogeneous cylinder having a radius of a = 9mm and a height of l_z = 12mm. The cylinder barrel is illuminated at z ₀ = 9mm and the transmittance is detected at z = 9mm on the opposite site of the cylinder (ϕ = π) (black curve) and at the middle of the cylinder bottom (red curve).

Fig. 2. Steady-state radial transmittance around the cylinder barrel at z = 10mm and ρ = a. The cylinder is illuminated perpendicular to the cylinder barrel at z ₀ = 8mm, see inset. The optical parameters are μ′_s = 0.9mm^-1, μ_a = 0.005 (black),0.01 (green),0.015^-1(red)mm , n ₀ = 1.0, and n ₁ = 1.4. The geometrical data of the cylinder are a = 5mm and l_z = 20mm.

Fig. 3. Relative difference of two steady-state solutions for the cylinder considered in Fig. 2.

Fig. 4. Time resolved transmittance from a cylinder that is illuminated at the cylinder barrel at z ₀ = 9mm. The transmittance is detected in z-direction at z = l_z , ρ = 0mm (red curve) and in radial direction at z = 9mm, ρ = a (black curve). The optical and geometrical parameters are μ′_s = 0.8mm^-1, μ_a = 0.03mm^-1, n ₀ = 1.0, n ₁ = 1.4, a = 9mm, and l_z = 12mm. The dashed curves show the transmittance from a homogeneous cube.

Figure 4 shows the results calculated with Eq. (48) (solid curves). The relative difference between the three solutions in the time domain is smaller than 10^-10 (data not shown). For comparison, in Fig. 4 also the time resolved transmittance from a rectangular parallelepiped is shown for both detection positions (dashed curves) [8

A. Kienle, “Light Diffusion Through a Turbid Parallelepiped,” J. Opt. Soc. Am. A 22, 1883–1888 (2005). [CrossRef]

]. The lengths of the edges of the parallelepiped in x-, y-, and z-directions are 9mm, 9mm, and 12mm, respectively. As expected, the results diverge for longer times when the detected photons have experienced the region where the two geometries (cylinder and parallelepiped) are different.

For frequency domain measurements an intensity modulated beam, which consists of a continuous and a sinusoidally oscillating part, is incident onto the turbid medium. In the following we consider only the oscillating part. Figure 5 shows the results of the oscillating part (ω ₀ = 2π · 500MHz) for a pencil beam that is incident onto the cylinder barrel computed with the solution obtained by solving the modified Bessel differential equation. The light transmitted from the barrel at an angle of ϕ = 3π/4 at the same height is shown. The transmittance is computed for different cylinder radii a = 5mm (black curve), a= 6mm (green curve), a = 7mm (red curve).

In order to compare the solution shown in Fig. 5 to the solution obtained with the hyperbolic functions (version B) in the frequency domain, we calculated the transmittance by using the following time harmonic signal

S (\vec{r}, t) = δ (\vec{r} - {\vec{r}}_{0}) \cos (ω_{0 t}) .

(59)

The transmittance is given by

T (\vec{r}, t) = Re {T (\vec{r}, ω = ω_{0}) e^{i ω_{0} t}} .

(60)

The disagreement between the two solutions is small. Thus, we give the numerical values for the above curve having a radius of a = 6mm:

Fig. 5. The transmittance of a sinusoidally oscillating source incident at z = 12mm, ρ = a and detected at the same height at ϕ = 3π/4. The optical and geometrical parameters are μ′_s = 0.7mm^-1, μ_a = 0.02mm^-1, n ₀ = 1.0, n ₁ = 1.4, a = 5,6,7mm, l_z = 22mm, z ₀ = 11mm. The modulation frequency is ω ₀ = 2π · 500MHz.

T (\vec{r}, ω_{0}) = 10^{- 4} {mm}^{- 2} \times {\begin{matrix} 2.5774906387 - 1.2540621701 i \\ 2.5774906366 - 1.2540621700 i \end{matrix} .

(61)

The upper complex number is calculated with the solution derived by solving the modified Bessel differential equation and the lower number with the hyperbolic functions. Thus, the relative differences are in the order of 10^-10.

3.2. Comparison with Monte Carlo simulations

The solutions of the diffusion equation are, first, compared to Monte Carlo simulations for an incident pencil beam, and, then, for an incident flat beam. Figure 6 and Fig. 7 show the spatially resolved reflectance and transmittance from the cylinder top and bottom, respectively. The radius of the cylinders is a = 14mm, whereas the height of the cylinders is l_z = 5mm (blue curves) and l_z = 10mm (red curves). The optical properties are μ′_s = 1.3mm^-1 and μ_a = 0.008mm^-1.

In general, the results of the Monte Carlo simulations (circles) agree well with those obtained from the diffusion theory (solid curves) both for the spatially resolved reflectance and transmittance. For small distances of the spatially resolved reflectance the differences are larger due to the well-known break down of the diffusion theory. At large distances the influence of the border on the light propagation can be seen causing a stronger decrease in the remitted and transmitted intensities. As expected, the increase of the height of the cylinder results in an increase of the remitted light especially at large distances from the source, whereas the transmittance is decreased at small distances and increased at large distances.

Figure 8 shows the spatially resolved reflectance from a homogeneous cylinder that is illuminated in the middle of the top side by a flat beam having radii of ρ_w = 0mm (blue curve), ρ_w = 5mm (green curve), and ρ_w = 10mm (red curve). We note that the beam with ρ_w = 0mm corresponds to a pencil beam. The optical properties of the cylinder are μ′_s = 1.2mm^-1 and μ_a = 0.01mm^-1. The radius and height are a = 20mm and l_z = 5mm, respectively.

Fig. 6. Comparison of the spatially resolved reflectance calculated with the solution of the diffusion equation and with the Monte Carlo method for a turbid cylinder illuminated with a pencil beam in the middle of the cylinder top. The optical and geometrical parameters of the cylinders are μ′_s = 1.3mm^-1, μ_a = 0.008mm^-1, a = 14mm, n ₁ = n ₀ = 1.0, l_z = 5 and 10mm.

Fig. 7. Comparison of the spatially resolved transmittance calculated with the solution of the diffusion equation and with Monte Carlo simulations for a turbid cylinder illuminated with a pencil beam in the middle of the cylinder top. The optical and geometrical parameters of the cylinders are μ′_s = 1.3mm^-1, μ_a = 0.008mm^-1, a = 14mm, n ₁ = n ₀ = 1.0, l_z = 5 and 10mm.

Fig. 8. Spatially resolved reflectance from a homogeneous cylinder that is illuminated by a flat beam with ρ_w = 0,5,10mm. The optical and geometrical parameters are μ′_s = 1.2mm^-1, μ_a = 0.01mm^-1, a = 20mm, l2 = 5mm, and n ₁ = n ₀ = 1.0.

In general, Fig. 8 shows good agreement between diffusion theory and Monte Carlo simulations. Interestingly, the break down of the diffusion equation at small distances can only be seen for ρ_w = 0mm. For large beam diameters the agreement is even good for small distances from the source.

4. Discussion

Solutions of the diffusion equation for a homogeneous cylinder were derived in the steady-state, frequency, and time domains. Comparison of the reflectance and trans-mittance calculated with these solutions showed excellent agreement among themselves in each domain. The speed for calculation of the solutions, however, are different. In the time domain we found for the general case of a beam incident at an arbitrary position that the solution obtained via solving the modified Bessel differential equation is about hundred times faster than the other two solutions. For example, the time needed for obtaining the time resolved transmittance using Eq. (58) is in the order of 10ms on a state-to-the-art computer programmed in Pascal (Delphi).

The derived solutions permit the calculation of the fluence rate for a pencil beam that is incident at an arbitrary position onto the surface of the cylinder. In addition, we derived the solutions for a circular flat beam incident onto the plane surfaces of the turbid cylinder for all three domains.

The solutions of the diffusion equation were compared to Monte Carlo simulations. Good agreement was found as expected from results calculated for turbid media having other geometries. Interestingly, quite a good agreement was obtained for the spatially resolved reflectance at small distances from the center of the incident circular beam having finite beam diameters, although it is expected that diffusion theory breaks down in this case. Probably, the distance dependent errors of the solution of the diffusion equation for an incident pencil beam (at small distances diffusion theory gives larger and at very small distances smaller values than Monte Carlo simulations) cancel out for finite beam diameters, compare Fig. 8.

Acknowledgement

We acknowledge the support by the European Union (nEUROPt, grant agreement no. 201076).

References and links

1.	F. Voit, J. Schäfer, and A. Kienle, “Light Scattering by Multiple Spheres: Comparison between Maxwell Theory and Radiative-Transfer-Theory Calculations,” Opt. Lett. 34, 2593–2595 (2009). [CrossRef] [PubMed]
2.	A. Ishimaru, Wave Propagation and Scattering in Random Media , (Academic Press, New York, 1978).
3.	M. S. Patterson, B. Chance, and B. C. Wilson, “Time Resolved Reflectance and Transmittance for the Noninvasive Measurement of Tissue Optical Properties,” Appl. Opt. 28, 2331–2336 (1989). [CrossRef] [PubMed]
4.	T. J. Farrell, M. S. Patterson, and B. C. Wilson, “A Diffusion Theory Model of Spatially Resolved, Steady-State Diffuse Reflectance for the Noninvasive Determination of Tissue Optical Properties in Vivo ,” Med. Phys. 19, 879–888 (1992). [CrossRef] [PubMed]
5.	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]
6.	S. R. Arridge, M. Cope, and D. T. Delpy, “The Theoretical Basis for the Determination of Optical Pathlength in Tissue: Temporal and Frequency Analysis,” Phys. Med. Biol. 37, 1531–1560 (1992). [CrossRef] [PubMed]
7.	B.W. Pogue and M.S. Patterson, “Frequency Domain Optical Absorption Spectroscopy of Finite Tissue Volumes Using Diffusion Theory,” Phys. Med. Biol. 39, 1157–1180 (1994). [CrossRef] [PubMed]
8.	A. Kienle, “Light Diffusion Through a Turbid Parallelepiped,” J. Opt. Soc. Am. A 22, 1883–1888 (2005). [CrossRef]
9.	A. Sassaroli, F. Martelli, D. Imai, and Y. Yamada, “Study on the Propagation of Ultra-Short Pulse Light in Cylindrical Optical Phantoms,” Phys. Med. Biol. 44, 2747–2763 (1999). [CrossRef] [PubMed]
10.	A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. II. Layered case,” accepted (2010).
11.	http://www.uni-ulm.de/ilm/index.php?id=10020200.
12.	A. Kienle and M. S. Patterson, “Improved Solutions of the Steady-State and the Time-Resolved Diffusion Equations for Reflectance from a Semi-Infinite Turbid Medium,” J. Opt. Soc. Am. A 14, 246–254 (1997). [CrossRef]
13.	R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. McAdams, and B. J. Tromberg, “Boundary Conditions for the Diffusion Equation in Radiative Transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). [CrossRef]
14.	E. Meissel and G.B. Mathews, A Treatise on Bessel Functions and Their Applications to Physics , (Bibliobazaar, 2008).
15.	A. Kienle, M. S. Patterson, N. Dögnitz, R. Bays, G. Wagnières, and H. van den Bergh, “Nonin-vasive Determination of the Optical Properties of Two-Layered Turbid Media,” Appl. Opt. 37, 779–791 (1998). [CrossRef]
16.	H. S. Carslaw and J. C. Jaeger, Conduction of Heat in Solids , (Clarendon, Oxford, 1959).
17.	M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions , (Dover Publications, New York, 1971).
18.	A. Liemert and A. Kienle, “Light Diffusion in N-layered Turbid Media: Frequency and Time Domains,” J. Biomed. Opt. 15, 025002 (2010). [CrossRef] [PubMed]

OCIS Codes
(170.0170) Medical optics and biotechnology : Medical optics and biotechnology
(170.3660) Medical optics and biotechnology : Light propagation in tissues
(170.5280) Medical optics and biotechnology : Photon migration

ToC Category:
Medical Optics and Biotechnology

History
Original Manuscript: February 12, 2010
Revised Manuscript: April 7, 2010
Manuscript Accepted: April 14, 2010
Published: April 21, 2010

Virtual Issues
Vol. 5, Iss. 9 Virtual Journal for Biomedical Optics

Citation
André Liemert and Alwin Kienle, "Light diffusion in a turbid cylinder. I. Homogeneous case," Opt. Express 18, 9456-9473 (2010)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-18-9-9456

Sort: Author | Year | Journal | Reset

References

F. Voit, J. Schäfer, and A. Kienle, “Light Scattering by Multiple Spheres: Comparison between Maxwell Theory and Radiative-Transfer-Theory Calculations,” Opt. Lett. 34, 2593–2595 (2009). [CrossRef] [PubMed]
A. Ishimaru, Wave Propagation and Scattering in Random Media, (Academic Press, New York, 1978).
M. S. Patterson, B. Chance, and B. C. Wilson, “Time Resolved Reflectance and Transmittance for the Noninvasive Measurement of Tissue Optical Properties,” Appl. Opt. 28, 2331–2336 (1989). [CrossRef] [PubMed]
T. J. Farrell, M. S. Patterson, and B. C. Wilson, “A Diffusion Theory Model of Spatially Resolved, Steady-State Diffuse Reflectance for the Noninvasive Determination of Tissue Optical Properties in Vivo,” Med. Phys. 19, 879–888 (1992). [CrossRef] [PubMed]
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]
S. R. Arridge, M. Cope, and D. T. Delpy, “The Theoretical Basis for the Determination of Optical Path length in Tissue: Temporal and Frequency Analysis,” Phys. Med. Biol. 37, 1531–1560 (1992). [CrossRef] [PubMed]
B. W. Pogue, and M. S. Patterson, “Frequency Domain Optical Absorption Spectroscopy of Finite Tissue Volumes Using Diffusion Theory,” Phys. Med. Biol. 39, 1157–1180 (1994). [CrossRef] [PubMed]
A. Kienle, “Light Diffusion through a Turbid Parallelepiped,” J. Opt. Soc. Am. A 22, 1883–1888 (2005). [CrossRef]
A. Sassaroli, F. Martelli, D. Imai, and Y. Yamada, “Study on the Propagation of Ultra-Short Pulse Light in Cylindrical Optical Phantoms,” Phys. Med. Biol. 44, 2747–2763 (1999). [CrossRef] [PubMed]
A. Liemert, and A. Kienle, “Light diffusion in a turbid cylinder. II. Layered case,” accepted (2010).
http://www.uni-ulm.de/ilm/index.php?id=10020200.
A. Kienle, and M. S. Patterson, “Improved Solutions of the Steady-State and the Time-Resolved Diffusion Equations for Reflectance from a Semi-Infinite Turbid Medium,” J. Opt. Soc. Am. A 14, 246–254 (1997). [CrossRef]
R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. McAdams, and B. J. Tromberg, “Boundary Conditions for the Diffusion Equation in Radiative Transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). [CrossRef]
E. Meissel, and G. B. Mathews, A Treatise on Bessel Functions and Their Applications to Physics, (Bibliobazaar, 2008).
A. Kienle, M. S. Patterson, N. Dögnitz, R. Bays, G. Wagnières, and H. van den Bergh, “Noninvasive Determination of the Optical Properties of Two-Layered Turbid Media,” Appl. Opt. 37, 779–791 (1998). [CrossRef]
H. S. Carslaw, and J. C. Jaeger, Conduction of Heat in Solids, (Clarendon, Oxford, 1959).
M. Abramowitz, and I. A. Stegun, Handbook of Mathematical Functions, (Dover Publications, New York, 1971).
A. Liemert, and A. Kienle, “Light Diffusion in N-layered Turbid Media: Frequency and Time Domains,” J. Biomed. Opt. 15, 025002 (2010). [CrossRef] [PubMed]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper