When the telegrapher ’ s equation furnishes a better approximation to the transport equation than the diffusion approximation

It has been suggested that a solution to the transport equation which includes anisotropic scattering can be approximated by the solution to a telegrapher’s equation @A.J. Ishimaru, Appl. Opt. 28, 2210~1989!#. We show that in one dimension the telegrapher’s equation furnishes an exact solution to the transport equation. In two dimensions, we show that, since the solution can become negative, the telegrapher’s equation will not furnish a usable approximation. A comparison between simulated data in three dimensions indicates that the solution to the telegrapher’s equation is a good approximation to that of the full transport equation at the times at which the diffusion equation furnishes an equally good approximation. @S1063-651X~97!04205-0#


I. INTRODUCTION
Problems related to the transport of matter in disordered media are common to a number of areas of science and technology, and have given rise to a considerable amount of research aimed at delineating properties of motion in such media ͓1,2͔.When the medium is a continuum, a natural starting point for any analysis is based on solving an appropriate transport equation.Since there are no general solutions other than purely numerical ones for such equations it is often difficult to fit experimental data to theory, as would be necessary to determine physically significant parameters.Hence a number of approximation schemes have been used to derive more tractable mathematical models.
The simplest such model which is adequate in many applications approximates the full transport equation by a diffusion equation.A related model which has recently been applied to problems arising in tissue optics is based on the theory of lattice random walks ͓3-5͔.Both of these are deficient in failing to properly account for anisotropic scattering, although this problem tends to vanish at increasingly long times, which is essentially the regime in which the central limit theorem is equivalent to the diffusion approximation.Nevertheless the diffusion equation is appealing because it is trivial to solve in unbounded spaces and easy to specify boundary conditions associated with finite media.
Several compromises have been suggested to incorporate at least some aspects of anisotropic scattering into the diffusion picture.A strategy used in several optical applications is to correct the diffusion constant by using a characteristic parameter related to the degree of anisotropic scattering ͓6͔.
A recent suggestion has been made by Ishimaru in the context of tissue optics, that anisotropy can be incorporated into the analysis by replacing the diffusion equation by a telegrapher's equation ͑TE͒ ͓7͔, which has the form where T is a parameter with the dimensions of time and c is a speed.Two motivating factors for Ishimaru's suggestion are to be found in an observation by Goldstein that the continuum limit of a persistent random walk is a telegrapher's equation ͓8͔, and from an approximation based on the transport equation derived by Ishimaru ͓9͔.
In the present paper we investigate whether the TE provides a better approximation to the solution of the full transport equation over a physically relevant time span than does the diffusion equation.We show that this is true for the single-speed model in one dimension, but that in two dimensions the TE has a solution that is negative in some region of the plane, and therefore cannot be identified as a probability density.In the case of three dimensions we show, by comparing the theoretical results to data simulated using a physically plausible model, that the solution to the TE is not a good approximation to the simulated results except at comparatively long times, at which the solution of the TE essentially coincides with the solution to the diffusion equation.

II. TRANSPORT EQUATION
In modeling transport we will consider the single-speed model only ͓10͔.Let be the scattering rate.Absorption in the medium will be modeled in terms of Beer's law characterized by a rate parameter ͓5͔.Only the case of a pulse propagating in a translationally invariant medium will be considered.Let the direction in which a single particle moves be denoted by the solid angle ⍀.In three dimensions ⍀ϭ (,), where is the polar angle and is the azimuthal angle.Further, let ␤(⍀͉⍀Ј) denote the scattering kernel, i.e., the probability that a single scattering changes the direction of the photon path from ⍀Ј to a direction falling somewhere in the infinitesimal interval (⍀,⍀ϩd⍀) is equal to ␤(⍀͉⍀Ј)d⍀.The speed between collisions will be taken to be a constant v, and the corresponding velocity, which takes the direction into account, will be denoted by v(⍀).For example, in three dimensions v͑ ⍀͒ϭ͑vsin cos ,vsin sin ,vcos ͒, where and are the polar and azimuthal angles, respectively.The object of our investigation is to find the probabil- ity density for the particle to be at r moving in the direction ⍀ at time t.This function will be denoted by p(r,⍀,t).In this notation the full transport equation is ‫ץ‬p/‫ץ‬t ϭϪ͑ϩ͒pϪv͑⍀͒•"p ϩ ͵ p͑r,⍀Ј,t͒␤͑⍀͉⍀Ј͒d⍀Ј. ͑2͒ Only the case in which the initial condition is a single pulse will be analyzed.In three dimensions this allows us to write p͑r,⍀Ј,0͒ϭ ͑ 1/4͒ ␦͑r͒ϭ ␦͑r͒/16 2 r 2 .

