text:
# Multiple scattering

## 1 Applications and overview

As a quantitative analytical technique multiple scattering is used to study clouds , pollution, chemichals in solution , radiation in the human body , etc… A great effort has been made to study the scattering of laser beams by clouds. To this end some authors  have developped numerical stochastic models, based on the same general philosophy as ours, which have been shown to be equivalent to the transport equation which is more commonly used. Simulations are a major tool in the study of multiple scattering , but real situations are computationally so demanding that numerical methods based on analitical understanding of the problem can be necessary.

In Physics multiple scattering is obviously related to Lorentz gases  and to diffusion of neutrons in nuclear reactors . It has many applications in solid state . In most situations the particle being scattered is a photon or an electron. Here we shall work with a scalar particle, but the methods can be generalized for the case of internal degrees of freedom as long as interference effects remain negligible.

In the 2nd section we find the directions taken by the particles of a beam which encounters randomly distributed diffusing centers. Conversely the differential cross section of the particle-scatterer interaction can be solved in terms of the distribution of directions. Asymptotics are also discussed. The results of the 2nd section are experimentally applicable as long as the detectors are placed at a distance to the medium to be studied which is much larger than the diameter of the medium. In the 3rd section the motion of the center of mass is found. In the 4th section we study the case of isotropic cross section. The probability density for the particles is found explicitely in 2 dimensions and in 3 dimensions its Fourier-Laplace transform is found. In the 5th section we show how to calculate the Fourier-Laplace transform of the probability density with arbitrary precision in the general anisotropic case for both 2 and 3 dimensions.

The initial condition for all the results in this article is a beam of particles. This is the general case in the sense that any initial condition can be decomposed in terms of beams. The mathematical setting of this work is that of integral equations. For some instances this work provides an alternative to the Boltzmann equation for the study of transport processes.

## 2 Diffusion of directions

Consider the repeated elastic scattering of a point-like particle by randomly distributed recoilless diffusing centers. We assume that the scattering itself takes place in a region of negligible volume (or area) in the sense that the trajectory of the particle can be well approximated by a sequence of straight line segments. Let $`f`$ be the angular cross section, where $`f`$ is a function of the angle. If $`f`$ is not equal to a constant, i. e., if the scattering is not isotropic, then a statistical description of the multiple scattering becomes difficult because the process is not Markovian in the physical space. However, the stochastic process for the directions is Markovian. In fact if, in $`d`$ dimensions, we represent the direction along which the particle is moving by a point on $`S^{d1}`$, then the probability of transition from one point of $`S^{d1}`$ to another will depend only on the angle between them, so that convolution techniques can be applied. We will show this in 2 and then in 3 dimensions.

In 2 dimensions let $`F(t,\phi )`$ be the probability density for the particle to be moving along the direction $`\phi `$ at time $`t`$. $`\stackrel{~}{}`$ will denote Fourier transform and $`\overline{\phi }`$ will be the Fourier frequency associated to the angle $`\phi `$. If $`\lambda `$ is the probability of collision per unit time then in an infinitesimal time $`dt`$ the probability density $`F`$ will change to:

$$F(t+dt,\phi )=(1\lambda dt)F(t,\phi )+\lambda dt𝑑\phi ^{}F(t,\phi ^{})f(\phi \phi ^{}),$$
(1)
The Fourier transform of the above equation is:

$$\stackrel{~}{F}(t+dt,\overline{\phi })=(1\lambda dt)\stackrel{~}{F}(t,\overline{\phi })+\lambda dt\stackrel{~}{F}(t,\overline{\phi })\stackrel{~}{f}(\overline{\phi })=\stackrel{~}{F}(t,\overline{\phi })+\lambda dt\stackrel{~}{F}(t,\overline{\phi })(\stackrel{~}{f}(\overline{\phi })1).$$
(2)
This yields the differential equation:

$$\frac{}{t}\mathrm{ln}\stackrel{~}{F}(t,\overline{\phi })=\lambda [\stackrel{~}{f}(\overline{\phi })1],$$
(3)
whose solution is:

$$\stackrel{~}{F}(t,\overline{\phi })=\stackrel{~}{F}(0,\overline{\phi })\mathrm{exp}\lambda [\stackrel{~}{f}(\overline{\phi })1]t.$$
(4)
An analogous reasoning shows that when $`\lambda `$ or $`f`$ are time dependent:

$$\stackrel{~}{F}(t,\overline{\phi })=\stackrel{~}{F}(0,\overline{\phi })\mathrm{exp}_0^t𝑑t\lambda (t)[\stackrel{~}{f}(t,\overline{\phi })1].$$
(5)
These formulae are valid for any cross section $`f`$. In particular left-right symmetry ($`f(\phi )=f(2\pi \phi )`$) is not required.

When all particles move along the direction $`\phi =0`$ at $`t=0`$, $`\stackrel{~}{F}(0,\overline{\phi })=1`$ and the above formulae simplify. This kind of initial condition, which we shall call stream initial condition, is important because an arbitrary initial condition can be constructed in terms of it. In 3 dimensions the stream initial condition will be a beam of particles moving towards the positive direction of the $`z`$ axis, due to the definition of polar coordinates in 3 dimensions.

In 3 dimensions the equation corresponding to (1) is:

$$F(t+dt,\mathrm{\Omega })=(1\lambda dt)F(t,\mathrm{\Omega })+\lambda dt𝑑\mathrm{\Omega }^{}F(t,\mathrm{\Omega }^{})f(\mathrm{\Omega }\mathrm{\Omega }^{}),$$
(6)
where $`\mathrm{\Omega }\mathrm{\Omega }^{}`$ is the angle between the directions $`\mathrm{\Omega }`$ and $`\mathrm{\Omega }^{}`$. Here we are going to assume that the scattering probability $`f`$ depends only on the polar angle $`\theta `$, or, in other words, that if the particle undergoing scattering has internal (pseudo)vectorial degrees of freedom (such as the polarization of a photon), these are not involved in the scattering or are averaged over. We also assume, for simplicity’s sake, that the initial condition $`F(0,\mathrm{\Omega })`$ is independent of $`\phi `$. Then $`F(t,\mathrm{\Omega })`$ is also going to be $`\phi `$ independent and both $`F(t,\mathrm{\Omega })`$ and $`f(\mathrm{\Omega })`$ can be expanded in terms of Legendre polynomials $`P_l`$ as follows (see  or ):

$$F(t,\mathrm{\Omega })=\underset{l=0}{\overset{\mathrm{}}{}}\sqrt{\frac{2l+1}{4\pi }}F_l(t)Y_l^0(\theta )$$
(7)
and

$$f(\mathrm{\Omega })=\underset{l=0}{\overset{\mathrm{}}{}}\sqrt{\frac{2l+1}{4\pi }}f_lY_l^0(\theta ),$$
(8)
where

$$Y_l^0(\theta )=\sqrt{\frac{2l+1}{4\pi }}P_l(\mathrm{cos}\theta )$$
(9)
are the spherical harmonics and

$$F_l(t)=𝑑\mathrm{\Omega }\sqrt{\frac{4\pi }{2l+1}}F(t,\mathrm{\Omega })Y_l^0(\mathrm{\Omega }),$$
(10)
and a similar formula for $`f_l`$.

