# Diffusion of particles moving with constant speed

## I Introduction

The propagation of light through a stochastic medium is traditionally described in the context of astrophysics by a Boltzmann transport equation for the specific intensity $`I(\stackrel{}{r},\stackrel{}{\mathrm{\Omega }},t)`$ in a heuristic Radiative transfer theory . However, since the general analytic solutions are unknown, one resorts to the diffusion approximation which can be shown to arise out of the Radiative transport equation in the limit of large length scales $`L>>l^{}`$, where $`l^{}`$ is the transport mean free path of light in the medium . Recently there has been considerable interest in the description of multiple light scattering at small length scales $`(Ll^{})`$ and small time scales ($`tt^{}`$ where $`t^{}`$ is the transport mean free time), both from the point of fundamental physics  and from the point of medical imaging, where the early arriving ‘snake’ photons are used to image through human tissues . It has been experimentally shown that the diffusion approximation fails to describe phenomena at distances of $`L<8l^{}`$ . Moreover the diffusion approximation which is strictly a Wiener process for the spatial co-ordinates of a particle is physically unrealistic. It holds in the limit of the mean free path $`l^{}0`$ and the speed of propagation $`c\mathrm{}`$ while keeping the diffusion coefficient $`D_o=\frac{cl^{}}{3}`$ constant. Thus the diffusion approximation neither accounts for a finite mean free path nor for a finite and constant speed of the particle which is charecteristic of light propagation in a stochastic medium. While approximately describing light as a particle the constancy of speed should be preserved at the very least. Hence it is of importance to develop better and alternative schemes to the diffusion approximation and also address the difficult question of the process of randomization of a directional beam in such media.   

For a particle moving with fixed speed in a one-dimensional disordered medium, it has be shown that the probability distribution function $`P(x,t)`$ for the displacement satisfies the telegrapher equation exactly. However, generalizations of the telegrapher equation to higher dimensions  have been shown not to yield better results than the diffusion approximation . Recently there have been a few attempts to overcome the shortcomings of the diffusion approximations and attack this problem using the concept of photon paths. In Ref., a Monte-Carlo approach was used to simulate photon paths and calculate their probabilities. An important advance was made in Refs., where the propagator for photons in highly forward scattering media was expressed as a Feynman path integral. However, this attempt has had only limited success in that it was possible to calculate the probability distribution subject to the constraint of constant photon speed only in the weaker (average) sense i.e., $`_0^t\left[(\frac{d\stackrel{}{r}}{dt})^2c^2\right]𝑑t=0`$. Moreover in addressing the backscattering from a semi-infinite medium  and reflection/transmission from a finite slab , the absorbing boundary conditions have not been rigorously implemented and it would be inappropriate to compare these to experimental data. It should be mentioned here that the Ornstein-Uhlenbeck (O-U) theory of Brownian motion  would also be able to incorporate the finiteness of the mean free path and a well defined root-mean-squared (rms) velocity but assuming, of course a distribution of speeds. This process has been compared with Monte-Carlo simulations  and used to explain the lowering of the effective diffusion coefficient measured in pulse transmission experiments through thin slabs . It can be shown that the finite r.m.s speed defined by the fluctuation-dissipation theorem for the O-U process is a stronger global constraint than the average speed constraint imposed in Ref..   

The next important step in describing these photon random walks with a constant speed was undertaken in Ref., where the authors describe this process as a non-Euclidean diffusion on the velocity sphere and intuitively put down a kind of a general Boltzmann equation for photons in a highly forward scattering medium. The solution to this equation was expressed as a path integral, which was then evaluated by a standard cumulant decomposition  truncated after the second cumulant. This yields a gaussian distribution similiar to the Ornstein-Uhlenbeck process. More recently, an explicit derivation of the Feyman path integral representation for the propagator of the radiative transfer equation has been given . Here it was again evaluated by truncating the cumulant expansion after the second term. This was justified by declaring that photons are massless and non-interacting. However, the imposition of the speed constraint would not allow this gaussian approximation.   

