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], whereas on the mesoscopic scale transport equation and on the macroscopic scale diffusion equation are normally applied [2].

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], further for a sphere [6, 7], for a parallel-epiped [8], and for a cylinder [6, 7, 9].

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

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

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

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

 

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

Download Full Size | PDF

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

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

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

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

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

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]. For the finite Hankel transform the relation

is valid [14]. 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, ω)

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

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]. 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]. Equation (8) is solved by using the solution of Eq. (9) given below in version A and B. Thus, Φ can be obtained with

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

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

By applying the relation

we get

and, finally, we obtain

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

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

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

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

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

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

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

is used to obtain

For the special case of central incidence we get

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)

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

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

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

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

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

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

We note that Eq. (31) is also used for solving the diffusion equation for N-layered turbid media, see the accompanying paper [10]. 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

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

where

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

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

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

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:

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

we obtain from Eq. (8)

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

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

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

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

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

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

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 $\alpha =\sqrt{{\mu }_a/D+s_n^2+\mathrm{i\omega }/\mathrm{Dc}}$. By employing the Fourier transform

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

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

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)

where $\alpha =\sqrt{{\mu }_a/D+\mathrm{i\omega }/\mathrm{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

I_m and K_m denote the modified Bessel functions of first and second kind of m^th order. By setting $\alpha =\sqrt{{\mu }_a/D+{\lambda }_k^2+\mathrm{i\omega }/\mathrm{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

where

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

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

and the transmittance is obtained with

For the solution in the time domain, the frequency domain solution, Φ(r⃗, ω), is used and split into real and imaginary parts [17]. 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]. The result is

and the transmittance is obtained with

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.

Download Full Size | PDF

 

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

Download Full Size | PDF

 

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.

Download Full Size | PDF

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

The transmittance is given by

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.

Download Full Size | PDF

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.

Download Full Size | PDF

 

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.

Download Full Size | PDF

 

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.

Download Full Size | PDF

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.