If we now substitute these expansions in the convolution product in (6) and use the sum formula for Legendre polynomials,

$$Y_l^0(\mathrm{\Omega }\mathrm{\Omega }^{})=\sqrt{\frac{4\pi }{2l+1}}\underset{m=l}{\overset{l}{}}Y_l^m(\mathrm{\Omega })Y_l^m(\mathrm{\Omega }^{}),$$
(11)
and

$$𝑑\mathrm{\Omega }Y_l^m(\mathrm{\Omega })Y_l^{}^m^{}(\mathrm{\Omega })=\delta _{ll^{}}\delta _{mm^{}},$$
(12)
we obtain:

$$𝑑\mathrm{\Omega }^{}F(t,\mathrm{\Omega }^{})f(\mathrm{\Omega }\mathrm{\Omega }^{})=\underset{l=0}{\overset{\mathrm{}}{}}\sqrt{\frac{2l+1}{4\pi }}F_l(t)f_lY_l^0(\mathrm{\Omega }).$$
(13)
Proceeding as in the 2-dimensional case we obtain the time evolution of the coefficients $`F_l(t)`$:

$$F_l(t)=F_l(0)\mathrm{exp}\lambda (f_l1)t,l=0,1,2,\mathrm{},$$
(14)
as well as an equation analogous to (5). For the $`z`$-stream initial condition, $`F_l(0)=1`$ for $`l=0,1,2,\mathrm{}`$. Different derivations of eq. (14) have been given by Goudsmit and Saunderson (,) and by Lewis  (these derivations are also reviewed in ).

Note that the fact that $`f(\theta )`$ is arbitrary allows us to treat cross sections like the one of the Rayleigh scattering, which a Fokker-Planck equation on the sphere (e. g. the Debye rotational diffusion equation ) can not approximate.

It is intuitive to think that asymptotically the directions will be uniformly distributed. From

$$F(t,\phi )=\frac{1}{2\pi }\underset{\overline{\phi }ϵZ}{}\stackrel{~}{F}(0,\overline{\phi })e^{\lambda [\stackrel{~}{f}(\overline{\phi })1]t}e^{i\overline{\phi }\phi },$$
(15)
and

$$\stackrel{~}{f}(\overline{\phi })=_0^{2\pi }𝑑\phi f(\phi )e^{i\overline{\phi }\phi },$$
(16)
we see that unless $`f(\phi )`$ is a set of Dirac deltas placed at the vertices of a regular polygon, our final distribution in 2 dimensions is indeed going to be uniform.

In 3 dimensions normalization implies $`f_0=1`$. We now show that the other coefficients, $`f_l,l=1,2,\mathrm{}`$, are smaller than 1:

$$f_l=𝑑\mathrm{\Omega }\sqrt{\frac{4\pi }{2l+1}}f(\mathrm{\Omega })Y_l^0(\mathrm{\Omega })=2\pi _1^1d(\mathrm{cos}\theta )f(\theta )P_l(\mathrm{cos}\theta ).$$
(17)
The Legendre polynomials take absolute values always smaller than 1 except for $`\mathrm{cos}\theta =\pm 1`$, where the even ones take the value 1 and the odd ones take the value $`\pm 1`$. Thus if $`f(\theta )\alpha \delta (\theta )+\overline{\alpha }\delta (\theta \pi )`$, where $`\alpha +\overline{\alpha }=1`$, then

$$f_l<2\pi _1^1d(\mathrm{cos}\theta )f(\theta )=𝑑\mathrm{\Omega }f(\mathrm{\Omega })=1,l=1,2,\mathrm{},$$
(18)
and the time evolution (14) implies that in 3 dimensions the final distribution $`F(\mathrm{},\mathrm{\Omega })`$ will be the uniform distribution.

To measure quantitatively how the distribution approaches its final state, we choose the $`L^2`$ distance between the current distribution and the final one. In 2 dimensions this is:

$$_0^{2\pi }𝑑\phi \left(F(t,\phi )\frac{1}{2\pi }\right)^2.$$
(19)
The Fourier transform is an isometry of the $`L^2`$ norm, so that the above distance is equal to:

$$\frac{1}{2\pi }\underset{\overline{\phi }0}{}e^{2\lambda [\stackrel{~}{f}(\overline{\phi })1]t}.$$
(20)
The time dependence of the second cumulant of the probability density for finding a particle which is suffering multiple scattering tends to be linear as $`t\mathrm{}`$. By studying the approach to this behaviour Boguñá et alii  found the above formula with the sum restricted to $`\overline{\phi }=\pm 1`$. Higher terms in (20) are related to higher moments.

In 3 dimensions the $`L^2`$ distance between the current distribution, $`F(t,\mathrm{\Omega })`$, and the final one, $`\frac{1}{4\pi }`$, can be expressed in a form analogous to (20) using expansion (7), the orthonormality relation (12) and the evolution equation (14). The final expression for the $`L^2`$ distance is:

$$\underset{l=1}{\overset{\mathrm{}}{}}\frac{2l+1}{4\pi }F_l(0)^2e^{2\lambda (f_l1)t}$$
(21)
## 3 Motion of the center of mass