In this paper, we describe the light propagation in stochastic media as the motion of a kind of Brownian particle on which the fluctuating forces act only perpendicular to the direction of its velocity. This is effective in strictly and dynamically preserving the speed of the particle. This process is shown to correspond to a diffusion in the angular co-ordinate in the velocity space for a white noise disorder. Exact expressions for the moments of the space variables are presented and the second cumulant approximation is shown to yield a gaussian expression similiar to the traditional Ornstein-Uhlenbeck theory of Brownian motion. An expression is derived for the probability distribution for large force strengths which preserves the light cone. The exact Fokker-Planck equation for the probability distribution is derived from the stochastic Langevin equations for a white noise process. Numerical solutions of this equation are presented . It is shown that the probability distribution in infinite media is strongly forward peaked for short times and randomizes only at times of about $`8t^{}`$ to $`10t^{}`$. We have also solved numerically the equation for a semi-infinite geometry and obtained the persistence exponent of $`0.435i\pm 0.005`$ in 2-dimensions for this process. Solutions for a finite geometry are also given, showing that the effective diffusion coefficient as measured in a pulse transmission experiment through very thin slabs ($`Ll^{}`$) would be lowered. The effect of light amplification in the slab is examined briefly.   

## II The modified Ornstein Uhlenbeck process