͑3͒
We will assume that the kernel ␤(⍀͉⍀Ј) depends only on the deflection angle ␥, i.e., the angle between directions ⍀ and ⍀Ј.In this case, the deviation from isotropic scattering will be denoted by gϭ͗cos␥͘.The solution to Eq. ͑2͒ can be expanded in terms of spherical harmonics, although in practice the series is generally truncated.The so-called P 1 approximation, commonly used in nuclear reactor analysis ͓10,11͔, is equivalent to the assumption that p(r,⍀,t) is only weakly dependent on angle.In three dimensions this is equivalent to writing ͓11͔ p͑r,⍀,t ͒Ϸ ͑ 1/4͒ P͑r,t͒ϩ ͑3/4v 2 ͒ v"⍀)•J͑r,t͒.͑4͒ Here P(r,t) is the angle-averaged value of p(r,⍀,t), i.e., P͑r,t ͒ϭ ͵ p͑r,⍀,t ͒d⍀, and the flux J(r,t) is related to p(r,⍀,t) by J͑r,t ͒ϭ ͵ v͑ ⍀͒p͑r,⍀,t ͒d⍀.

͑5͒
The angle variables appear linearly in Eq. ͑4͒ as indicated.
On substituting this expansion into Eq.͑2͒, we find a coupled set of equations for P(r,t) and J(r,t), which is ‫ץ‬J/‫ץ‬t ϭϪ͑Јϩ͒JϪ ͑v 2 /3͒"P, ͑6b͒ where is an anisotropy-corrected scattering rate.On the further assumption that ‫ץ‬J/‫ץ‬t can be neglected in comparison to the terms on the right-hand side of Eq. ͑6b͒, one finds a damped diffusion equation for P(r,t), ͑‫ץ‬P/‫ץ‬t͒ ϩPϭD" 2 P, with an effective diffusion constant If the time derivative of the flux is not considered negligible then the function J(r,t) can be eliminated from Eq. ͑6b͒, and the equation satisfied by P(r,t) is the damped telegrapher's equation where the velocity c is equal to v/ͱ3.This equation has been obtained by taking the divergence on both sides of Eq. ͑6b͒, commuting the divergence operator and ‫‪t‬ץ/ץ‬ on the left-hand side, and finally using Eq.͑6a͒ to write "•J in terms of P. We note that the substitution P(r,t)ϭQ(r,t)e Ϫt transforms Eq. ͑8͒ into the standard form of the TE in Eq. ͑1͒ with the time parameter Tϭ1/Ј.We also observe that the solution to the telegrapher's equation at very long times reduces to a Gaussian ͓12͔, which is the solution to the diffusion equation in free space.The question to be addressed is whether the solution of the telegrapher's equation more faithfully reproduces the solution of the full transport equation than does the diffusion equation.We consider this problem in one, two, and three dimensions.

A. One dimension
In one dimension a particle can move in only two directions, the coordinate of its location at any time t either increasing or decreasing with a constant speed.Thus a description of the system requires consideration of the evolution of two probability densities, p i (x,t), iϭϮ1, where iϭϮ1 refers to a particle having a speed Ϯv.If P(x,t) and J(x,t) are defined by then the full transport equation coincides exactly with the P 1 approximation.Hence it follows from our earlier remark that the equation satisfied by P(x,t) is the damped telegrapher's equation, Eq. ͑8͒, with cϭv.Hence in one dimension the TE provides an exact, rather than approximate, solution to the full transport equation.This is also implied in the work of Goldstein ͓8͔.

B. Two dimensions
The problem of the appropriateness of the TE as a model accounting for anisotropy is much more interesting in dimensions greater than one, since the P 1 approximation is not equivalent to the full transport equation.In two dimensions the angular coordinate ⍀ consists of a single polar angle .
The P 1 approximation is found by retaining only the terms nϭϪ1, 0, and 1 in the Fourier expansion of p(r,,t), p͑r,,t