If the cross section $`f`$ is not ballistic then at every collision some of the initial momentum of the particle is transfered to the scatterer. As a result the velocity of the center of mass goes exponentially to zero. It is interesting to find how far the center of mass will go before (at $`t=\mathrm{})`$ it stops. To study this note that in 2 dimensions for a particle (which we assume of unit mass) impinging along the direction $`\phi `$ upon a scatterer, its momentum changes from $`(v\mathrm{cos}\phi ,v\mathrm{sin}\phi )`$ to, on the average, $`(v_0^{2\pi }𝑑\phi ^{}f(\phi ^{}\phi )\mathrm{cos}\phi ^{},v_0^{2\pi }𝑑\phi ^{}f(\phi ^{}\phi )\mathrm{sin}\phi ^{})`$. Thus Newton’s law for the center of mass reads:

$$\frac{d\stackrel{}{p}}{dt}=\lambda v$$

$$_0^{2\pi }𝑑\phi F(t,\phi )(_0^{2\pi }𝑑\phi ^{}f(\phi ^{}\phi )\mathrm{cos}\phi ^{}\mathrm{cos}\phi ,_0^{2\pi }𝑑\phi ^{}f(\phi ^{}\phi )\mathrm{sin}\phi ^{}\mathrm{sin}\phi ),$$
(22)
or, using (15),

$$\frac{d\stackrel{}{p}}{dt}=\frac{\lambda v}{2\pi }\underset{\overline{\phi }ϵZ}{}\stackrel{~}{F}(0,\overline{\phi })e^{\lambda [\stackrel{~}{f}(\overline{\phi })1]t}$$

$$_0^{2\pi }𝑑\phi e^{i\mu \phi }(_0^{2\pi }𝑑\phi ^{}f(\phi ^{}\phi )\mathrm{cos}\phi ^{}\mathrm{cos}\phi ,_0^{2\pi }𝑑\phi ^{}f(\phi ^{}\phi )\mathrm{sin}\phi ^{}\mathrm{sin}\phi )=$$

$$\frac{\lambda v}{2\pi }\underset{\overline{\phi }ϵZ}{}\stackrel{~}{F}(0,\overline{\phi })e^{\lambda [\stackrel{~}{f}(\overline{\phi })1]t}$$

$$(\stackrel{~}{f}(\overline{\phi })\pi [\delta _{\overline{\phi },1}+\delta _{\overline{\phi },1}]\pi [\delta _{\overline{\phi },1}+\delta _{\overline{\phi },1}],\stackrel{~}{f}(\overline{\phi })\frac{\pi }{i}[\delta _{\overline{\phi },1}\delta _{\overline{\phi },1}]\frac{\pi }{i}[\delta _{\overline{\phi },1}\delta _{\overline{\phi },1}])=$$

$$\frac{v}{2}\underset{\overline{\phi }ϵZ}{}\stackrel{~}{F}(0,\overline{\phi })\frac{d}{dt}e^{\lambda [\stackrel{~}{f}(\overline{\phi })1]t}([\delta _{\overline{\phi },1}+\delta _{\overline{\phi },1}],i[\delta _{\overline{\phi },1}\delta _{\overline{\phi },1}])$$
(23)
Since the initial condition $`F(0,t)`$ and the scattering cross section $`f(\phi )`$ are real valued functions, $`\stackrel{~}{F}(0,1)=\stackrel{~}{F}(0,1)^{}`$ and $`\stackrel{~}{f}(1)=\stackrel{~}{f}(1)^{}`$. Therefore

$$\frac{d\stackrel{}{p}}{dt}=\mathrm{Re}\left[v\stackrel{~}{F}(0,1)\frac{d}{dt}e^{\lambda [\stackrel{~}{f}(1)1]t}(1,i)\right]$$
(24)
For the initial condition $`\stackrel{}{p}(0)=(v_0^{2\pi }𝑑\phi F(0,\phi )\mathrm{cos}\phi ,v_0^{2\pi }𝑑\phi F(0,\phi )\mathrm{sin}\phi )`$ the solution to this differential equation is:

$$\stackrel{}{p}(t)=\mathrm{Re}\left[v\stackrel{~}{F}(0,1)e^{\lambda [\stackrel{~}{f}(1)1]t}(1,i)\right]$$
(25)
In particular for the stream initial condition and symmetric cross section,

$$\stackrel{}{p}(t)=ve^{\lambda [\stackrel{~}{f}(1)1]t}\stackrel{}{i}.$$
(26)
Finally, for $`\stackrel{}{x}(0)=\stackrel{}{0}`$ a last integration yields:

$$\stackrel{}{x}(t)=\mathrm{Re}\left[v\stackrel{~}{F}(0,1)\frac{(e^{\lambda [\stackrel{~}{f}(1)1]t}1)}{\lambda [\stackrel{~}{f}(1)1]}(1,i)\right].$$
(27)
For the stream case with symmetric cross section this implies that the center of mass penetrates a distance

$$|\stackrel{}{x}(\mathrm{})|=\frac{v}{\lambda [1\stackrel{~}{f}(1)]}.$$
(28)
into the scattering medium.

In 3 dimensions for a particle of unit mass impinging along the direction $`\mathrm{\Omega }(\theta ,\phi )`$ upon a scatterer, its momentum changes from $`(v\mathrm{sin}\theta \mathrm{cos}\phi ,v\mathrm{sin}\theta \mathrm{sin}\phi ,v\mathrm{cos}\theta )`$ to, on the average, $`(v_0^\pi 𝑑\theta ^{}\mathrm{sin}\theta ^{}_0^{2\pi }𝑑\phi ^{}f(\theta ^{}\theta ,\phi ^{}\phi )\mathrm{sin}\theta ^{}\mathrm{cos}\phi ^{},v_0^\pi 𝑑\theta ^{}\mathrm{sin}\theta ^{}_0^{2\pi }𝑑\phi ^{}f(\theta ^{}\theta ,\phi ^{}\phi )\mathrm{sin}\theta ^{}\mathrm{sin}\phi ^{},v_0^\pi 𝑑\theta ^{}\mathrm{sin}\theta ^{}_0^{2\pi }𝑑\phi ^{}f(\theta ^{}\theta ,\phi ^{}\phi )\mathrm{cos}\theta ^{})`$. Thus Newton’s law for the center of mass reads:

$$\frac{d\stackrel{}{p}}{dt}=\lambda v_0^\pi 𝑑\theta \mathrm{sin}\theta _0^{2\pi }𝑑\phi F(t,\theta ,\phi )$$

$$(_0^\pi d\theta ^{}\mathrm{sin}\theta ^{}_0^{2\pi }d\phi ^{}f(\theta ^{}\theta ,\phi ^{}\phi )\mathrm{sin}\theta ^{}\mathrm{cos}\phi ^{}\mathrm{sin}\theta \mathrm{cos}\phi ,$$

$$_0^\pi 𝑑\theta ^{}\mathrm{sin}\theta ^{}_0^{2\pi }𝑑\phi ^{}f(\theta ^{}\theta ,\phi ^{}\phi )\mathrm{sin}\theta ^{}\mathrm{sin}\phi ^{}\mathrm{sin}\theta \mathrm{sin}\phi ,$$

$$_0^\pi d\theta ^{}\mathrm{sin}\theta ^{}_0^{2\pi }d\phi ^{}f(\theta ^{}\theta ,\phi ^{}\phi )\mathrm{cos}\theta ^{}\mathrm{cos}\theta ).$$
(29)
As in section 2 we restrict ourselves to $`F`$ and $`f`$ functions of $`\theta `$ only. Then only the $`z`$ component of the above vector is non zero. Using the expression (14) for the time evolution of the Legendre polynomials coefficients of $`F(t,\theta )`$ we obtain:

$$\frac{d\stackrel{}{p}}{dt}=\lambda v_0^\pi 𝑑\theta \mathrm{sin}\theta _0^{2\pi }𝑑\phi \underset{l=0}{\overset{\mathrm{}}{}}\sqrt{\frac{2l+1}{4\pi }}F_l(0)e^{\lambda [f_l1]t}Y_l^0(\theta )$$

$$(0,0,_0^\pi 𝑑\theta ^{}\mathrm{sin}\theta ^{}_0^{2\pi }𝑑\phi ^{}f(\theta ^{}\theta )\mathrm{cos}\theta ^{}\mathrm{cos}\theta )$$
(30)
Using the summation formula for Legendre polynomials (11) and the orthonormality of the spherical harmonics (12) we obtain

$$\stackrel{}{p}(t)=vF_1(0)e^{\lambda [f_11]t}\stackrel{}{k}.$$
(31)
For the $`z`$-stream initial condition, $`F(0,\mathrm{\Omega })=\frac{\delta (\mathrm{cos}\theta 1)}{2\pi }`$, $`F_1(0)=1`$, and $`\stackrel{}{p}(t)=ve^{\lambda [f_11]t}\stackrel{}{k}`$.

From (31),

$$\stackrel{}{x}(t)=vF_1(0)\frac{(e^{\lambda [f_11]t}1)}{\lambda [f_11]}\stackrel{}{k}.$$
(32)
and

$$|\stackrel{}{x}(\mathrm{})|=\frac{vF_1(0)}{\lambda [1f_1]}\stackrel{}{k}.$$
(33)
For the $`z`$-stream case we obtain a formula similar to the 2-dimensional one:

$$|\stackrel{}{x}(\mathrm{})|=\frac{v}{\lambda [1f_1]}\stackrel{}{k}.$$
(34)
## 4 Isotropic scattering

When the scattering is isotropic there is a direct way of obtaining the $`x`$-stream probability density ($`\rho _x`$) in terms of the probability density for isotropic initial conditions ($`\rho `$):

$$\rho _x(t,\stackrel{}{x})=e^{\lambda t}\delta (xvt)\delta (y)+\lambda _0^t𝑑t^{}e^{\lambda t^{}}\rho (tt^{},|\stackrel{}{x}vt^{}\stackrel{}{i}|).$$
(35)
The rationale for the above equation is that at each collision the stream gives rise to an isotropic probability cloud.

For the 2-dimensional case $`\rho `$ is known :

$$\rho (t,\stackrel{}{x})=e^{\lambda t}\left[\frac{\delta (rvt)}{2\pi r}+\frac{\lambda }{2\pi v\sqrt{(vt)^2r^2}}\mathrm{exp}\left(\frac{\lambda }{v}\sqrt{(vt)^2r^2}\right)H(vtr)\right],$$
(36)
where $`H`$ is the Heaviside function. The integral in (35) can be done analitically (see the Appendix). The result for $`\rho _x`$ is:

$$\rho _x(t,\stackrel{}{x})=e^{\lambda t}\left[\delta (xvt)\delta (y)+\frac{\lambda }{2\pi v(vtx)}\mathrm{exp}\left(\frac{\lambda }{v}\sqrt{(vt)^2r^2}\right)H(vtr)\right].$$
(37)
$`\rho _x(t,\stackrel{}{x})`$ is, in 2 dimensions, the probability density for finding the particle in $`\stackrel{}{x}`$ at time $`t`$ provided that at time $`t=0`$ it left the origin moving to the right. As long as the particle doesn’t collide its $`\rho _x(t,\stackrel{}{x})`$ will be $`\delta (vtx)\delta (y)`$, where $`v`$ is the speed of the particle. Since the particle does suffer, on the average, $`\lambda `$ collisions per unit time, the density of probability contributed by the possibility of no scattering is $`e^{\lambda t}\delta (vtx)\delta (y)`$. Once the particle is scattered at a time $`t`$’ and to a direction $`\phi ^{}`$ \- which will happen with a probability density $`\lambda f(\phi ^{})e^{\lambda t^{}}`$ \- its contribution to $`\rho _x(t,\stackrel{}{x})`$ will be $`\rho _x(tt^{},R_\phi ^{}(\stackrel{}{x}vt^{}\stackrel{}{i}))`$, where $`R_\phi ^{}`$ is the rotation that takes the $`\phi ^{}`$ direction back to the positive $`x`$-axis direction and $`\stackrel{}{i}`$ is the unit vector that points to the right. All this paragraph is the following integral equation:

$$\rho _x(t,\stackrel{}{x})=e^{\lambda t}\delta (vtx)\delta (y)+\lambda _0^t𝑑t^{}e^{\lambda t^{}}_0^{2\pi }𝑑\phi ^{}f(\phi ^{})\rho _x(tt^{},R_\phi ^{}(\stackrel{}{x}vt^{}\stackrel{}{i}))$$
(38)
The Fourier transform that we use for continuous variable is:

$$_{\mathrm{}}^{\mathrm{}}𝑑x\mathrm{function}(x)e^{i2\pi \nu x},$$
(39)
and the same transformation except for the sign of the exponent for its inverse. The Fourier transform of the integral equation is:

$$\stackrel{~}{\rho _x}(t,\stackrel{}{\nu })=e^{\lambda t}e^{i2\pi (\stackrel{}{\nu })_xvt}+\lambda _0^t𝑑t^{}_0^{2\pi }𝑑\phi ^{}f(\phi ^{})e^{\lambda t^{}}e^{i2\pi (\stackrel{}{\nu })_xvt^{}}\stackrel{~}{\rho _x}(tt^{},R_\phi ^{}\stackrel{}{\nu }),$$
(40)
where $`(\stackrel{}{\nu })_x`$ stands for the $`x`$ component of $`\stackrel{}{\nu }`$. We now take the Laplace transform (denoted by a hat) and obtain:

$$(\omega +\lambda +i2\pi (\stackrel{}{\nu })_xv)\widehat{\stackrel{~}{\rho _x}}(\omega ,\stackrel{}{\nu })=1+\lambda _0^{2\pi }𝑑\phi ^{}f(\phi ^{})\widehat{\stackrel{~}{\rho _x}}(\omega ,R_\phi ^{}\stackrel{}{\nu }),$$
(41)
where $`\omega `$ is the Laplace conjugate variable to time.

Consider the case of isotropic scattering ($`f(\phi ^{})=\frac{1}{2\pi }`$). Then the integral in the preceding equation is the Fourier-Laplace transform of the solution to the problem of isotropic scattering with isotropic initial conditions. For this case the following integral equation can be written for its probability density $`\rho `$ :

$$\rho (t,\stackrel{}{x})=\eta (t,\stackrel{}{x})+\lambda _0^t𝑑t^{}𝑑\stackrel{}{x}^{}\eta (t^{},\stackrel{}{x}^{})\rho (tt^{},\stackrel{}{x}\stackrel{}{x}^{})$$
(42)
where

$$\eta (t,\stackrel{}{x})=\frac{e^{\lambda t}}{2\pi |\stackrel{}{x}|}\delta (|\stackrel{}{x}|vt).$$
(43)
Using the Fourier and Laplace convolution theorems one can solve for the Fourier-Laplace transform of $`\rho `$ :

$$\widehat{\stackrel{~}{\rho }}=\frac{\widehat{\stackrel{~}{\eta }}}{1\lambda \widehat{\stackrel{~}{\eta }}}$$
(44)
in terms of the Fourier-Laplace transform of $`\eta `$:

$$\widehat{\stackrel{~}{\eta }}(\omega ,\stackrel{}{\nu })=\frac{1}{\sqrt{(\lambda +\omega )^2+(2\pi v\nu )^2}}$$
(45)
Substitution into (24) yields:

$$\widehat{\stackrel{~}{\rho _x}}(\omega ,\stackrel{}{\nu })=\frac{1}{(\omega +\lambda +i2\pi (\stackrel{}{\nu })_xv)(1\lambda \widehat{\stackrel{~}{\eta }}(\omega ,\stackrel{}{\nu }))}=$$

$$\frac{\sqrt{(\lambda +\omega )^2+(2\pi v\nu )^2}}{(\omega +\lambda +i2\pi (\stackrel{}{\nu })_xv)\left(\sqrt{(\lambda +\omega )^2+(2\pi v\nu )^2}\lambda \right)},$$
(46)
which is the Fourier-Laplace transform of (37).

In 3 dimensions the analogous procedure starts from the equation:

$$\rho _z(t,\stackrel{}{x})=e^{\lambda t}\delta (x)\delta (y)\delta (vtz)+\lambda _0^t𝑑t^{}e^{\lambda t^{}}𝑑\mathrm{\Omega }^{}f(\theta ^{})\rho _z(tt^{},R_\mathrm{\Omega }^{}(\stackrel{}{x}vt^{}\stackrel{}{k})).$$
(47)
We are assuming here that the scattering is independent of the azimuthal angle $`\phi `$. $`R_\mathrm{\Omega }^{}`$ is the rotation which maps the $`\mathrm{\Omega }^{}`$ direction into the z-direction and has its axis in the $`xy`$-plane. The scattering cross section is normalized as follows:

$$𝑑\mathrm{\Omega }f(\theta )=2\pi _0^\pi 𝑑\theta \mathrm{sin}\theta f(\theta )=1.$$
(48)
The Fourier-Laplace transform of the integral equation is now:

$$(\omega +\lambda +i2\pi (\stackrel{}{\nu })_zv)\widehat{\stackrel{~}{\rho _z}}(\omega ,\stackrel{}{\nu })=1+\lambda 𝑑\mathrm{\Omega }^{}f(\theta ^{})\widehat{\stackrel{~}{\rho _z}}(\omega ,R_\mathrm{\Omega }^{}\stackrel{}{\nu }),$$
(49)
Consider the case of isotropic scattering ($`f(\theta )=\frac{1}{4\pi }`$). Then the integral in the preceding equation is the Fourier-Laplace transform of the solution to the problem of isotropic scattering with isotropic initial conditions. For the this case the integral equation (42) holds and its solution is again given by (44). The difference is that $`\eta `$ is now:

$$\eta (t,\stackrel{}{r})=\frac{e^{\lambda t}}{4\pi r^2}\delta (rvt),$$
(50)
and $`\widehat{\stackrel{~}{\eta }}(\omega ,\stackrel{}{\nu })`$ is:

$$\widehat{\stackrel{~}{\eta }}(\omega ,\stackrel{}{\nu })=\frac{1}{2\pi v\nu }\mathrm{arc}\mathrm{tg}\frac{2\pi v\nu }{\lambda +\omega }.$$
(51)
Finally, the Fourier-Laplace transform of $`\widehat{\stackrel{~}{\rho _z}}(\omega ,\stackrel{}{\nu })`$ in 3 dimensions for the case of isotropic scattering is:

$$\widehat{\stackrel{~}{\rho _z}}(\omega ,\stackrel{}{\nu })=\frac{1}{(\omega +\lambda +i2\pi (\stackrel{}{\nu })_xv)}\frac{2\pi v\nu }{2\pi v\nu \lambda \mathrm{arc}\mathrm{tg}\frac{2\pi v\nu }{\lambda +\omega }}.$$
(52)
## 5 Anisotropic scattering

The form of the integral in (41) is similar to a convolution. This suggests the following procedure to factorize the unknown function out of the integral. First, write $`\stackrel{}{\nu }`$ as $`R_\phi \stackrel{}{\nu }`$:

$$(\omega +\lambda +i2\pi (R_\phi \stackrel{}{\nu })_xv)\widehat{\stackrel{~}{\rho _x}}(\omega ,R_\phi \stackrel{}{\nu })=1+\lambda _0^{2\pi }𝑑\phi ^{}f(\phi ^{})\widehat{\stackrel{~}{\rho _x}}(\omega ,R_{\phi \phi ^{}}\stackrel{}{\nu }).$$
(53)
Define now the following integral transform for a function $`g`$ of $`\stackrel{}{x}`$ ($`\stackrel{}{x}`$ is a 2-dimensional vector):

$$\overline{g}(\stackrel{}{x},\overline{\phi })_0^{2\pi }𝑑\phi g(R_\phi \stackrel{}{x})e^{i\overline{\phi }\phi },\overline{\phi }ϵZ.$$
(54)
From Fourier series we know that this transform can be inverted:

$$g(R_\phi \stackrel{}{x})=\frac{1}{2\pi }\underset{\overline{\phi }ϵZ}{}\overline{g}(\stackrel{}{x},\overline{\phi })e^{i\overline{\phi }\phi }.$$
(55)
It has the following property:

$$\overline{g}(R_\alpha \stackrel{}{x},\overline{\phi })=e^{i\alpha \overline{\phi }}\overline{g}(\stackrel{}{x},\overline{\phi })$$
(56)
Thus the transform that we have defined is, up to a phase factor, the Fourier transform of the restriction of the function to a circumference about the origin.

Application of $`_0^{2\pi }𝑑\phi e^{i\overline{\phi }\phi }`$ to both sides of the integral equation (53) yields:

$$2\pi \delta _{0,\overline{\phi }}+\lambda \stackrel{~}{f}(\overline{\phi })\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,\stackrel{}{\nu },\overline{\phi })=(\omega +\lambda )\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,\stackrel{}{\nu },\overline{\phi })+i2\pi v_0^{2\pi }𝑑\phi \widehat{\stackrel{~}{\rho _x}}(\omega ,R_\phi \stackrel{}{\nu })(R_\phi \stackrel{}{\nu })_xe^{i\overline{\phi }\phi }.$$
(57)
Let $`\phi _\stackrel{}{\nu }`$ denote the angle between $`\stackrel{}{\nu }`$ and the $`x`$ axis. Then the integral in the above equation can be developped as follows:

$$_0^{2\pi }𝑑\phi \widehat{\stackrel{~}{\rho _x}}(\omega ,R_\phi \stackrel{}{\nu })(R_\phi \stackrel{}{\nu })_xe^{i\overline{\phi }\phi }=$$

$$_0^{2\pi }𝑑\phi \widehat{\stackrel{~}{\rho _x}}(\omega ,R_\phi \stackrel{}{\nu })|\stackrel{}{\nu }|\frac{e^{i(\phi +\phi _\stackrel{}{\nu })}+e^{i(\phi +\phi _\stackrel{}{\nu })}}{2}e^{i\overline{\phi }\phi }=$$

$$\frac{|\stackrel{}{\nu }|}{2}\left[e^{i\phi _\stackrel{}{\nu }}\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,\stackrel{}{\nu },\overline{\phi }1)+e^{i\phi _\stackrel{}{\nu }}\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,\stackrel{}{\nu },\overline{\phi }+1)\right].$$
(58)
Substitution of this result back into (57) yields a set of recursive equations for the transform $`\overline{\widehat{\stackrel{~}{\rho }}}_x`$ of the probability density:

$$2\pi \delta _{0,\overline{\phi }}+[\lambda (\stackrel{~}{f}(\overline{\phi })1)\omega ]\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,\stackrel{}{\nu },\overline{\phi })=$$