Light scattering in a stochastic medium is treated as a probabilistic process where each scattering event only changes the direction of the photon. The wave nature and polarization effects are ignored and light is treated as a particle in a medium which exerts transverse fluctuating forces on the particle. It should be remarked here that while the actual disorder is maybe in space (quenched disorder), all current treatments including ours, are in terms of a Brownian motion (temporal disorder i.e., a stochastic process). This is a valid approximation for incoherent transport in the weak scattering limit ($`kl^{}>>1`$ where $`k=\frac{2\pi }{\lambda }`$, $`\lambda `$ being the wavelength of light in the medium. The equation for the motion of a randomly accelerated particle with the special condition that the random forces always act only perpendicular to the velocity can be written as,

$$\ddot{\stackrel{}{r}}=\dot{\stackrel{}{r}}\times \stackrel{}{f}(t)$$
(1)   
This we term as the modified Ornstein-Uhlenbeck process. We will consider two dimesions for simplicity, and write

$`\ddot{x}`$ $`=`$ $`f(t)\dot{y}`$ (2)   
$`\ddot{y}`$ $`=`$ $`f(t)\dot{x}`$ (3)   
where the force term $`f(t)`$ is a random funtion of time. We will assume a delta-correlated force with gaussian distribution. i.e.,

$`f(t)`$ $`=`$ $`0`$ (4)   
$`f(t)f(t^{})`$ $`=`$ $`\mathrm{\Gamma }\delta (tt^{})`$ (5)   
and all higher moments of $`f(t)`$ being zero. This makes our treatment most valid for a very dense collection of highly forward scattering weak anisotropic scatterers. This set of stochastic Langevin equations yield on integration a first constant of integration $`\dot{x}^2+\dot{y}^2=c^2`$, where $`c`$ is the constant speed. So we can choose $`\dot{x}=c\mathrm{cos}\theta (t)`$ and $`\dot{y}=c\mathrm{sin}\theta (t)`$ where $`\theta (t)`$ is some function of $`t`$. $`\theta `$ is recognised to be the angular co-ordinate in the velocity space. Substituting these expressions back into equation(2,3), we obtain $`\dot{\theta }=f(t)`$ or

$$\theta (t)\theta _0=_0^tf(t)𝑑t$$
(6)   
Hence $`\theta (t)`$ follows a Wiener process and we can write the probability for $`\theta (t)`$ as

$$P_t(\theta )=\left(\frac{1}{2\pi \mathrm{\Gamma }t}\right)^{\frac{1}{2}}\mathrm{exp}\left\{\frac{(\theta \theta _0)^2}{2\mathrm{\Gamma }t}\right\}$$
(7)   
This is the result for a diffusion in $`\theta `$ the angular co-ordinate in the velocity space, and we recognise this modified O-U process to be a random walk on the circle of radius $`c`$ in the velocity space. Constraining $`\theta `$ to the range $`[0,2\pi ]`$, we get the marginal probability distribution for $`\theta `$.

$$P_t(\theta )=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\left(\frac{1}{2\pi \mathrm{\Gamma }t}\right)^{\frac{1}{2}}\mathrm{exp}\left\{\frac{(\theta \theta _0+2n\pi )^2}{2\mathrm{\Gamma }t}\right\}$$
(8)   
The value of $`\theta _0`$ can be conveniently chosen to be zero.   

Now we will derive the probability distribution function in the phase space. Consider the system of Stochastic Langevin equations

$`\dot{x}=u`$ (9)   
$`\dot{y}=v`$ (10)   
$`\dot{u}=f(t)v`$ (11)   
$`\dot{v}=f(t)u`$ (12)   
Let $`\mathrm{\Pi }(x,y,u,v)`$ be the phase space density of points for the given system and $`𝐔`$ be the vector $`(x,y,u,v)`$. Now, $`\mathrm{\Pi }`$ satisfies the stochastic Liouville equation.

$$\frac{\mathrm{\Pi }}{t}+_𝐔\left(\dot{𝐔}\mathrm{\Pi }\right)=0$$
(13)   
where $`_𝐔=(\frac{}{x},\frac{}{y},\frac{}{u},\frac{}{v})`$. Substituting for $`𝐔`$ and averaging over all possible configurations of disorder, by the van Kampen Lemma , the probability distribution $`P(x,y,u,v)=\mathrm{\Pi }(x,y,u,v)`$ and satisfies,

$$\frac{P}{t}+u\frac{P}{x}+v\frac{P}{y}v\frac{}{u}f(t)\mathrm{\Pi }+u\frac{}{v}f(t)\mathrm{\Pi }=0$$
(14)   
By the Novikov theorem  for a white noise process $`f(t)`$,

$`f(t)\mathrm{\Pi }[f(t)]`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }}{2}}{\displaystyle \frac{\delta \mathrm{\Pi }[f]}{\delta f(t)}}`$ (15)   
$`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }}{2}}\left(v{\displaystyle \frac{P}{u}}u{\displaystyle \frac{P}{v}}\right)`$ (16)   
Using the above, we obtain for P(x,y,u,v) the differential equation

$$\frac{P}{t}+u\frac{P}{x}+v\frac{P}{y}=\frac{\mathrm{\Gamma }}{2}\left(u\frac{}{v}v\frac{}{u}\right)^2P$$
(17)   
Now expressing u and v in terms of the angular co-ordinate $`\theta `$, we finally get

$$\frac{P}{t}+c\mathrm{cos}\theta \frac{P}{x}+c\mathrm{sin}\theta \frac{P}{y}=2\mathrm{\Gamma }\frac{^2P}{\theta ^2}$$
(18)   
This differential equation explicitly preserves the constancy of the speed of the photon. This Fokker-Planck equation is the same equation (in two dimensions) which was written down in Ref.. It is rigorously proved therein that this has a path integral solution and the two approaches are equivalent. It appears that this equation has solutions in terms of Mathieu functions. However we have not been able to analytically solve the equation.   

The moments of the displacements can however be calculated analytically. The displacements can be written in terms of $`\theta `$ as

$`xx_0`$ $`=`$ $`{\displaystyle _0^t}c\mathrm{cos}\theta d\theta `$ (19)   
$`yy_0`$ $`=`$ $`{\displaystyle _0^t}c\mathrm{sin}\theta d\theta `$ (20)   
Using these and a gaussian distribution for $`f(t)`$, we get

$`xx_0`$ $`=`$ $`{\displaystyle \frac{2c}{\mathrm{\Gamma }}}\mathrm{cos}\theta _0\left(1e^{\frac{\mathrm{\Gamma }t}{2}}\right)`$ (21)   
$`yy_0`$ $`=`$ $`{\displaystyle \frac{2c}{\mathrm{\Gamma }}}\mathrm{sin}\theta _0\left(1e^{\frac{\mathrm{\Gamma }t}{2}}\right)`$ (22)   
$`(xx_0)^2`$ $`=`$ $`c^2\left[{\displaystyle \frac{2t}{\mathrm{\Gamma }}}{\displaystyle \frac{2}{3}}\left({\displaystyle \frac{2}{\mathrm{\Gamma }}}\right)^2\left(1e^{\frac{\mathrm{\Gamma }t}{2}}\right){\displaystyle \frac{1}{12}}\left({\displaystyle \frac{2}{\mathrm{\Gamma }}}\right)^2\left(1e^{\mathrm{\Gamma }t}\right)\right]`$ (23)   
$`(yy_0)^2`$ $`=`$ $`c^2\left[{\displaystyle \frac{2t}{\mathrm{\Gamma }}}{\displaystyle \frac{4}{3}}\left({\displaystyle \frac{2}{\mathrm{\Gamma }}}\right)^2\left(1e^{\frac{\mathrm{\Gamma }t}{2}}\right)+{\displaystyle \frac{1}{12}}\left({\displaystyle \frac{2}{\mathrm{\Gamma }}}\right)^2\left(1e^{\mathrm{\Gamma }t}\right)\right]`$ (24)   
$`(xx_0)(yy_0)`$ $`=`$ $`0`$ (25)   
This reproduces the result of the traditional Ornstein-Uhlenbeck process in that the first moment saturates at a mean free path $`l^{}`$ and the second moment increases linearly with time at long times ($`\mathrm{\Gamma }t/21`$). For short times ($`\mathrm{\Gamma }t/21`$), the longitudinal spread $`\mathrm{\Delta }x^2t^2`$ and the lateral spread $`\mathrm{\Delta }y^2t^3`$ which are considerably slower than the diffusive linear behaviour. From these relations, we identify the mean free time $`t^{}`$ to be $`\frac{2}{\mathrm{\Gamma }}`$ and the transport mean free path $`l^{}=ct^{}`$. The diffusion coefficient is identified as the coefficient of the linear term of the second moment i.e., $`\frac{c^2}{\mathrm{\Gamma }}`$   

It is of interest to note that an analytic expression for moments of all orders for the displacements can be obtained. This expression is given in the Appendix. The marginal probability distribution function $`P(x,y,t;x_0,y_0,0)`$ can be written in terms of a cumulant expansion (See the Appendix). Truncation of the cumulant series after the second term yields the result of Ref. for the probability distribution.

$`P(x,y,t;x_0,y_0,0)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi det(M)}}\mathrm{exp}\left\{{\displaystyle \frac{M_{ij}^1}{2}}(\stackrel{}{r}\stackrel{}{r_0}\stackrel{}{a})_i(\stackrel{}{r}\stackrel{}{r_0}\stackrel{}{a})_j\right\}`$ (26)   
$`\stackrel{}{a}`$ $`=`$ $`{\displaystyle \frac{2c}{\mathrm{\Gamma }}}(1e^{\mathrm{\Gamma }t/2})(\mathrm{cos}\theta _0,\mathrm{sin}\theta _0)`$ (27)   
$`M_{ij}`$ $`=`$ $`(\stackrel{}{r}\stackrel{}{r_0})_i(\stackrel{}{r}\stackrel{}{r_0})_j(\stackrel{}{r}\stackrel{}{r_0})_i(\stackrel{}{r}\stackrel{}{r_0})_j`$ (28)   
The distribution is gaussian in this approximation and similiar to the distribution for the traditional O-U process . Thus it does not exactly preserve the light cone and would appear to constrain the speed only in an average sense. Higher cumulants would be required to describe this feature of fixed speed.   

An approximate solution which preserves the light cone can be obtained under the assumption that $`\theta `$ is completely randomized in time $`t^{}`$, so that $`\theta `$ has a uniform distribution over $`[0,2\pi ]`$. This can be justified in the limit of large force strength ($`\mathrm{\Gamma }`$), when the scattering events change the momentum by a large amount. Now the time can be discretized on this time-scale and the probability distribution can be written as

$$P(x,t;x_0,0)=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑\omega e^{i\omega (xx_0)}\mathrm{exp}[i\omega c\underset{j=1}{\overset{n}{}}\mathrm{cos}\theta _jt^{}]$$
(29)   
where we have used that at $`t=0`$, the angle $`\theta `$ was uniformly distributed. To evaluate the average, we will use the fact that in this approximation, each $`\theta _j`$ is independent of all others, giving

$$\mathrm{exp}[i\omega c\underset{j=1}{\overset{n}{}}\mathrm{cos}\theta _jt^{}]=\left[J_0(\omega ct^{})\right]^n$$
(30)   
where $`J_0`$ is the ordinary Bessel’s function of order zero. Using $`n=t/t^{}`$, we have

$$P(x,t;x_0,0)=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑\omega e^{i\omega (xx_0)}\left[J_0(\omega ct^{})\right]^{t/t^{}}$$
(31)   
Using the fact that the Fourier transform of $`J_0(\omega ct^{})`$ is zero for $`\left|xx_0\right|>ct^{}`$ and the fact that $`P(x,t,x_0,0)`$ is an $`n^{th}`$ convolution of $`J_0(\omega ct^{})`$, it is seen that $`P(x,t,x_0,0)`$ is zero for $`\left|xx_0\right|>nct^{}=ct`$. Thus the light cone is preserved. In fig-1, we plot the $`P(x,t,x_0,0)`$ obtained by numerical evaluation for different times. It is seen that for $`t=t^{}`$, the probability is accumulated at the light front, and all the curves show a cut-off at $`\left|xx_0\right|=ct`$. At long times, using the Laplace approximation, we have (for large $`n`$);

$`\left[J_0(\omega ct^{})\right]^n`$ $``$ $`\mathrm{exp}\left\{{\displaystyle \frac{c^2t_{}^{}{}_{}{}^{2}\omega ^2n}{4}}\right\}`$ (32)   
$`P(x,t;x_0,0)`$ $``$ $`\left({\displaystyle \frac{1}{\pi c^2t^{}t}}\right)^{1/2}\mathrm{exp}\left\{{\displaystyle \frac{(xx_0)^2}{c^2t^{}t}}\right\}`$ (33)   
Thus we recover the diffusion limit at long times.   

## III Numerical Solutions and Results

In this section , numerical solutions for the differential equation (18) are presented. The particle is released in the $`xy`$ plane at the origin (generally) along an initial direction $`\theta _0`$. Here $`\theta `$ is the angle made by the velocity vector with the $`x`$-axis. Let us first further simplify by assuming invarience with respect to $`y`$ i.e., we have a line source along the $`y`$ axis. Then the derivative with respect to $`y`$ drops out and and we have a partial differential equation in three variables. This is essentially a parabolic equation with an advective term. To numerically propagate the probability distribution in time, we use an alternating direction implicit -explicit method  for $`x`$ and $`\theta `$. A local von Neumann stability analysis  shows that this differencing scheme is unconditionally stable. The initial condition is a $`\delta `$-function at $`x=0,\theta =0`$ which is approximated by a sharp gaussian for numerical purposes. For infinite media, the boundary condition $`P(x,t)=0`$ for $`\left|x\right|>ct`$ is used. For a semi-infinite medium $`\mathrm{}<x<L`$ with an absorbing boundary at $`x=L`$, the appropriate boundary condition is given by $`P(L,\theta ,t;x_0,\theta _0,0)=0`$ for $`\pi <\theta <\pi /2`$ and $`\pi /2<\theta <\pi `$, corresponding to no flux entering the medium from free space. Also, we can write the Fokker-Plank equation in the form of a continuity equation.

$`{\displaystyle \frac{P}{t}}+\stackrel{}{ȷ}=0`$ (34)   
$`=\widehat{e}_x{\displaystyle \frac{}{x}}+\widehat{e}_\theta {\displaystyle \frac{}{\theta }}`$ (35)   
$`\stackrel{}{ȷ}=\widehat{e}_x\mathrm{cos}\theta P+\widehat{e}_\theta \mathrm{\Gamma }{\displaystyle \frac{P}{\theta }}`$ (36)   
Since $`\mathrm{\Gamma }=0`$ outside the medium, we can conclude that the current density $`\stackrel{}{ȷ}`$ in the real ($`x`$) space is conserved across the boundary in the forward direction $`(\pi /2<\theta <\pi /2)`$ while the current density in the velocity ($`\theta `$) space is not conserved. This explains why the output flux at the boundary is proportional to the value of the probability distribution function at the boundary itself ( rather than the space derivative of the probability distribution $`(\frac{P}{x})`$ given by Fick’s law ) as observed in experiments . For a finite slab we use a similiar boundary condition at the other boundary.   

In Fig-2, we show the probability distributions in an infinite medium with the initial condition, $`P(x,\theta ,t=0)=\delta (x)\delta (\theta )`$. It is clearly seen that the probability distribution for times upto $`5t^{}`$ is peaked in the forward direction $`\theta 0`$ for $`x>0`$, with a tail in the backward direction $`(\theta \pm \pi )`$ at $`x<0`$. There is also a clear cut off at $`\left|x\right|=ct`$, which is prominently noticeable for positive $`x`$. The small amount of tailing arises from the finite width of the gaussian by which the $`\delta `$-function was approximated. One can also note that the probability distribution becomes almost flat along the $`\theta `$-axis only at times of about eight times the mean free time ($`8t^{}`$). In Fig-3, the first and second moments of the $`x`$ co-ordinate are shown. The solid lines show the analytical results of equation(21 and 23) and the symbols($``$) represent the results of the numerical solutions. Excellent agreement is found between them. In Fig-4, we show the marginal probability distribution for $`x`$ i.e., $`P(x,t;x_0,\theta _0,0)=_\pi ^\pi 𝑑\theta P(x,\theta ,t;x_0,\theta _0,0)`$. At short times ($`t3t^{}`$), there is a clear ballistic peak, separate from the more randomized tail. The probability distribution for these times is also clearly forward-peaked. One can also note that the probability distribution randomizes and becomes almost gaussian, centred at $`xl^{}`$ only at times $`t8t^{}`$. As noted above, this is also the time by which the angular coordinate $`\theta `$ randomizes. This is when the diffusion approximation becomes valid. This can be understood by noting that, by equation(7), the time required for $`P_t(\theta )`$ to attain an angular width of $`2\pi `$ is $`T`$ where $`T`$ is given by $`\mathrm{\Delta }\theta ^2=(2\pi )^22\mathrm{\Gamma }T`$. This yields (using $`\mathrm{\Gamma }/2=t^{}`$) a value of $`T=\pi ^2t^{}10t^{}`$ for the randomization time. Thus we now have a clear picture of the reason for this long known experimental fact . This forward peaked behaviour at short times also illustrates the deficiency of the second cumulant approximation where the probability distribution is a Gaussian and symmetric about the first moment. Higher cumulants are clearly required to describe these asymmetric features.   

The probability distribution functions for a semi-infinite medium are shown in Fig-5. Here the particle is released at the origin inside the random medium and the initial direction is towards the boundary (in this case at $`x=4l^{}`$) For times lesser than $`4t^{}`$, there is no difference in the probability distribution from the case of the infinite medium. This is because the wave front has not propagated upto the boundary and the effect of the boundary is not felt. This is to be contrasted with the diffusion approximation where the effect of the boundary is felt everywhere simultaneously and causality is violated. At long times the probability distributions attain a typical shape with a long tail at negative $`x`$ within the medium and a sharp cut-off at the boundary. In Fig-6, we show the marginal probability distribution for $`x`$ i.e., $`P(x,t;x_0,\theta _0,0)=_\pi ^\pi 𝑑\theta P(x,\theta ,t;x_0,\theta _0,0)`$. The value of $`P(x,t;x_0,\theta _0,0)`$ is finite at the boundary and zero outside. As seen in Fig-6b, if the points near the boundary are linearly extrapolated outside the boundary, they all roughly cross the x-axis at about $`0.7l^{}`$ which is the value of the extrapolation length used in the diffusion approximation . In Fig-7, the surviving probability inside the medium $`P_s=𝑑x𝑑\theta P(x,\theta ,t;x_0,\theta _0,0)`$ is plotted with time. For long times, this quantity should scale as $`t^\vartheta `$ where $`\vartheta `$ is the persistence exponent for this process . We have performed these calculations for several source-boundary distances and obtained a value of $`0.435\pm 0.005`$ as the persistence exponent for this process in two dimensions.   

Finally we present solutions for a finite slab with absorbing boundaries at $`x=\pm L`$. The particle is released from the origin at $`t=0`$ along the positive $`x`$ direction. Fig-8a shows the first and second moments of the probability with time in a thin slab of thickness $`2l^{}`$. The first and second moments initially increase as in an unbounded medium until the photon-front hits the boundary and dip before increasing again and saturating at an almost constant value. The dips occur because just after the ballistic and near ballistic components exit the slab, only the photons which are effectively moving in the opposite directions are left behind. In fact, the first moment is seen to become negative, implying that the net transport is in the backward direction for some time. The dip in the second moment implies that the photon cloud is effectively expanding at a slower rate. This would cause a lowered ’effective diffusion coefficient’ to be measured in a pulse transmission measurement. This reinforces the conclusions reached in Ref. based on Monte-carlo simulations and explains the experimental results of Ref. on a more rigorous footing. Fig.8b shows the survival probability for the case of a finite slab. This decays considerably faster than the in case of the semi-infinite slab, though at early times ($`tt^{}`$) the decay rates are comparable. The initial rates of decay are comparable because of the forward peaked nature of the probability distribution at early times, when the effect of the boundary at the back is hardly felt. This is to be compared with the mirror-image solution in the diffusion approximation, where equal weightage is given to both boundaries at all times.   

Finally we turn to the case of an amplifying stochastic medium. The effect of medium gain can be incorporated straight forwardly by noting that in our treatment, the time of exit from the slab directly translates into path-length traversed within the medium because speed is kept absolutely fixed. In the presence of amplification in the medium therefore, the net gain is directly proportional to the time. Thus the output flux at the boundary in a given direction is simply $`P(L,\theta ,t)\mathrm{cos}\theta \mathrm{exp}(\alpha t)`$, where $`\alpha `$ is the gain coefficient in the medium. It is thus simple to obtain a picture of amplified emission from such a medium. In Fig-9, we show the total light emitted by slab with boundaries at $`x=\pm 2l^{}`$ for several amplification factors. The photon is released from the origin in the positive $`x`$ \- direction. For large times, the output increases exponentially because of the presence of a exponential gain in the medium with no saturation. It is seen that the ballistic part is only slightly amplified while the output in the tail regions are increased considerably. To obtain a more realistic picture of lasing in random media  however, one would have to consider the lasing level population depletion and saturation effects.

## acknowledgement

SAR would like to thank Prof. Rajaram Nityananda for very helpful discussions.

## Appendix : Expression for the moments $`x^n`$

The $`n^{th}`$ order moment is given by

$$(xx_0)^n=c^n_0^t_0^t\mathrm{}_0^t𝑑t_n𝑑t_{n1}\mathrm{}𝑑t_1\mathrm{cos}\theta (t_1)\mathrm{cos}\theta (t_2)\mathrm{}\mathrm{cos}\theta (t_n)$$
(37)   
Writing $`\theta (t_i)`$ as $`\theta _i`$, the quantity within the angular brackets can be expressed as follows,

$`\mathrm{cos}\theta _1\mathrm{cos}\theta _2\mathrm{}\mathrm{cos}\theta _n`$ $`=`$ $`2^n(e^{i\theta _1}+e^{i\theta _1})(e^{i\theta _2}+e^{i\theta _2})\mathrm{}(e^{i\theta _n}+e^{i\theta _n})`$ (38)   
$`=`$ $`2^n{\displaystyle \underset{\stackrel{\sigma _1,\sigma _2\mathrm{}\sigma _n}{\sigma _i=\pm 1}}{}}\mathrm{exp}[i{\displaystyle \underset{j=1}{\overset{n}{}}}\sigma _j\theta _j]`$ (39)   
This can be expressed as a path integral using a gaussian distribution for $`f(t)`$.

$`\mathrm{exp}[i{\displaystyle \underset{j=1}{\overset{n}{}}}\sigma _j\theta _j]`$ $`=`$ $`{\displaystyle 𝒟[f(t)]\mathrm{exp}\left\{_0^t\left[\frac{f^2(t^{})}{2\mathrm{\Gamma }}+i\underset{j=1}{\overset{n}{}}\sigma _jf(t^{})\right]𝑑t^{}\right\}}`$ (40)   
$`=`$ $`exp\left\{{\displaystyle \frac{\mathrm{\Gamma }}{2}}{\displaystyle \underset{k=1}{\overset{n}{}}}\left({\displaystyle \underset{j=k}{\overset{n}{}}}\sigma _j\right)^2(t_kr_{k1})\right\}`$ (41)   
where $`t_0=0`$ and we assumed a time ordering of $`t_1<t_2<\mathrm{}<t_n`$. Thus,

$$(xx_0)^n=c^n(n!)2^n_0^t𝑑t_n_0^{t_n}𝑑t_{n1}\mathrm{}_0^{t_2}𝑑t_1\underset{\stackrel{\sigma _1,\sigma _2\mathrm{}\sigma _n}{\sigma _i=\pm 1}}{}\mathrm{exp}\left\{\frac{\mathrm{\Gamma }}{2}\underset{k=1}{\overset{n}{}}\left(\underset{j=k}{\overset{n}{}}\sigma _j\right)^2(t_kr_{k1})\right\}$$
(42)   
A similiar expression can be obtained for the $`(yy_0)^n`$ by noting that $`\mathrm{sin}\theta =\mathrm{cos}(\pi /2\theta )`$.   

Now we can obtain the joint probability distribution of $`x`$ and $`y`$ as

$$P(x,y,t;x_0,y_0,0)=\delta (xx_0c_0^t\mathrm{cos}\theta (t^{})𝑑t^{})\delta (yy_0c_0^t\mathrm{sin}\theta (t^{})𝑑t^{})$$
(43)   
Expressing the $`\delta `$-functions in terms of the Fourier transforms,

$$P(x,y,t;x_0,y_0,0)=\left(\frac{1}{2\pi }\right)^2_{\mathrm{}}^{\mathrm{}}𝑑\omega _x_{\mathrm{}}^{\mathrm{}}𝑑\omega _ye^{i[\omega _x(xx_0)+\omega _y(yy_0)]}\mathrm{exp}\left[ic_0^t[\omega _x\mathrm{cos}\theta (t^{})+\omega _y\mathrm{sin}\theta (t^{})]𝑑t^{}\right]$$
(44)   
This statistical average can be evaluated by a cumulant expansion  and since we have an expression for moments of all orders, we can in principle evaluate the cumulant expansion to any desired order.

FIGURE CAPTIONS   

Figure-1 : The marginal probability distributions $`P(x,t;x_0,0)`$ predicted by the approximate solution given by equation(28) at different times indicated in the figure. There is a clear cut-off at the light front and initially the probability accumulates at the light front (for $`t=t^{}`$).   

Figure-2 : The Probability distributions in the phase space of a particle in an infinite medium at different times obtained by numerically propagating equation(17). The particle is released at $`x=0`$ along the positive $`x`$ direction($`\theta =0`$) at $`t=0`$. The probability distribution is clearly forward peaked and and becomes almost flat along the $`\theta `$-axis only at times of about $`8t^{}`$.   

Figure-3 : The first and second moments of the displacement for the probability distribution of a particle in an infinite medium. The solid lines show the anlytical result of equation(20) and equation(22) while the symbol ( $``$ ) show the result obtained from the numerical solutions.   

Figure-4 : The marginal probability distribution $`P(x,t;x_0,\theta _0,0)=_\pi ^\pi P(x,\theta ,t;x_0\theta _0,0)𝑑\theta `$ at different times. The marginal probability distribution becomes almost a Gaussian at times of $`8t^{}`$.   

Figure-5 : The Probability distributions in the phase space of a particle in a semi-infinite medium at different times. The particle is released at $`x=0`$ along the positive $`x`$ direction($`\theta =0`$) at $`t=0`$. The absorbing boundary is located at $`4l^{}`$. The probability distribution is zero in the range $`\pi <\theta <\pi /2`$ and $`\pi /2<\theta <\pi `$ at the boundary implying there is no incoming flux into the medium.   

Figure-6 : The marginal probability distributions in a semi-infinite medium with an absorbing boundary at $`x=4l^{}`$. The plot on the right shows an expanded view of the distributions near the boundary. The solid straight lines are the linear extrapolations of the behaviour near the boundary. All of them are seen to cross the $`x`$-axis roughly at $`0.7l^{}`$ outside the boundary.   

Figure-7 : The surviving probability of the particle inside the semi-infinite medium for an absorbing boundary at $`4l^{}()`$ and $`2l^{}()`$. The persistence exponent $`\vartheta `$ is obtained from the long time behaviour of the survival probability. The lines ($`\mathrm{}`$) and ($``$) show the linear fits and give a persistence exponent of 0.4309 and 0.4364 respectively.   

Figure-8 : The first moments of the displacement of a particle in a finite slab of thickness $`2l^{}`$ (left plot). The right plot shows the survival probability in a semi-infinite medium and a finite slab. The distance between the point where the particle is released and the boundary is same in both case. ($`2l^{}`$).   

Figure-9 : Total light emitted (from both sides) by a disordered slab with amplification for different values of the gain coefficient $`A`$ in the medium. The output increases exponentially at long times.