͑13͒
In this equation H(z) is the Heaviside step function, H(z)ϭ1, zϾ0, and H(z)ϭ0 for zϽ0.
Figure 1 compares profiles of the angular density 2rP(r,t) for the isotropic model where gϭ0.The three values of P(r,t) that are plotted are P exact (r,t), the exact solution, P TE (r,t), the solution obtained from Eq. ͑13͒; and P diffusion (r,t), the solution obtained from the diffusion approximation.The exact solution was obtained in ͓13͔ and reads Note that the exact result, Eq. ͑14͒, and the solution to TE, Eq. ͑13͒, both converge when t→ϱ to the diffusion approximation given by P diffusion ͑ r,t ͒ϭ͑ /2v 2 t͒ exp͑Ϫ r 2 /2v 2 t͒.͑15͒ We note that the divergence of P(r,t) at rϭvt/ͱ2 ͓c.f.Eq. ͑13͒ is due to the time propagation of the singular initial condition Eq. ͑11͒.The same happens to P exact (r,t) at the boundary rϭvt.We also observe that the divergence is integrable, and both P(r,t) and P exact (r,t) are normalized to 1.
It is evident that the solution for P TE (r,t) becomes negative somewhere in the interval (0,vt/ͱ2), and therefore cannot be a probability density.Such behavior is a characteristic of solutions found in even-dimensional spaces ͓12͔.In consequence, the TE cannot furnish an accurate approximation to the solution of the full transport equation in even dimensions, except in the limit of times so large that the solution of the TE is essentially equal to the solution provided by the diffusion equation.In addition, even in the regions where the solution of the TE is positive, numerical data indicate that the diffusion approximation is better than the TE solution ͑see Fig. 1͒.

C. Three dimensions
An expression is known for the Fourier-Laplace transform of the solution of the isotropic transport equation in Eq. ͑2͒ for the fully isotropic case ͓13͔, ␤͑⍀͉⍀Ј͒ϭ 1/4 .͑16͒ Since the transform cannot be inverted except numerically, we conducted a simulation study to test the quality of the approximation furnished by the TE.Our simulations were only for particles moving in an unbounded space.The time between two successive scattering events were random variables whose properties were described by a negative exponential probability density ͑t͒ϭe Ϫt .͑17͒ In our choice of the form of the kernel ␤(⍀͉⍀Ј), we assumed that this depended only on the deflection angle ␥, while the azimuthal angle was assumed to be uniformly distributed in the interval (Ϫ,).Between scatterings the particles were assumed to move in straight lines.The effect of anisotropic scattering was mimicked by choosing the scattering angle according to the widely used Henyey-Greenstein phase function ͓14͔ 1.The two-dimensional radial probability densities 2rP exact (r,t) ͑solid line͒, 2rP TE (r,t) ͑dashed line͒, and 2rP diffusion (r,t) ͑dot-dashed line͒ plotted as a function of r for ͑a͒ tϭ1 and ͑b͒ tϭ10.Both sets of curves were generated using the parameter values ϭvϭ1.
Other phase functions may have been used, because the only requirement is that the anisotropic factor gϭ͗cos␥͘ 0 is different than zero.
The data generated in this way were compared with results generated from the solution to the three-dimensional TE in free space.Analogously to the two-dimensional case, this solution is expressed in terms of the variable ϭ ͑ Ј/2v͒ͱv 2 t 2 Ϫ3r 2 ͑19͒ as where I n (z) are modified Bessel functions.At long times the solution given in Eq. ͑20͒ tends toward the Gaussian solution of the diffusion equation, i.e., P͑r,t ͒→͑ 3Ј/4v 2 t͒ 3/2 exp͑Ϫ 3Јr 2 /4v 2 t͒.͑21͒ Figure 2 shows the results of 10 5 simulated runs as described above and carried out for the parameter values gϭ0.5, vϭ1, and ϭ1.Two sets of data are shown that correspond to times tϭ5 and 25, which is in the regime in which the Gaussian form in Eq. ͑21͒ is expected to be valid.The agreement with the asymptotic Gaussian form of P(r,t) at the larger time is seen to be quite good, while the agreement with the TE at shorter times is not.The solution given in Eq. ͑20͒ is non-negative, but the fact that the rescaled speed of propagation cannot correspond to the true speed of the scattered particles prevents finding a more accurate approximation, at least in terms of the TE.This is apparently due to the contributions to P(r,t) from particles which have never been scattered.It is possible that this shortcoming can be overcome by explicitly decomposing the solution of the transport equation into a contribution from the unscattered particles, and one from those that were scattered at least once, as is sometimes done in other problems related to photon diffusion in turbid media ͓7͔.A study of this type of approximation is presently under consideration.

III. A FINAL COMMENT
Both simulations and analytical considerations suggest that, except in one dimension, there are significant difficulties in basing an approximate solution to the full transport equation on an ''equivalent'' telegrapher's equation.In two dimensions this occurs because the solution to the telegrapher's equation has regions in which the solution is negative.In three dimensions the difficulty can be traced to an inaccurate accounting for particles that remain unscattered at time t.At very long times, when the number of such particles is considerably reduced, the solution to the telegrapher's equation is quite accurately approximated by the solution to a diffusion equation.Our simulations suggest that in this regime the solution to the full transport equation can also be modeled quite accurately.