$$i\pi v|\stackrel{}{\nu }|\left[e^{i\phi _\stackrel{}{\nu }}\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,\stackrel{}{\nu },\overline{\phi }1)+e^{i\phi _\stackrel{}{\nu }}\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,\stackrel{}{\nu },\overline{\phi }+1)\right],\overline{\phi }ϵZ.$$
(59)
If the scattering has left-right symmetry ($`f(\phi )=f(2\pi \phi )`$), then,

$$\stackrel{~}{f}(\overline{\phi })=\stackrel{~}{f}(\overline{\phi })$$
(60)
For $`\overline{\widehat{\stackrel{~}{\rho }}}_x`$ this symmetry implies:

$$\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,|\stackrel{}{\nu }|\stackrel{}{i},\overline{\phi })=\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,|\stackrel{}{\nu }|\stackrel{}{i},\overline{\phi }).$$
(61)
From this and (56):

$$\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,\stackrel{}{\nu },\pm \overline{\phi })=e^{\pm i\overline{\phi }\phi _\stackrel{}{\nu }}\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,|\stackrel{}{\nu }|\stackrel{}{i},\overline{\phi }).$$
(62)
Now the set of equations (59) can be simplified to:

$$2\pi \omega \overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,|\stackrel{}{\nu }|,0)=i2\pi v|\stackrel{}{\nu }|\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,|\stackrel{}{\nu }|\stackrel{}{i},1)$$
(63)
and

$$[\lambda (\stackrel{~}{f}(\overline{\phi })1)\omega ]\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,|\stackrel{}{\nu }|\stackrel{}{i},\overline{\phi })=i\pi v|\stackrel{}{\nu }|\left[\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,|\stackrel{}{\nu }|\stackrel{}{i},\overline{\phi }1)+\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,|\stackrel{}{\nu }|\stackrel{}{i},\overline{\phi }+1)\right],$$

$$\mathrm{for}\overline{\phi }=1,2,\mathrm{},$$
(64)
or (bearing (61) in mind and remembering that $`\stackrel{~}{f}(0)=1`$) simply to

$$[\lambda (1\stackrel{~}{f}(\overline{\phi }))+\omega ]\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,|\stackrel{}{\nu }|\stackrel{}{i},\overline{\phi })+i\pi v|\stackrel{}{\nu }|\left[\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,|\stackrel{}{\nu }|\stackrel{}{i},\overline{\phi }1)+\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,|\stackrel{}{\nu }|\stackrel{}{i},\overline{\phi }+1)\right]=2\pi \delta _{0,\overline{\phi }},$$

$$\mathrm{for}\overline{\phi }=0,1,2,\mathrm{}.$$
(65)
This more compact form brings out the matricial structure of the equations:

$$\left(\begin{array}{ccccc}& & & & \\ & i\pi v|\stackrel{}{\nu }|& [\lambda (1\stackrel{~}{f}(2))+\omega ]& i\pi v|\stackrel{}{\nu }|& 0\\ & 0& i\pi v|\stackrel{}{\nu }|& [\lambda (1\stackrel{~}{f}(1))+\omega ]& i\pi v|\stackrel{}{\nu }|\\ & 0& 0& 2i\pi v|\stackrel{}{\nu }|& [\lambda (1\stackrel{~}{f}(0))+\omega ]\end{array}\right)\left(\begin{array}{c}\\ \overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,|\stackrel{}{\nu }|\stackrel{}{i},2)\\ \overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,|\stackrel{}{\nu }|\stackrel{}{i},1)\\ \overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,|\stackrel{}{\nu }|\stackrel{}{i},0)\end{array}\right)$$

$$=\left(\begin{array}{c}\\ 0\\ 0\\ 2\pi \end{array}\right)$$
(66)
From (55), (61) and (62) the Fourier-Laplace transform of the probability density is:

$$\widehat{\stackrel{~}{\rho _x}}(\omega ,\stackrel{}{\nu })=\frac{1}{\pi }\left(\frac{\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,|\stackrel{}{\nu }|\stackrel{}{i},0)}{2}+\underset{\overline{\phi }=1,2,..}{}\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,|\stackrel{}{\nu }|\stackrel{}{i},\overline{\phi })\mathrm{cos}(\overline{\phi }\phi _\stackrel{}{\nu })\right)$$
(67)
Thus by approximating equation (66) by a finite dimensional truncation we can find the Fourier-Laplace transform of the probability density with arbitrary precision. There are recursive formulae that allow to invert finite tridiagonal matrices at a low computational cost .

We assume that at $`t=0`$ all particles leave from the same point and have the same speed. If this is not the case, summation over speeds and positions is necessary. Let $`I(\phi )=\frac{1}{2\pi }_{\overline{\phi }ϵZ}\stackrel{~}{I}(\overline{\phi })e^{i\overline{\phi }\phi }`$ be the initial angular distribution of velocities. Let $`\overline{\widehat{\stackrel{~}{\rho }}}_\phi (\omega ,\stackrel{}{\nu },\overline{\phi })`$ be the function corresponding to $`\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,\stackrel{}{\nu },\overline{\phi })`$ when the initial direction of the particle is $`\phi `$. $`R_\phi `$ is the adjoint operator of $`R_\phi `$, therefore:

$$\overline{\widehat{\stackrel{~}{\rho }}}_\phi (\omega ,\stackrel{}{\nu },\overline{\phi })=\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,R_\phi \stackrel{}{\nu },\overline{\phi }).$$
(68)
Then, using standard properties of Fourier sums and property (56) of the $`\overline{}`$ transform,

$$_0^{2\pi }𝑑\phi I(\phi )\overline{\widehat{\stackrel{~}{\rho }}}_\phi (\omega ,\stackrel{}{\nu },\overline{\phi })=\stackrel{~}{I}(\overline{\phi })\overline{\widehat{\stackrel{~}{\rho }}}_x(\omega ,\stackrel{}{\nu },\overline{\phi }).$$
(69)
In 3 dimensions we can proceed analogously. The starting equation is now (49). In it we write $`\stackrel{}{\nu }`$ as $`R_\mathrm{\Omega }\stackrel{}{\nu }`$.

$$(\omega +\lambda +i2\pi (R_\mathrm{\Omega }\stackrel{}{\nu })_zv)\widehat{\stackrel{~}{\rho _z}}(\omega ,R_\mathrm{\Omega }\stackrel{}{\nu })=1+\lambda 𝑑\mathrm{\Omega }^{}f(\theta ^{})\widehat{\stackrel{~}{\rho _z}}(\omega ,R_{\mathrm{\Omega }\mathrm{\Omega }^{}}\stackrel{}{\nu }),$$
(70)
3-dimensional rotations can not be added commutatively, but $`R_{\mathrm{\Omega }\mathrm{\Omega }^{}}`$ is well defined. It is the rotation that maps the $`\mathrm{\Omega }^{}`$ direction into the $`\mathrm{\Omega }`$ direction and has its axis perpendicular to the plane defined by these directions and the origin.

Define now the following integral transform for a function $`g`$ of $`\stackrel{}{x}`$ ($`\stackrel{}{x}`$ is a 3-dimensional vector):

$$\overline{g}(\stackrel{}{x},l)𝑑\mathrm{\Omega }g(R_\mathrm{\Omega }\stackrel{}{x})Y_l^0(\mathrm{\Omega }),lϵN.$$
(71)
If $`g`$ does not depend on the azimuthal angle $`\phi `$ the following 2 properties hold:

a)

$$\overline{g}(R_\mathrm{\Omega }\stackrel{}{x},l)=\overline{g}(\stackrel{}{x},l)\sqrt{\frac{4\pi }{2l+1}}Y_l^0(\mathrm{\Omega })$$
(72)
and b) the above transform is one-to-one and its inverse is:

$$g(\stackrel{}{x})=\underset{lϵN}{}\overline{g}(\stackrel{}{x},l)Y_l^0(0,0)=\underset{lϵN}{}\sqrt{\frac{2l+1}{4\pi }}\overline{g}(\stackrel{}{x},l),$$
(73)
where $`(0,0)`$ means $`(\theta ,\phi )=(0,0)`$. Similarly to (10) we define:

$$f_l=\sqrt{\frac{4\pi }{2l+1}}𝑑\mathrm{\Omega }f(\mathrm{\Omega })Y_l^0(\mathrm{\Omega }).$$
(74)
Application of $`𝑑\mathrm{\Omega }Y_l^0(\mathrm{\Omega })`$ to both sides of the integral equation (70) and use of the convolution theorem (13) yields:

$$(\omega +\lambda )\overline{\widehat{\stackrel{~}{\rho }}}_z(\omega ,\stackrel{}{\nu },l)+i2\pi v𝑑\mathrm{\Omega }\widehat{\stackrel{~}{\rho _z}}(\omega ,R_\mathrm{\Omega }\stackrel{}{\nu })(R_\mathrm{\Omega }\stackrel{}{\nu })_zY_l^0(\mathrm{\Omega })=\sqrt{4\pi }\delta _{0,l}+\lambda f_l\overline{\widehat{\stackrel{~}{\rho }}}_z(\omega ,\stackrel{}{\nu },l).$$
(75)
If $`\mathrm{\Omega }(\theta ,\phi )`$ and the polar coordinates of $`\stackrel{}{\nu }`$ are $`(\theta _\stackrel{}{\nu },\phi _\stackrel{}{\nu })`$, then geometry yields:

$$(R_\mathrm{\Omega }\stackrel{}{\nu })_z=|\stackrel{}{\nu }|(\mathrm{cos}\theta \mathrm{cos}\theta _\stackrel{}{\nu }\mathrm{sin}\theta \mathrm{sin}\theta _\stackrel{}{\nu }\mathrm{cos}(\phi \phi _\stackrel{}{\nu })).$$
(76)
Due to the assumed cilindrical symmetry of the scattering cross section we can take $`\phi _\stackrel{}{\nu }=0`$. Then, writing the trigonometric functions as linear combinations of spherical harmonics, we obtain:

$$(R_\phi \stackrel{}{\nu })_z=\sqrt{\frac{4\pi }{3}}|\stackrel{}{\nu }|\left[Y_1^0(\mathrm{\Omega })\mathrm{cos}\theta _\stackrel{}{\nu }\left(Y_1^1(\mathrm{\Omega })Y_1^1(\mathrm{\Omega })\right)\frac{\mathrm{sin}\theta _\stackrel{}{\nu }}{\sqrt{2}}\right].$$
(77)
When we substitute this expression back into (75) we are going to have products of spherical harmonics. By means of Clebsch-Gordan coefficients  the products $`Y_1^{\pm 1}(\mathrm{\Omega })Y_l^0(\mathrm{\Omega })`$ can be written as linear combinations of $`Y_{l\pm 1}^{\pm 1}(\mathrm{\Omega })`$. But it easy to check that for a cilindrically symmetric function $`g(\stackrel{}{x}),𝑑\mathrm{\Omega }g(R_\mathrm{\Omega }\stackrel{}{x})Y_l^m(\mathrm{\Omega })=0`$ if $`m0`$. Therefore only the product $`\sqrt{\frac{4\pi }{3}}Y_l^0(\mathrm{\Omega })Y_1^0(\mathrm{\Omega })`$ remains. Using the recurrence relations between Legendre polynomials  this product can be expanded as follows:

$$\sqrt{\frac{4\pi }{3}}Y_l^0(\mathrm{\Omega })Y_1^0(\mathrm{\Omega })=\frac{l+1}{\sqrt{(2l+1)(2l+3)}}Y_{l+1}^0(\mathrm{\Omega })+\frac{l}{\sqrt{(2l+1)(2l1)}}Y_{l1}^0(\mathrm{\Omega }).$$
(78)
Substituting this result back into (75) we obtain the following set of recursive equations for $`\overline{\widehat{\stackrel{~}{\rho }}}_z`$:

$$[\lambda (1f_l)+\omega ]\overline{\widehat{\stackrel{~}{\rho }}}_z(\omega ,\stackrel{}{\nu },l)+\frac{i2\pi v|\stackrel{}{\nu }|(l+1)}{\sqrt{(2l+1)(2l+3)}}\mathrm{cos}\theta _\stackrel{}{\nu }\overline{\widehat{\stackrel{~}{\rho }}}_z(\omega ,\stackrel{}{\nu },l+1)+\frac{i2\pi v|\stackrel{}{\nu }|l}{\sqrt{(2l+1)(2l1)}}\mathrm{cos}\theta _\stackrel{}{\nu }\overline{\widehat{\stackrel{~}{\rho }}}_z(\omega ,\stackrel{}{\nu },l1)$$

$$=\sqrt{4\pi }\delta _{0,l},$$

$$\mathrm{for}l=0,1,2,\mathrm{}.$$
(79)
In matricial form:

$$\left(\begin{array}{ccccc}& & & & \\ & \frac{i2\pi v|\stackrel{}{\nu }|3}{\sqrt{35}}\mathrm{cos}\theta _\stackrel{}{\nu }& [\lambda (1f_2)+\omega ]& \frac{i2\pi v|\stackrel{}{\nu }|2}{\sqrt{15}}\mathrm{cos}\theta _\stackrel{}{\nu }& 0\\ & 0& \frac{i2\pi v|\stackrel{}{\nu }|2}{\sqrt{15}}\mathrm{cos}\theta _\stackrel{}{\nu }& [\lambda (1f_1)+\omega ]& \frac{i2\pi v|\stackrel{}{\nu }|}{\sqrt{3}}\mathrm{cos}\theta _\stackrel{}{\nu }\\ & 0& 0& \frac{i2\pi v|\stackrel{}{\nu }|}{\sqrt{3}}\mathrm{cos}\theta _\stackrel{}{\nu }& [\lambda (1f_0)+\omega ]\end{array}\right)\left(\begin{array}{c}\\ \overline{\widehat{\stackrel{~}{\rho }}}_z(\omega ,\stackrel{}{\nu },2)\\ \overline{\widehat{\stackrel{~}{\rho }}}_z(\omega ,\stackrel{}{\nu },1)\\ \overline{\widehat{\stackrel{~}{\rho }}}_z(\omega ,\stackrel{}{\nu },0)\end{array}\right)$$

$$=\left(\begin{array}{c}\\ 0\\ 0\\ \sqrt{4\pi },\end{array}\right)$$
(80)
where, for symmetry’s sake, we have not simplified $`[\lambda (1f_0)+\omega ]`$ to $`\omega `$.

Truncation of the equations to a finite $`l`$ yields a set of equations that can be solved. As in the 2-dimensional case the matrix is tridiagonal and it can be inverted using recursive formulae . By means of the inversion formula (73) we can thus find the Fourier-Laplace transform of the probability density with arbitrary precision.

As in the 2-dimensional case we assume that at $`t=0`$ all particles leave from the same point and have the same speed. Let $`I(\mathrm{\Omega })=_{l,m}I_{l,m}Y_l^m(\mathrm{\Omega })`$ be the initial angular distribution of velocities. Let $`\overline{\widehat{\stackrel{~}{\rho }}}_\mathrm{\Omega }(\omega ,\stackrel{}{\nu },l)`$ be the function corresponding to $`\overline{\widehat{\stackrel{~}{\rho }}}_z(\omega ,\stackrel{}{\nu },l)`$ when the initial direction of the particle is $`\mathrm{\Omega }`$. $`R_\mathrm{\Omega }`$ is the adjoint operator of $`R_\mathrm{\Omega }`$, therefore:

$$\overline{\widehat{\stackrel{~}{\rho }}}_\mathrm{\Omega }(\omega ,\stackrel{}{\nu },l)=\overline{\widehat{\stackrel{~}{\rho }}}_z(\omega ,R_\mathrm{\Omega }\stackrel{}{\nu },l).$$
(81)
Then, using standard properties of spherical harmonics and property (72) of the $`\overline{}`$ transform,

$$𝑑\mathrm{\Omega }I(\mathrm{\Omega })\overline{\widehat{\stackrel{~}{\rho }}}_\mathrm{\Omega }(\omega ,\stackrel{}{\nu },l)=\sqrt{\frac{4\pi }{2l+1}}\overline{\widehat{\stackrel{~}{\rho }}}_z(\omega ,\stackrel{}{\nu },l)I_{l,0}.$$
(82)
This result does not depend on the value of the $`I_{l,m}`$’s for $`m0`$. As noted earlier $`𝑑\mathrm{\Omega }I(\mathrm{\Omega })\overline{\widehat{\stackrel{~}{\rho }}}_\mathrm{\Omega }(\omega ,\stackrel{}{\nu },l)`$ can be inverted by means of (73) to yield the Fourier-Laplace transform of $`\rho (t,\stackrel{}{x})`$ only when $`I`$ is cilidrically symmetric. Thus to apply formula (82) one has to decompose a general initial condition as a sum of cilindrically symmetric (with respect to different axis) components and then apply (82) to each component separately.

## Acknowledgments

Conversations with Andrés Hombría Maté are happily remembered.

Appendix

Substitution of (36) into (35) yields 3 terms, which are:

$$e^{\lambda t}\delta (xvt)\delta (y),$$
(83)
$$\lambda e^{\lambda t}_0^t𝑑t^{}\frac{\delta (|\stackrel{}{x}vt^{}\stackrel{}{i}|v(tt^{}))}{2\pi |\stackrel{}{x}vt^{}\stackrel{}{i}|}$$
(84)
and

$$\lambda ^2e^{\lambda t}_0^t𝑑t^{}\frac{\mathrm{exp}\left(\frac{\lambda }{v}\sqrt{(v(tt^{}))^2|\stackrel{}{x}vt^{}\stackrel{}{i}|^2}\right)}{2\pi v\sqrt{(v(tt^{}))^2|\stackrel{}{x}vt^{}\stackrel{}{i}|^2}}H(v(tt^{})|\stackrel{}{x}vt^{}\stackrel{}{i}|).$$
(85)
Since the integration variable is not the argument of the Dirac delta function in the second term, the following formula will be needed:

$$𝑑x\delta (u(x))=\underset{u(x)=0}{}\frac{1}{|u^{}(x)|}.$$
(86)
It is convenient to write the integral in the second term (84) as follows:

$$\frac{1}{2\pi v}_0^{vt}𝑑y\frac{\delta (|\stackrel{}{x}y\stackrel{}{i}|vt+y))}{2\pi |\stackrel{}{x}y\stackrel{}{i}|}.$$
(87)
The derivative that appears in the r.h.s. of (86) is then:

$$\frac{}{y}(|\stackrel{}{x}y\stackrel{}{i}|vt+y)=1+\frac{yx}{|\stackrel{}{x}y\stackrel{}{i}|}.$$
(88)
Using of the last 3 formulae the second term becomes

$$\frac{\lambda e^{\lambda t}}{2\pi v}\frac{1}{yx+|\stackrel{}{x}y\stackrel{}{i}|}$$
(89)
The argument of the Dirac delta function in the second term vanishes for

$$y=\frac{v^2t^2r^2}{2(vtx)}.$$
(90)
Substituting this in (88) yields

$$\frac{\lambda e^{\lambda t}}{2\pi v}\frac{1}{vtx}.$$
(91)
Defining

$$Av^2t^2r^2$$
(92)
and

$$B2v(xvt),$$
(93)
the third term (85) can be written as

$$\frac{\lambda ^2e^{\lambda t}}{2\pi v}_0^{\frac{v^2t^2r^2}{2(vtx)}}𝑑t^{}\underset{j=0}{\overset{\mathrm{}}{}}\frac{\lambda ^j(A+Bt^{})^{\frac{j1}{2}}}{j!v^j}=$$

$$\frac{\lambda e^{\lambda t}}{2\pi v(xvt)}\left[\mathrm{exp}\frac{\lambda \sqrt{A+Bt^{}}}{v}1\right]_0^{\frac{v^2t^2r^2}{2(vtx)}}=$$

$$\frac{\lambda e^{\lambda t}}{2\pi v(vtx)}\left[\mathrm{exp}\left(\frac{\lambda \sqrt{v^2t^2r^2}}{v}\right)1\right].$$
(94)
Collecting results we obtain (37).

## References

> A. V. Starkov, M. Noormohammadian, U. G. Oppel, ”A stochastic model and a variance-reduction Monte-Carlo method for the calculation of light transport”, Applied Physics B, Lasers and Optics, B 60, (1995) 335-340.
>
>  Milton Kerker, ”The scattering of light” (Academis Press, 1969).
>
>  R. F. Bonner and R. Nossal, Appl. Opt. 20 (1981) 2097.
>
>  A. P. Sheperd and P. A. Oberg, Laser-Doppler Blood Flowmetry (Kluwer, New York, 1989)
>
>  B. Chance et alii, Pro. Natl. Acad. Sci. USA 85 (1988) 4971.
>
>  J. M. Fernández-Varea et alii, Nucl. Instr. and Meth. in Phys. Res. B73 (1993), 447-473.
>
>  R. García-Pelayo, Physica A, 216, (1995) 299-315.
>
>  G. I. Bells and S. Glasstone, ”Nuclear Reactor Theory” (Van Nostrand Reinhold, 1970)
>
>  J. H. Ferziger et al., ”The theory of neutron slowing down in nuclear reactors” (Pergamon Press, New York, 1966).
>
>  Harry Hochstadt, ”The Functions of Mathematical Physics”, Dover (1986).
>
>  Chun Wa Wong, ”Introduction to Mathematical Physics”, Oxford (1991).
>
>  S. Goudsmit and J. L. Saunderson, Phys. Rev. 57 (1940) 24.
>
>  S. Goudsmit and J. L. Saunderson, Phys. Rev. 58 (1940) 36.
>
>  H. W. Lewis, Phys. Rev. 78 (1950) 526.
>
>  C. S. Johnson Jr., Don A. Gabriel, ”Laser light scattering”, Dover, (1995) p. 53.
>
>  M. Boguñá et alii, Physica A 230 (1996) 149-155.
>
>  J. Masoliver, J. M. Porrá and G. H. Weiss, Physica A 193 (1993) 469-482.
>
>  Y. Huang and W. F. McColl, J. Phys. A (1997) 7919-7933.
>
>  C. Cohen-Tannoudji, B. Diu and F. Laloë, ”Quantum Mechanics”, Wiley-Interscience (1977).