**Journal of Modern Physics**

Vol.08 No.11(2017), Article ID:79809,41 pages

10.4236/jmp.2017.811108

Brownian Motion of Decaying Particles: Transition Probability, Computer Simulation, and First-Passage Times

M. P. Silverman^{ }

Department of Physics, Trinity College, Hartford, CT, USA

Copyright © 2017 by author and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

Received: September 18, 2017; Accepted: October 21, 2017; Published: October 24, 2017

ABSTRACT

Recent developments in the measurement of radioactive gases in passive diffusion motivate the analysis of Brownian motion of decaying particles, a subject that has received little previous attention. This paper reports the derivation and solution of equations comparable to the Fokker-Planck and Langevin equations for one-dimensional diffusion and decay of unstable particles. In marked contrast to the case of stable particles, the two equations are not equivalent, but provide different information regarding the same stochastic process. The differences arise because Brownian motion with particle decay is not a continuous process. The discontinuity is readily apparent in the computer-simulated trajectories of the Langevin equation that incorporate both a Wiener process for displacement fluctuations and a Bernoulli process for random decay. This paper also reports the derivation of the mean time of first passage of the decaying particle to absorbing boundaries. Here, too, particle decay can lead to an outcome markedly different from that for stable particles. In particular, the first-passage time of the decaying particle is always finite, whereas the time for a stable particle to reach a single absorbing boundary is theoretically infinite due to the heavy tail of the inverse Gaussian density. The methodology developed in this paper should prove useful in the investigation of radioactive gases, aerosols of radioactive atoms, dust particles to which adhere radioactive ions, as well as diffusing gases and liquids of unstable molecules.

**Keywords:**

Brownian Motion, Random Walk, Diffusion, Radioactivity, Transition Probability, Fokker-Planck Equation, Langevin, Equation, First Passage, Fick’s Law, Wiener Process

1. Introduction: Two Approaches to Brownian Motion

Brownian motion, one of the simplest examples of a random walk, is a nonequilibrium statistical process the mathematics of which serves to model a wide variety of stochastic processes throughout the physical and social sciences. From the earliest applications of Einstein to the random motion of small particles in a fluid [1] and of Bachelier to price fluctuations of the stock market [2] up to the present day, diffusive processes have been investigated intensively and reported in numerous monographs, textbooks, and articles of which the following provide representative examples of methods and applications [3] - [9] .

Despite their great diversity, nearly all such studies known to the author have in common the feature that the diffusing particles maintain their identity throughout the stochastic process. There are a few exceptions, such neutron diffusion with beta decay in reactor materials [10] [11] or the chemical transformation of diffusing reactants leading to pattern formation [12] [13] [14] , but the content and objectives of these investigations are very different from those of this paper, which was motivated by recent experimental studies of the diffusion of radioactive gases [15] . In contrast to previous treatments of diffusion, the subject matter of this paper concerns the transition probabilities, statistical moments, and computer simulation of the Brownian motion of particles that can decay randomly. To underscore the distinction further, note that the decay terms that may appear in the equations of motion of well-studied stochastic processes such as the Ornstein-Uhlenbeck process [16] arise from the frictional interaction of a stable particle with the environment and involve only the process variables, like velocity, of the diffusing particle. In the Brownian motion treated in this paper, there is no frictional interaction with the environment, and it is the diffusing particles themselves that decay.

Under the assumption, usually justified by the Fermi Golden Rule of time- dependent perturbation theory in quantum mechanics [17] , that the probability of particle decay within a short time interval $\Delta t$ takes the form

$p\left(\Delta t\right)=\lambda \Delta t\ll 1$ , (1)

where $\lambda $ is the intrinsic decay rate, it is readily deducible [18] that the probability density of survival to time t satisfies the differential equation

$\text{d}p\left(t\right)/\text{d}t=-\lambda p\left(t\right)$ (2)

with normalized solution

$p\left(t\right)=\lambda {\text{e}}^{-\lambda t}$ . (3)

From Equation (3) it follows that the mean lifetime $\tau $ is the reciprocal of $\lambda $

$\tau ={\displaystyle {\int}_{0}^{\infty}tp\left(t\right)\text{d}t}={\lambda}^{-1}$ , (4)

and that the half-life, i.e. duration within which the survival probability is 50%, is

${\tau}_{1/2}={\lambda}^{-1}\mathrm{ln}2=\tau \mathrm{ln}2$ . (5)

The premise that the probability (1), with $\lambda $ independent of environmental interactions, accurately characterizes the transmutation of radioactive nuclei has been challenged on various grounds during the past 20 years. However, in each case the theoretical consequences of relation (1) were validated experimentally and shown to be in accord with currently known laws of physics [19] [20] .

Because unstable particles are removed from the system at random instances during the time they are moving randomly within a specified sample space, the familiar standard equations (e.g. Langevin and Fokker-Planck) for Brownian motion must be re-derived on the basis of a conservation law that takes particle loss into account. This more general relation, presented in Section 2, leads to significant differences in the analytical structure of the stochastic equations and statistical moments compared with corresponding expressions derived for the random walk of a stable particle.

In general, there are two mathematically different approaches to analyzing Brownian motion. One method, to be referred to as the Langevin approach, is to examine the process variables directly. For a particle undergoing a one-dimensional random walk, the process variable of interest in this paper is the displacement $X\left(t\right)$ as a function of time t. Conventional symbolic notation, which is used in this paper unless stated otherwise, employs an upper-case letter (e.g. X) to designate a random variable and a corresponding lower-case letter (e.g. x) to represent a realization or measurement of the variable (referred to as a variate in statistical terminology). The Langevin equation of motion for the time-evolution of the coordinate variable usually takes the form of an Ito stochastic differential equation [21]

$\text{d}X\left(t\right)\equiv X\left(t+\text{d}t\right)-X\left(t\right)=A\left(X,t\right)\text{d}t+\sqrt{B\left(X,t\right)}\text{\hspace{0.05em}}\text{d}W\left(t\right)$ (6)

in which $A\left(X,t\right)$ is the drift function, $B\left(X,t\right)$ is the diffusivity, and $\text{d}W\left(t\right)=\sqrt{\text{d}t}N\left(0,1\right)$ is the differential Wiener process in which $N\left(0,1\right)$ symbolizes a Gaussian distribution of mean $\mu =0$ and variance ${\sigma}^{2}=1$ . The drift $A\left(X,t\right)$ is a deterministic function, whereas the Wiener term is the source of fluctuations.

It is clear from Equation (6) that the analytical solution, if one can be found, is not a deterministic function of a dynamical variable, but a probability distribution. A distribution is said to be stable if two independent random variables characterized by this distribution sum to form a random variable governed by a distribution of the same kind [22] . Gaussian and Lorentzian processes are stable, but most stochastic processes are not. For example, in regard to a Gaussian distribution whose probability density function (PDF) takes the general form

${p}_{X}\left(x|\mu ,\sigma \right)=\frac{1}{\sqrt{2\text{\pi}{\sigma}^{2}}}\mathrm{exp}\left(-{\left(x-\mu \right)}^{2}/2{\sigma}^{2}\right)$ , (7)

one can express the corresponding random variable as

$X=N\left(\mu ,{\sigma}^{2}\right)=\mu +\sigma N\left(0,1\right)$ (8)

and the sum of two independent Gaussian random variables as

$\begin{array}{c}{X}_{1}+{X}_{2}={N}_{1}\left({\mu}_{1},{\sigma}_{1}^{2}\right)+{N}_{2}\left({\mu}_{2},{\sigma}_{2}^{2}\right)=N\left({\mu}_{1}+{\mu}_{2},{\sigma}_{1}^{2}+{\sigma}_{2}^{2}\right)\\ ={\mu}_{1}+{\mu}_{2}+\sqrt{{\sigma}_{1}^{2}+{\sigma}_{2}^{2}}N\left(0,1\right).\end{array}$ (9)

It follows from the stability relation (9) that the solution to (6) is itself a Gaussian distribution.

The second method, to be referred to as the Fokker-Planck approach, is to solve for the transition probability function (TPF) $p\left(x,t|{x}_{0},{t}_{0}\right)$ of finding the particle at location x at time t, given that the particle was initially at ${x}_{0}$ at time ${t}_{0}<t$ . The partial differential equation corresponding to the stochastic differential Equation (6) is the forward Fokker-Planck equation [Ref [3] , pp. 117-122]

$\frac{\partial p\left(x,t|{x}_{0},{t}_{0}\right)}{\partial t}=-\frac{\partial \left(A\left(x,t\right)p\left(x,t|{x}_{0},{t}_{0}\right)\right)}{\partial x}+\frac{1}{2}\frac{{\partial}^{2}\left(B\left(x,t\right)p\left(x,t|{x}_{0},{t}_{0}\right)\right)}{\partial {x}^{2}}$ . (10)

For non-decaying particles, Equations ((6) and (10)) have the same information content, but provide different perspectives on the same stochastic process. For example, the stochastic differential Equation (6) is particularly suitable for numerical solution by an iterative up-dating algorithm, thereby permitting computer simulation of particle displacement and visualization of the process variables in real time. On the other hand, the forward Fokker-Planck Equation (10) yields a probability density function from which all statistical moments can be calculated. The solution $p\left(x,t|{x}_{0},{t}_{0}\right)$ also solves the backward Fokker- Planck equation [Ref [8] , pp. 168-171]

$\frac{\partial p\left(x,t|{x}_{0},{t}_{0}\right)}{\partial {t}_{0}}=-A\left({x}_{0},{t}_{0}\right)\frac{\partial p\left(x,t|{x}_{0},{t}_{0}\right)}{\partial {x}_{0}}-\frac{1}{2}B\left({x}_{0},{t}_{0}\right)\frac{{\partial}^{2}p\left(x,t|{x}_{0},{t}_{0}\right)}{\partial {x}_{0}^{2}}$ , (11)

in which derivatives are taken with respect to the initial coordinates. Equation (11) provides the same information as the forward Fokker-Planck Equation (10), but is of particular utility in the analysis of Brownian motion within a region confined by boundaries. The structure of Equation (11) facilitates treatment of the problem of first-passage times to be analyzed in Section 3 by alternative, simpler methods.

Generally speaking, one solves a Brownian motion problem to answer two kinds of questions: 1) How far has a particle randomly walked in a given time? and 2) In how much time will a particle randomly walk to a given location? It is to be understood, of course, that answers to questions like these are probabilistic, not deterministic, statements. Questions of the first kind are perhaps more familiar, but there are numerous circumstances that call for answers to questions of the second kind. These are usually referred to as first-passage processes [23] , and in statistical terminology the time to reach a designated location has been called first-passage time, first hitting time, first return time, exit time, and possibly other names depending on the exact circumstances.

This paper is concerned primarily with one-dimensional Brownian motion of a particle of finite statistical lifetime in free (unconstrained) space and in a space with absorbing boundaries. Derivation and analysis of the equations corresponding to the Fokker-Planck equation and Langevin equation will show that these equations do not lead to equivalent descriptions of the Brownian motion of a decaying particle, in marked contrast to the case of stable particles. This inequivalence is traceable to the fact that particle decay is a discontinuous process. The Fokker-Planck equation, which yields a transition probability density, remains a continuous differential equation, but the Langevin equation, which describes a process variable, and not a probability density, must directly manifest the discontinuous nature of decay.

The distinction between the Brownian motion of stable and unstable particles becomes very apparent in the problem of first-passage time to absorbing boundaries. Upon reaching an absorbing boundary for the first (and therefore only) time, the particle is removed from the system and the process is terminated. Likewise, upon decay at some unpredictable moment, the particle is also removed and the process is terminated. Thus, the instability of the particle introduces into a first-passage problem an additional time element that can radically modify expectation values. From the perspective of physics, this modification extends the applicability of first-passage time theory to a broader class of physical systems. As a physical model, the solution to a first-passage time problem with particle decay has been applied to the diffusion of the radioactive gas radon-222 in the atmosphere, and should likewise prove useful in the study of diffusion of other radioactive gases, radioactive ions that form as daughter products in radioactive decay, as well as unstable molecules that change their identity by chemical transformation.

As a purely stochastic problem, the analysis undertaken in this paper

・ derives and provides an exact solution to the Fokker-Planck and Langevin equations of Brownian motion of an unstable particle,

・ extends to unstable particles the two principal methods of calculating first-passage times,

・ demonstrates how to simulate by computer the Brownian motion of an unstable particle, and

・ clarifies a number of confusing issues that arise in the case of unstable particles (but not stable particles) regarding Fokker-Planck and Langevin equations, expectation values, probability density functions, and transition probability functions.

This paper is organized as follows. Section 2 is concerned with spatial aspects of a decaying particle undergoing a one-dimensional random walk. Derivation and solution of the Fokker-Planck equation to obtain the transition probability density and associated statistical moments are given in Section 2.1. In Section 2.2 the mean-square displacement of the decaying particle is reconsidered from the perspective of a random walk on a discrete lattice and shown to coincide in the appropriate limit with the result obtained in Section 2.1. The derivation, numerical solution, and computer simulation of the Langevin equation to obtain Brownian motion trajectories of the decaying particle are given in Section 2.3. In Section 2.4 the Langevin equation is solved analytically to obtain the distribution function of the particle displacement. Section 3 is concerned with temporal aspects of a decaying particle undergoing a one-dimensional random walk. The mean first-passage time to absorbing boundaries is solved in Section 3.1 by the method of image functions and in Section 3.2 by solution of a screened Poisson equation. Section 3.3 illustrates the critical role of particle decay in leading to first-passage time results that differ markedly from corresponding results for a stable particle. Section 4 examines the validity of the stochastic model, based on a Wiener process or Fick’s law, to account for fluctuations in the spatial displacement of a decaying particle. Finally, conclusions drawn from these analyses are summarized in Section 5.

2. Diffusion with Decay

2.1 Fokker-Planck Equation and Transition Probability

Consider a quantity $Q\left(t\right)={\displaystyle \underset{V}{\iiint}n\left(x,t\right)\text{d}V}$ of particles with decay constant $\lambda $ and number density $n\left(x,t\right)$ within a volume V bound by surface S. Since loss of Q can occur either by intrinsic decay at the rate $-\lambda Q$ or by diffusion of a current density $j\left(x,t\right)$ outward across surface S, macroscopic mass balance requires that

$\frac{\partial}{\partial t}{\displaystyle \underset{V}{\iiint}n\text{d}V}=-{\displaystyle \underset{S}{\u222f}j\cdot \text{d}S}-\lambda {\displaystyle \underset{V}{\iiint}n\text{d}V}$ . (12)

Application of the divergence theorem then leads from the integral relation (12) to the differential equation

$\frac{\partial n}{\partial t}=-\nabla \cdot j-\lambda n$ (13)

for the conservation law of a disintegrating quantity.

Upon dividing Equation (13) by the initial number of particles ${N}_{0}\gg 1$ and relating the current density to the gradient of the particle density by Fick’s law

$j=-D\nabla n$ , (14)

with diffusion coefficient D (here taken to be a constant), one obtains an equation of the form

$\frac{\partial w}{\partial t}=D{\nabla}^{2}w-\lambda w$ (15)

in which $w\left(x,t\right)=n\left(x,t\right)/{N}_{0}$ is interpretable as the probability density for diffusion of a single decaying particle. However, it is demonstrable that the solution $w\left(x,t\right)$ under the special initial condition at ${t}_{0}$

$w\left(x,{t}_{0}\right)=\delta \left(x-{x}_{0}\right)$ (16)

is identical to the transition probability density $p\left(x,t|{x}_{0},{t}_{0}\right)$ , which is the conditional probability for an unstable particle to have reached position $x$ at time t given that it was at ${x}_{0}$ at ${t}_{0}<t$ . The equivalence is readily established by substitution of condition (16) into the Chapman-Kolmogorov equation [24]

$\begin{array}{c}w\left(x,t\right)={\displaystyle {\int}_{-\infty}^{\infty}p\left(x,t|{x}_{1},{t}_{0}\right)w\left({x}_{1},{t}_{0}\right)\text{d}{x}_{1}}\\ ={\displaystyle {\int}_{-\infty}^{\infty}p\left(x,t|{x}_{1},{t}_{0}\right)\delta \left({x}_{1}-{x}_{0}\right)\text{d}{x}_{1}}=p\left(x,t|{x}_{0},{t}_{0}\right).\end{array}$ (17)

In the case of one-dimensional Brownian motion treated in this paper, Equation (15) then takes the form

$\frac{\partial p\left(x,t|{x}_{0},{t}_{0}\right)}{\partial t}=D\frac{{\partial}^{2}p\left(x,t|{x}_{0},{t}_{0}\right)}{{\partial}^{2}x}-\lambda p\left(x,t|{x}_{0},{t}_{0}\right)$ , (18)

which, together with the initial delta-function condition (16)

$p\left(x,{t}_{0}|{x}_{0},{t}_{0}\right)=\delta \left(x-{x}_{0}\right)$ (19)

comprises the equation for the Green’s function of a one-dimensional random walk with decay [25] .

The solution to Equations ((18) and (19)) can be obtained in several ways. One method is to take the Fourier transform with respect to the spatial coordinate, which converts the partial differential Equation (18) of two variables (space, time) into an ordinary differential equation in one variable (time). The simplest method, however, which has been used to solve the Schroedinger equation for transitions between excited atomic states [26] and the rate equations for a sequence of nuclear transformations [27] , is to eliminate the decay term in Equation (18) by the substitution $p\left(x,t|{x}_{0},{t}_{0}\right)={p}_{0}\left(x,t|{x}_{0},{t}_{0}\right)\mathrm{exp}\left(-\lambda \left(t-{t}_{0}\right)\right)$ , thereby transforming Equation (18) into the equation for the Green’s function ${p}_{0}\left(x,t|{x}_{0},{t}_{0}\right)$ of a non-decaying diffusing particle, the solution of which is known [Ref [6] , p. 12]. Either way, one obtains the transition probability function (TPF)

$p\left(x,t|{x}_{0},{t}_{0}\right)=\frac{1}{\sqrt{4\text{\pi}D\left(t-{t}_{0}\right)}}\mathrm{exp}\left(-{\left(x-{x}_{0}\right)}^{2}/4D\left(t-{t}_{0}\right)\right){\text{e}}^{-\lambda \left(t-{t}_{0}\right)}$ . (20)

From the form of relation (20), one sees that the transition probability is a function of the time interval $t-{t}_{0}$ , and not the separate time coordinates $\left(t,{t}_{0}\right)$ . From this point on, it will be assumed that ${t}_{0}=0$ and t will represent a time interval.

Although the functions $w\left(x,t\right)$ and $p\left(x,t|{x}_{0},0\right)$ are mathematically identical in form, they serve different purposes and are used to calculate different quantities. For example, the probability density function (PDF) is employed to calculate the statistical moments $\left\{{m}_{k}\right\}$ of a distribution, where the ${k}^{\text{th}}$ moment (expectation value) with respect to the origin is defined by

${m}_{k}={\displaystyle \int {x}^{k}w\left(x,t\right)\text{d}x}/{\displaystyle \int w\left(x,t\right)\text{d}x}$ , (21)

and the range of integration covers the defined sample space. If the PDF is normalized, then the denominator in relation (21) is unity. However, the normalization of PDF (20) is not unity

${\int}_{-\infty}^{\infty}w\left(x,t\right)\text{d}x}={\displaystyle {\int}_{-\infty}^{\infty}p\left(x,t|{x}_{0}\right)\text{d}x}={\text{e}}^{-\lambda t$ , (22)

but yields, instead, the survival probability to time t of the unstable particle. In contrast to (21), the TPF is employed to calculate transition probabilities and statistical moments of a different nature. For example, the expectation value of the ${k}^{\text{th}}$ power of displacement of a particle in time interval t starting from position ${x}_{0}$ is

$\langle {x}^{k}\rangle ={\displaystyle \int {x}^{k}p\left(x,{t}_{0}|{x}_{0},0\right)\text{d}x}={m}_{k}{\text{e}}^{-\lambda t}$ . (23)

Because the TPF is a conditional probability, the expectation value (23) does not include a normalizing denominator as in Equation (21). For a stable particle, relations (21) and (23) yield the same result, but this is not the case for a decaying particle.

The distinction between PDF and TPF leads to different results for the mean- square displacement. Employing solution (20) as the Gaussian PDF $w\left(x,t\right)$ in the expectation value (21), one obtains the mean location ${m}_{1}={x}_{0}$ and variance about the mean

${m}_{2}-{m}_{1}^{2}=2Dt$ (24)

of a non-decaying diffusing particle. However, calculation of expectation values (23) with the TPF $p\left(x,t|{x}_{0},0\right)$ leads to the initial location $\langle x\rangle ={x}_{0}$ of the diffusing particle and its mean-square distance from the point of origin

${\sigma}_{x}^{2}\left(t\right)=2Dt\text{\hspace{0.05em}}{\text{e}}^{-\lambda t}$ . (25)

Thus, the same mathematical function (20) generates expectation values that differ in interpretation and mathematical form depending on what is sought by the analyst. In the context of understanding how decay affects the probability of displacement of a single particle in continuous Brownian motion, relation (25) is the relevant quantity, as shown in the following section.

2.2. Mean-Square Displacement

Consider a one-dimensional Gaussian random walk on a lattice with time step $\delta t$ and displacement $\delta {X}_{j}$ during the jth time step given by

$\delta {X}_{j}=\delta X{N}_{j}\left(0,1\right)$ (26)

where the lattice spacing $\delta X$ sets the scale of spatial displacement. Each displacement is taken to be an independent Gaussian random variable. However, a displacement can be made only if the particle has survived during that time step. From Equation (3) it is seen that the probability of survival during a time step $\delta t$ is ${p}_{s}=1-\lambda \delta t$ , and the probability that the process ends at a particular time step is $\lambda \delta t$ . Thus, although the distance traveled in the jth time step is determined by a unit normal distribution ${N}_{j}\left(0,1\right)$ , the probability that the displacement is made at all is determined by a Bernoulli random variable

${\epsilon}_{j}={B}_{j}\left(1,{p}_{s}\right)=\{\begin{array}{l}1\text{withprobability}{p}_{s}=1-\lambda \delta t\\ 0\text{withprobability}1-{p}_{s}=\lambda \delta t\end{array}$ . (27)

A Bernoulli distribution is the special case $n=1$ of the binomial distribution $B\left(n,p\right)$ of n trials with probability of success $p$ ; the corresponding discrete probability function is

${p}_{\text{B}}\left(k|n,p\right)=\left(\begin{array}{l}n\\ k\end{array}\right){p}^{k}{\left(1-p\right)}^{n-k}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(n\ge k\ge 0\right)$ . (28)

The subscript j in relations (26) and (27) denotes that each random sample, whether Gaussian or Bernoulli, is independent of all the others.

The displacement after n time steps is given by the random variable ${X}_{n}$

${X}_{n}={\displaystyle \underset{j=1}{\overset{n}{\sum}}\delta {X}_{j}}=\delta X{\displaystyle \underset{j=1}{\overset{n}{\sum}}{N}_{j}\left(0,1\right)}=\delta X\text{\hspace{0.05em}}N\left(0,n\right)$ . (29)

The expectation value of the mean square displacement at time $t=n\delta t$ is therefore

$\langle {X}_{n}^{2}\rangle ={\left(\delta X\right)}^{2}\langle {\left(N\left(0,n\right)\right)}^{2}\rangle {p}_{s}^{n}=n{\left(\delta X\right)}^{2}{\left(1-\lambda \delta t\right)}^{n}$ , (30)

which can be re-expressed in the form

$\langle {X}_{n}^{2}\rangle =t\left(\frac{{\left(\delta X\right)}^{2}}{\delta t}\right){\left(1-\lambda \frac{t}{n}\right)}^{n}$ . (31)

Upon defining the diffusion constant D in the standard way

$2D={\left(\delta X\right)}^{2}/\delta t$ (32)

and taking the limit $\delta t\to 0$ , $\delta X\to 0$ , with requirement that D remain finite, Equation (31) becomes

$\langle {X}_{t}^{2}\rangle =2D\text{\hspace{0.05em}}t{\text{e}}^{-\lambda t}$ (33)

in accord with the expectation value (25) obtained directly from the TPF (20).

As a point of clarification, the reason for the factor of 2 in the conventional definition (32) of the diffusion coefficient is that the quantity $2D$ corresponds to the fluctuation function $B\left(x,t\right)$ in the Wiener term of the forward Fokker- Planck Equation (10). Under the conditions that $B\left(x,t\right)=2D$ is constant and there is no drift, $A\left(x,t\right)=0$ , Equation (10) then reduces to the standard one- dimensional diffusion equation

$\frac{\partial p\left(x,t|{x}_{0},{t}_{0}\right)}{\partial t}=D\frac{{\partial}^{2}p\left(x,t|{x}_{0},{t}_{0}\right)}{\partial {x}^{2}}$ (34)

in which D is the physically measurable diffusion constant first introduced by Einstein in his theory of Brownian motion.

2.3. Langevin Equation: Update Algorithm and Computer Simulation

The Langevin Equation (6) for one-dimensional Brownian motion of a non-decaying particle in the absence of drift is expressible in a form

$x\left(t+\text{d}t\right)=x\left(t\right)+\sqrt{2D\text{d}t}{n}_{t}$ (35)

that facilitates numerical solution by an update algorithm. The lower-case letters x in Equation (35) signify numerical realizations of the random variable X in Equation (6); dt is the numerical time step; and ${n}_{t}$ is a random sample from the unit normal distribution ${N}_{t}^{t+\text{d}t}\left(0,1\right)$ , where the subscript t and superscript $t+\text{d}t$ explicitly denote the temporal range with respect to which ${n}_{t}$ is associated. Thus, two samples ${n}_{{t}_{1}}$ and ${n}_{{t}_{2}}$ , corresponding to distributions ${N}_{{t}_{1}}^{{t}_{1}+\text{d}t}\left(0,1\right)$ and ${N}_{{t}_{2}}^{{t}_{2}+\text{d}t}\left(0,1\right)$ , are independent for $\left|{t}_{2}-{t}_{1}\right|>\text{d}t$ . Given an initial value $x\left(0\right)$ , a sequence of points $\left\{x\left(0\right),x\left(\text{d}t\right),x\left(2\text{d}t\right),\cdots ,x\left(n\text{d}t\right)\right\}$ is generated by iterative use of Equation (35) up to time $t=n\text{d}t$ . (Note: The symbol n without subscript is the number of time steps, not a sample from a normal distribution.) The sequence of points is an approximation to the true Brownian motion in the limit $\text{d}t\to 0$ . In that theoretical limit, the trajectory of Brownian motion is a curve that is everywhere continuous, but nowhere differentiable; in other words, a particle trajectory for which the particle velocity is undefined.

Although Equation (35) is simple enough to be solved analytically, a numerical solution provides a graphical visualization of Brownian motion paths. Moreover, starting from the same initial condition and generating numerous Brownian trajectories for a specified time interval t provides a Gibbs ensemble [28] [29] , from which ensemble and time averages of moments, correlation functions, and other statistical quantities can be obtained by Monte Carlo methods [30] [31] . Such methods are particularly useful when the Langevin or Fokker-Planck equations cannot be solved analytically.

The question addressed in this section is this: What is the Langevin equation for Brownian motion of a decaying particle? Since the Langevin and Fokker- Planck equations of a stable particle ordinarily provide equivalent information, one can in principle start with either one and obtain the functions $A\left(x,t\right)$ and $B\left(x,t\right)$ needed for the other. Thus, for a stable particle the Langevin Equation (35) leads to the corresponding Fokker-Planck Equation (34), and vice-versa. It is to be stressed, however, that only a continuous Markov process is completely and equivalently defined by either the Langevin equation or Fokker-Planck equation [Ref [8] , p. 158]. In the case of a decaying particle, the Brownian motion is not a continuous process because it is abruptly terminated by disintegration of the particle. Equation (18)―although referred to in this paper as a Fokker- Planck equation in deference to conventional usage―is actually not in the form of a Fokker-Planck equation. The term $\lambda p\left(x,t|{x}_{0},{t}_{0}\right)$ is not a drift term since it does not involve the spatial derivative of $p\left(x,t|{x}_{0},{t}_{0}\right)$ . Here, then, is another important difference in the Brownian motion of a particle that decays randomly compared to one that is stable.

Theoretically, it is possible to transform Equation (18) into a Fokker-Planck equation. One merely needs to find a drift function $A\left(x,t\right)$ that satisfies the differential equation

$\lambda \frac{\partial}{\partial x}\left(A\left(x,t\right)p\left(x,t|{x}_{0},{t}_{0}\right)\right)=\lambda p\left(x,t|{x}_{0},{t}_{0}\right)$ . (36)

Integration of Equation (36) yields

$A\left(x,t\right)=\sqrt{\text{\pi}Dt}\text{erf}\left(\frac{x}{\sqrt{4Dt}}\right)\mathrm{exp}\left(\frac{{x}^{2}}{4Dt}\right)$ (37)

with boundary and initial conditions

$A\left(0,t\right)=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}A\left(x,0\right)\sim \delta \left(x\right)$ . (38)

The error function is defined by the integral

$\text{erf}\left(x\right)\equiv \frac{2}{\sqrt{\text{\pi}}}{\displaystyle {\int}_{0}^{x}{\text{e}}^{-{u}^{2}}\text{d}u}=-\text{erf}\left(-x\right)$ . (39)

However, the preceding solution (37)―or, indeed, any transformation that generates a Fokker-Planck equation from Equation (18) by finding a drift term―seriously misrepresents the physics of the problem. This is a stochastic process, as emphasized in the preceding sections, in which the physical particle (or the probability of particle existence), and not a process variable, is decaying. There is neither drift nor friction in this process.

The key point to recognize in constructing an appropriate stochastic differential equation―which, for the sake of conventional nomenclature, is referred to in this paper as a Langevin equation, even though rigorously it is not―is that the unstable particle, as long as it exists, undergoes Brownian motion as described by the standard diffusion Equation (35). The Brownian motion does not continue indefinitely, however, but is interrupted randomly by particle decay. The first instance of decay terminates the process; there is no further diffusion of that particle. (One could then introduce another particle and follow its Brownian motion if it is desired to acquire an ensemble of trajectories.) The stochastic process, therefore, entails sequences of two paired independent distributions: (a) the unit normal distribution ${N}_{t}^{t+\text{d}t}\left(0,1\right)$ , which determines the direction and extent of displacement during the time period $\left[t,t+\text{d}t\right]$ , and (b) the Bernoulli distribution ${B}_{t}^{t+\text{d}t}\left(1,{p}_{s}\right)$ , which determines whether or not the particle survives the interval dt with a survival probability ${p}_{s}$ given by Equation (27).

The appropriate Langevin equation would then take the update form

$x\left(t+\text{d}t\right)=\left[x\left(t\right)+\sqrt{2D\text{d}t}{n}_{t}\right]{\epsilon}_{t}$ (40)

where the Bernoulli variate ${\epsilon}_{t}$ is defined in relation (27). Note that an occurrence of ${\epsilon}_{t}=0$ terminates the entire process by setting $x\left(t+\text{d}t\right)$ to zero. Thus, just before the moment of decay, the particle will have reached location $x\left(t\right)$ and no further. Although the Langevin Equation (40) and the TPF Equation (18) both describe Brownian motion of a decaying particle, the descriptions are not equivalent. The TPF is a continuous probability function for all displacements x of a given particle; there is no built-in mechanism to terminate the process at decay. The Langevin Equation (40) has such a mechanism.

To get a sense of the progression of the process, the first few iterations of (40) for a particle starting at the origin ${x}_{0}=0$ are shown explicitly below:

$\begin{array}{l}x\left(\text{d}t\right)=\sqrt{2D\text{d}t}{n}_{1}{\epsilon}_{1}\\ x\left(2\text{d}t\right)=\sqrt{2D\text{d}t}\left({n}_{1}{\epsilon}_{1}{\epsilon}_{2}+{n}_{2}{\epsilon}_{2}\right)\\ x\left(3\text{d}t\right)=\sqrt{2D\text{d}t}\left({n}_{1}{\epsilon}_{1}{\epsilon}_{2}{\epsilon}_{3}+{n}_{2}{\epsilon}_{2}{\epsilon}_{3}+{n}_{3}{\epsilon}_{3}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\vdots \\ x\left(n\text{d}t\right)=\sqrt{2D\text{d}t}\left({n}_{1}{\displaystyle \underset{j=1}{\overset{n}{\prod}}{\epsilon}_{j}}+{n}_{2}{\displaystyle \underset{j=2}{\overset{n}{\prod}}{\epsilon}_{j}}+{n}_{3}{\displaystyle \underset{j=3}{\overset{n}{\prod}}{\epsilon}_{j}}+\cdots +{n}_{n}{\epsilon}_{n}\right).\end{array}$ (41)

The mean displacement $\langle x\left(n\text{d}t\right)\rangle =0$ follows immediately from the properties of the unit Gaussian $\langle {n}_{t}\rangle =\langle {N}_{t}^{t+\text{d}t}\left(0,1\right)\rangle =0$ . The mean-square displacement is less obvious and leads to two results depending on whether one retains the structure of discrete time steps or takes the limit for continuous displacement in time.

Consider first the continuous-time limit of the mean-square displacement, where one sets $\text{d}t=t/n$ and eventually takes the limit $n\to \infty $ and $\text{d}t\to 0$ such that t is a fixed quantity. Equation (41) then yields

$\begin{array}{l}\langle {x}^{2}\left(n\text{d}t\right)\rangle \\ =\frac{2Dt}{n}\left(\langle {\left({n}_{1}\right)}^{2}\rangle {\displaystyle \prod _{j=1}^{n}\langle {\left({\epsilon}_{j}\right)}^{2}\rangle}+\langle {\left({n}_{2}\right)}^{2}\rangle {\displaystyle \prod _{j=2}^{n}\langle {\left({\epsilon}_{j}\right)}^{2}\rangle}+\cdots +\langle {\left({n}_{n}\right)}^{2}\rangle \langle {\left({\epsilon}_{n}\right)}^{2}\rangle \right)\end{array}$ (42)

where the Bernoulli variates ${\epsilon}_{t}=1$ at each time step (otherwise there would not be n steps). There are no cross-terms in the expectation values because all the Gaussian and Bernoulli samples are realizations of independent uncorrelated random variables. Upon insertion in Equation (42) of the expectation values

$\begin{array}{l}\langle {\left({n}_{t}\right)}^{2}\rangle =1\\ \langle {\left({\epsilon}_{t}\right)}^{2}\rangle ={p}_{s}=\left(1-\lambda \frac{t}{n}\right),\end{array}$ (43)

summing the terms in Equation (42), and taking the forementioned limit, one arrives at

$\langle {x}^{2}\left(t\right)\rangle =\frac{2Dt}{n}{{\displaystyle \underset{j=1}{\overset{n}{\sum}}\left(1-\frac{\lambda t}{n}\right)}}^{j}\underset{n\to \infty}{\to}\frac{2D}{\lambda}\left(1-{\text{e}}^{-\lambda t}\right)$ . (44)

It should not be surprising that expression (44) differs from the mean-square displacement (25) derived from the TPF, because Equation (18) provides different information than the stochastic differential Equation (40). As shown in Section 2.2, the expectation value (25) is equivalent to the mean-square displacement of a particle prior to decay multiplied by the survival probability to time t. Statistically, it may be thought of as a compound stochastic process $N\left(0,n\right)\times $ $B\left(n,{p}_{s}\right)$ . By comparison, one can think of the outcome (44) as resulting from a sum of n compound stochastic processes of the form $N\left(0,1\right)\times B\left(1,{p}_{s}\right)$ .

It is of interest to examine the two limiting cases of Equation (44):

$\langle {x}^{2}\left(t\right)\rangle \to \{\begin{array}{l}2Dt\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\left(\lambda t\ll 1\right)\\ 2D/\lambda \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(\lambda t\gg 1\right)\end{array}$ . (45)

The first limit in (45) shows that the mean-square displacement of a particle whose statistical lifetime is much longer than the diffusion time is the same as for Brownian motion of a stable particle. The second limit shows that the root- mean-square distance $\sqrt{\langle {x}^{2}\left(t\right)\rangle}$ reached by a particle very likely to decay during the diffusion time is equivalent, within a factor $\sqrt{2}$ , to the characteristic diffusion length [15]

$\zeta =\sqrt{D/\lambda}$ (46)

that occurs in the solution of the diffusion equation for radioactive gases [32] . Recall, however, that $2D$ , rather than $D$ , is the mathematical diffusivity equivalent to the function $B\left(x,t\right)$ appearing in the Fokker-Planck Equation (10). In this paper, therefore, the mathematical diffusion length will be defined as ${\zeta}_{\text{m}}\equiv \sqrt{2D/\lambda}$ .

Consider next the mean-square displacement as obtained from Equation (44) with n discrete time steps of duration $\delta t$

$\langle {x}^{2}\left(n\delta t\right)\rangle =2D\delta t{\displaystyle \underset{j=1}{\overset{n}{\sum}}{p}_{s}^{j}}=\frac{2D}{\lambda}\left(1-\lambda \delta t\right)\left[1-{\left(1-\lambda \delta t\right)}^{n}\right]$ , (47)

where use was made of the summation formula for a geometric series

$\underset{j=1}{\overset{n}{\sum}}{p}^{j}}={\displaystyle \underset{j=0}{\overset{n}{\sum}}{p}^{j}}-1=\frac{1-{p}^{n+1}}{1-p}-1$ . (48)

For $\lambda \delta t\ll 1$ , in accordance with the assumption underlying Equation (1), Equation (47) can be reduced to

$\langle {x}^{2}\left(n\delta t\right)\rangle \approx 2Dn\delta t=2Dt$ (49)

by application of the approximation ${\left(1-\lambda \delta t\right)}^{n}\approx 1-n\lambda \delta t$ and neglect of the term $\lambda \delta t\ll 1$ where it occurs alone, i.e. not multiplied by n.

The two methods of taking limits led to two different outcomes, (44) vs. (49), because each method held a different quantity constant. In the approach leading to (44), the total diffusion time t was fixed, and the number of time steps n was taken to an infinite limit as time interval $\delta t$ approached zero. In other words, n was merely an intermediary variable; finite values of n did not define the duration t of the process. However, in the approach leading to (49), the number of time steps n is the quantity of interest, and the duration $t=n\delta t$ is determined by n.

Because the decay of the diffusing particle can occur randomly at any time step, the number n in Equation (49) is itself a random quantity unknown at the outset of a Brownian random walk. One can therefore regard n as a realization of a discrete random variable N (not to be confounded with the symbol $N\left(0,1\right)$ for a unit normal distribution) that is subject to a geometric distribution with probability function

${p}_{N}\left(n|{p}_{s}\right)={p}_{s}^{n}{q}_{s}={p}_{s}^{n}\left(1-{p}_{s}\right)$ (50)

in which ${p}_{s}=1-\lambda \delta t$ is the probability of success (or survival) at each time step and ${q}_{s}=\lambda \delta t$ is the probability of failure (or decay). Thus, Equation (50) expresses the probability of a process with n successes in sequence followed by a single failure that terminates the process. The distribution is normalized

$\underset{n=0}{\overset{\infty}{\sum}}{p}_{N}\left(n|{p}_{s}\right)}=1$ (51)

with mean $\stackrel{\xaf}{n}$

$\stackrel{\xaf}{n}={\displaystyle \underset{n=0}{\overset{\infty}{\sum}}n{p}_{N}\left(n|{p}_{s}\right)}=\frac{{p}_{s}}{1-{p}_{s}}=\frac{1-\lambda \delta t}{\lambda \delta t}$ . (52)

The outcome of a large number of Brownian random walks of a decaying particle, each starting from the same initial condition ${x}_{0}=0$ , leads to a distribution of time steps n given by Equation (50). One can ask, therefore, for the ensemble-average of the mean-square displacement (49)

$\langle {x}^{2}\rangle \equiv {\langle \langle {x}^{2}\left(n\delta t\right)\rangle \rangle}_{n}=2D\stackrel{\xaf}{n}\delta t=\frac{2D}{\lambda}\left(1-\lambda \delta t\right)\approx 2D/\lambda $ (53)

where, again, the term $\lambda \delta t\ll 1$ was dropped. The result (53) of the ensemble average is a root-mean square displacement $\sqrt{\langle {x}^{2}\rangle}$ equal to the mathematical diffusion length ${\zeta}_{\text{m}}$ .

It is instructive to recalculate the ensemble average of the mean-square displacement starting with the second equality in (49), in which the continuous variable t, rather than the discrete time step n, is the variable of interest. The duration t of a Brownian random walk is a realization of a continuous random variable T governed by the exponential distribution $E\left(\lambda \right)$ with PDF (see Equation (3))

${p}_{T}\left(t|\lambda \right)=\lambda {\text{e}}^{-\lambda t}$ . (54)

The distribution is normalized

${\int}_{0}^{\infty}{p}_{T}\left(t|\lambda \right)\text{d}t}=1$ (55)

with mean $\stackrel{\xaf}{t}$

$\stackrel{\xaf}{t}={\displaystyle {\int}_{0}^{\infty}t{p}_{T}\left(t|\lambda \right)\text{d}t}={\lambda}^{-1}$ . (56)

The ensemble average of the mean-square displacement (49) is then

$\langle {x}^{2}\rangle \equiv {\langle \langle {x}^{2}\left(t\right)\rangle \rangle}_{t}=2D\stackrel{\xaf}{t}=2D/\lambda $ (57)

in necessary agreement with ensemble average (53) obtained from the geometric distribution. The geometric and exponential distributions are, respectively, the discrete and continuous distributions for the class of statistical problems referred to as waiting-time problems. For continuous or discrete Markov processes that are characterized by a lack of memory, i.e. that have the same probability of success or failure at each trial, the distribution of the duration of the process must take the form of either a geometric or exponential distribution [33] .

Figure 1 shows in graphical form a computer-simulated Brownian random walk obtained by iterative solution of the stochastic differential Equation (40). The figure displays n = 500 time steps of duration $\text{d}t=1$ unit of a decaying particle with diffusion constant $D=0.5$ units and survival probability ${p}_{s}=0.99$ (and therefore ${q}_{s}=0.01$ ). The parameters were chosen to facilitate graphical display of the process. If an actual physical example were illustrated, standard metrical units would have been employed such as: [time] = seconds, $\left[D\right]={\text{cm}}^{2}/\mathrm{sec}$ , $\left[\lambda \right]={\mathrm{sec}}^{-1}$ . For example, the parameters characterizing the radioactive isotope radon-222 are: $D\approx 11\text{\hspace{0.17em}}{\text{cm}}^{2}/\mathrm{sec}$ , $\lambda =2.1\times {10}^{-6}\text{\hspace{0.17em}}{\text{sec}}^{-1}$ , and $\zeta =\sqrt{D/\lambda}\approx 229\text{\hspace{0.17em}}\text{cm}$ . From the given decay rate, it follows from Equation (5) that the half-life of the isotope is about 3.8 days. One would therefore not expect a radon-222 atom to decay within the short span of 500 seconds. For $\text{d}t=1$ sec, the decay probability for radon-222 is ${q}_{s}=\lambda \text{d}t\approx 2.1\times {10}^{-4}\%$ . However, for a measurement interval $\text{d}t=1$ hour, one has ${q}_{s}\approx 0.76\%$ , and for $\text{d}t=1$ day, ${q}_{s}\approx 18.1\%$ .

The red trace in Figure 1 shows the Brownian trajectory of a stable particle over the duration of 500 time units. The light black curves, which constitute a horizontal parabola, are plots of $\pm \sigma \left(t\right)$ , i.e. one standard deviation (49) in each direction from the point of origin ${x}_{0}=0$ . The mostly horizontal brown line traces the discrete sequence of 500 outcomes $\epsilon =\left(1,0\right)$ of a binomial random generator, scaled by a factor 5 for better visibility in the plot. A sample of 500 Bernoulli trials with probability of decay of 1% is expected to yield $\sim 5\pm 2$ decays. In the illustrated run, 4 decays occurred where the brown trace plunges sharply from 5 to 0. The solid black trace is the Brownian trajectory of the decaying particle. At each occurrence of decay, the trajectory plunges to 0 and the process is terminated. A new process then begins as represented by the portion of the red trace recommencing at the origin and continuing until the next decay, whereupon the process of termination and commencement repeats. Whereas the TPF is a continuous conditional probability function that shows how the distribution of displacements of a single particle evolves in time, provided the particle has not decayed, the numerical solutions to the stochastic differential equation permit one to see in real time the discontinuous nature of diffusion with decay. The horizontal dashed blue lines in the figure represent boundaries that will be discussed briefly in connection with Figure 2 and in greater detail in Section 3.

A point of particular interest in Figure 1, exhibited by the Brownian trajectory (red) of the stable particle, is the strikingly unequal amount of time the particle has spent in the positive half-space (representing motion to the right of the origin) compared to the negative half-space (representing motion to the left of the origin). Since the probability to move right or left at each step is the same (50%), one might have thought that a Brownian trajectory would fluctuate so as to spend about half the time on each side of the origin in conformity with some intuitive “law of averages”. That this is not the case is a known feature of a binomial or Gaussian random walk as quantified by the arcsine law [Ref [33] , pp. 80-81].

Figure 1. Computer simulation of Brownian motion with decay. Parameters: $n=500$ time steps $\text{d}t=1$ with diffusivity $D=0.5$ and survival probability ${p}_{s}=0.99$ . (a) Displacement of stable particle (solid red). (b) Displacement of unstable particle (solid black); each decay event terminates the process, which recommences with a new particle. (c) Sequence of Bernoulli trials (solid brown) each with possible outcomes $\epsilon =\left(1,0\right)$ scaled up by factor of 5 for clarity; outcome $\epsilon =0$ signifies particle decay. (d) Parabolic branches (light black) of the root-mean-square displacement $\sigma \left(t\right)=\sqrt{2Dt}$ . (e) Absorbing boundaries (dashed blue) at ${x}_{b}=\pm 15$ .

Figure 2. Computer simulation of Brownian motion with decay, employing $n=1000$ time steps; other parameters and plot labels are the same as in Figure 1 except for scaling the Bernoulli outcomes in plot (c) by a factor 10. The distributions of net displacements and first-passage times, summarized in Table 1 and Table 2, are in accord with predictions of Equations ((53) and (80)).

Moreover, there is no formal “law of averages”; the closest rigorous principle would be the law of large numbers [Ref [33] , pp. 190-191]. What a correctly formulated law of large numbers does imply is that over many repetitions of a random walk starting from the same initial conditions the ensemble of Brownian trajectories will be found in the positive half-space to approximately the same extent as in the negative half-space [Ref [18] , p. 372]. The simulated Brownian trajectories (black) of the unstable particle in Figure 1 is consistent with this implication. Upon several repetitions of the random walk by a new particle after decay of a previous one, the ensemble of resulting trajectories, obtained by downward projection of the red trace, shows a more equal balance of time between the two half spaces.

Figure 2 shows a simulation of longer duration $n=1000$ for the same values of D and ${p}_{s}$ as in Figure 1. In this simulation the trajectory (red) of the stable particle starts out at the origin, moves again into the positive half-space, but crosses the origin at around 400 and then spends approximately 61% of the total duration, in the negative half-space. A decay probability of 1% is expected to yield $\sim 10\pm 3$ decays in 1000 Bernoulli trials; the simulation in Figure 2 yielded 11 events, as shown by the horizontal brown trace (scaled by a factor 10 for clarity) with sharp drops to 0 at each decay. The corresponding partition of the continuous trajectory into the disjointed trajectories of 11 sequentially decaying particles is still unequally distributed over the two half-spaces because 1000 time steps is too short a duration, and 11 decays lead to too small an ensemble of trajectories. Simulations of $n=5000$ or more, not reproduced here, show closer conformity to the law of large numbers.

The black trajectories of Figure 2 permit confirmation of a significant feature of the Brownian motion of an unstable particle that is addressed in Section 3, viz. the question of exit time, i.e. how much time, on average, is required for the particle to reach one of the two absorbing boundaries arbitrarily set at $\mathcal{l}=\pm 15$ as marked by dashed blue lines. Since the characteristics of Brownian motion are independent of where or when the motion begins, a new Brownian trajectory starts at the origin following each particle decay, and the random walk of the particle to either boundary terminates the process. In Section 3 the expectation of the exit time is (a) derived for a decaying particle, (b) compared with the empirical result deduced from Figure 2, and (c) shown to differ markedly from the theoretical result for a stable particle.

2.4. Langevin Equation: Analytical Solution

It is possible to solve stochastic Equation (40) analytically in closed form for the displacement $X\left(t\right)$ and PDF ${p}_{X}\left(x,t\right)$ of a particle undergoing Brownian motion with decay. It is worthwhile to do so for at least two reasons. First, it is easy to misconstrue the form of Equation (40) and thereby end up with a solution that is structurally incorrect. And second, the correct solution leads to a distribution with features that physicists do not ordinarily encounter.

Consider first the misleading (or at least incomplete) way to proceed. If one regards the numerical update algorithm (41) as a sum of n independent unit normal variates multiplied by Bernoulli variates all taken to have value $\epsilon =1$ , then $X\left(t\right)$ must also be expressible in closed form as a normal random variable of zero mean because of the stability of the normal distribution. Moreover, since a normal distribution is uniquely determined by the mean and variance, it follows from relation (44) that the solution (in the limit $n\text{d}t\to t$ ) must be

$X\left(t\right)=N\left(0,\frac{2D}{\lambda}\left(1-{\text{e}}^{-\lambda t}\right)\right)$ , (58)

whereupon the corresponding PDF, generalized on the basis of space and time translational symmetry to include initial conditions $\left({x}_{0},{t}_{0}\right)$ , takes the form

${p}_{X}^{\left(\text{L}\right)}\left(x,t|{x}_{0},{t}_{0}\right)=\frac{1}{\sqrt{4\text{\pi}D\left(1-{\text{e}}^{-\lambda \left(t-{t}_{0}\right)}\right)/\lambda}}\mathrm{exp}\left(-{\left(x-{x}_{0}\right)}^{2}/4D\left(1-{\text{e}}^{-\lambda \left(t-{t}_{0}\right)}\right)/\lambda \right)$ .(59)

The superscript L distinguishes solution (59) to the Langevin stochastic equation from solution (20) to the Fokker-Planck equation repeated below for ease of comparison

${p}_{X}^{\left(\text{FP}\right)}\left(x,t|{x}_{0},{t}_{0}\right)=\frac{1}{\sqrt{4\text{\pi}D\left(t-{t}_{0}\right)}}\mathrm{exp}\left(-{\left(x-{x}_{0}\right)}^{2}/4D\left(t-{t}_{0}\right)\right){\text{e}}^{-\lambda \left(t-{t}_{0}\right)}$ (60)

and denoted by a superscript FP.

The problem with solution (58), however, is that it no longer describes a process that can be interrupted at any time by the decay of the particle. The transformation from a discontinuous to a continuous stochastic process arose by confounding the Bernoulli random variables, which can be either 1 or 0, with the specific outcomes, or variates, all taken to be 1 in the previous calculation of expectation values. Re-examining the last line in Equation (41)―i.e. before expectation values are taken―shows, together with the stability of the normal distribution, that Langevin Equation (40), expressed in terms of random variables, can be written as

$X\left(n\text{d}t\right)\equiv X\left(t\right)=N\left(0,2D\text{d}t{\Sigma}_{n}^{2}\right)$ (61)

in which

${\Sigma}_{n}^{2}={\displaystyle \underset{j=1}{\overset{n}{\prod}}{\epsilon}_{j}^{2}}+{\displaystyle \underset{j=2}{\overset{n}{\prod}}{\epsilon}_{j}^{2}}+{\displaystyle \underset{j=3}{\overset{n}{\prod}}{\epsilon}_{j}^{2}}+\cdots +{\epsilon}_{n}^{2}$ (62)

is a random variable, not an expectation value. So as not to complicate the notation unnecessarily, the symbol for a Bernoulli random variable will remain $\epsilon $ , in departure from the conventional notation to use an upper-case letter.

The question then becomes: What kind of random variable is ${\Sigma}_{n}^{2}$ and what are its properties? Details of the analysis are left to Appendix 1, but the salient points are summarized as follows. Any random variable Z is uniquely characterized by its moment-generating function (MGF)

${g}_{Z}\left(\theta \right)\equiv \langle \mathrm{exp}\left(Z\theta \right)\rangle $ (63)

(if it exists), or its characteristic function (CF)

${h}_{Z}\left(\theta \right)\equiv \langle \mathrm{exp}\left(iZ\theta \right)\rangle $ (64)

(which always exists), or its probability function (for discrete outcomes) or probability density (for continuous outcomes) [34] . The variable $\theta $ in the argument of relations (63) and (64) has no physical significance, but merely serves as a dummy variable for purposes of differentiation, after which it is set equal to 0. In the analysis of ${\Sigma}_{n}^{2}$ in Appendix 1, the MGF was used progressively, starting from the relation $\epsilon =B\left(1,{p}_{s}\right)$ , to show that

${\epsilon}_{j}^{2}={\epsilon}_{j}=B\left(1,{p}_{s}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(j=1,2,\cdots ,n\right)$ (65)

${\omega}_{j}\equiv {\displaystyle \underset{k=j}{\overset{n}{\prod}}{\epsilon}_{k}^{2}}={\displaystyle \underset{k=j}{\overset{n}{\prod}}{\epsilon}_{k}}=B\left(1,{p}_{s}^{n-j+1}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(j=1,2,\cdots ,n\right)$ (66)

and

${\Sigma}_{n}^{2}\equiv {\displaystyle \underset{j=1}{\overset{n}{\sum}}{\omega}_{j}}$ . (67)

If the set of random variables $\left\{{\omega}_{j}\right\}$ were mutually independent, the MGF of the sum in (67) could be easily calculated. However, by virtue of the defining relation (66), any pair of the $\omega $ variables, e.g. ${\omega}_{j}={\epsilon}_{j}{\epsilon}_{j+1}\cdots {\epsilon}_{k}{\epsilon}_{k+1}\cdots {\epsilon}_{n}$ and ${\omega}_{k}={\epsilon}_{k}{\epsilon}_{k+1}\cdots {\epsilon}_{n}$ , could contain some identical factors and be highly correlated. Nevertheless, although less easily done, the MGF of ${\Sigma}_{n}^{2}$ can be shown to be

${g}_{{\Sigma}_{n}^{2}}\left(\theta \right)=\frac{1-{p}_{s}+{p}_{s}^{n+1}{\text{e}}^{n\theta}\left(1-{\text{e}}^{\theta}\right)}{\left(1-{p}_{s}{\text{e}}^{\theta}\right)}$ . (68)

Although MGF (68) does not correspond to any of the tabulated random variables known to the author, all the moments of ${\Sigma}_{n}^{2}$ can be determined by differentiation

$\langle {\left({\Sigma}_{n}^{2}\right)}^{k}\rangle ={\frac{{\text{d}}^{k}{g}_{{\Sigma}^{2}}\left(\theta \right)}{\text{d}{\theta}^{k}}|}_{\theta =0}\equiv {g}_{{\Sigma}^{2}}^{\left(k\right)}\left(0\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(k=0,1,2,\cdots \right)$ , (69)

thereby characterizing ${\Sigma}_{n}^{2}$ uniquely.

For example, moments $k=0,1,2$ of ${\Sigma}_{n}^{2}$ calculated from (68) and (69) are

$\begin{array}{l}\langle {\left({\Sigma}_{n}^{2}\right)}^{0}\rangle ={g}_{{\Sigma}^{2}}^{\left(0\right)}\left(0\right)=1\\ \langle {\left({\Sigma}_{n}^{2}\right)}^{1}\rangle ={g}_{{\Sigma}^{2}}^{\left(1\right)}\left(0\right)=\frac{{p}_{s}\left(1-{p}_{s}^{n}\right)}{1-{p}_{s}}\\ \langle {\left({\Sigma}_{n}^{2}\right)}^{2}\rangle ={g}_{{\Sigma}^{2}}^{\left(2\right)}\left(0\right)=\frac{{p}_{s}\left(1+{p}_{s}\right)\left(1-{p}_{s}^{n}\right)}{{\left(1-{p}_{s}\right)}^{2}}-\frac{2n{p}_{s}^{n+1}}{1-{p}_{s}}\end{array}$ (70)

The first line of Equation (70) expresses the completeness relation for probabilities, required of the MGF. The second line reproduces expression (48). Thus, calculation of the mean-square displacement from Equation (61)

$\langle {\left(X\left(t\right)\right)}^{2}\rangle =\langle {\left(N\left(0,2D\text{d}t{\Sigma}_{n}^{2}\right)\right)}^{2}\rangle =2D\text{d}t\langle {\Sigma}_{n}^{2}\rangle =2D\text{d}t\left(\frac{{p}_{s}\left(1-{p}_{s}^{n}\right)}{1-{p}_{s}}\right)$ , (71)

leads to

$\langle {\left(X\left(t\right)\right)}^{2}\rangle =\underset{\begin{array}{l}n\to \infty \\ \text{d}t\to 0\end{array}}{\mathrm{lim}}\left[2D\text{d}t\left(\frac{\left(1-\lambda \text{d}t\right)\left(1-{\left(1-\lambda \text{d}t\right)}^{n}\right)}{\lambda \text{d}t}\right)\right]=\frac{2D}{\lambda}\left(1-{\text{e}}^{-\lambda t}\right)$ (72)

upon substitution of relation (27) for
${p}_{s}$ and transforming from discrete time steps to continuous time. The expectation value (72) is precisely the result obtained by a different procedure in (44) and justifies the form of PDF (59). The third line of (70) enables one to calculate the variance of
${\Sigma}_{n}^{2}$ and the 4^{th} moment of the displacement

$\langle {\left(X\left(t\right)\right)}^{4}\rangle =\underset{\begin{array}{l}n\to \infty \\ \text{d}t\to 0\end{array}}{\mathrm{lim}}\left[4{D}^{2}{\left(\text{d}t\right)}^{2}\langle {\left({\Sigma}_{n}^{2}\right)}^{2}\rangle \right]=\frac{8{D}^{2}}{{\lambda}^{2}}\left(1-\left(1+\lambda t\right){\text{e}}^{-\lambda t}\right)$ (73)

by means of the same substitution and limiting process employed in Equation (72).

In short, therefore, the solution (61) takes the form of a normal distribution whose variance is itself a random variable of un-named variety (as far as the author is aware), but which is completely and uniquely specified by its MGF (68). In principle, the PDF of the distribution of ${\Sigma}_{n}^{2}$ can be calculated by taking the inverse Fourier transform of the CF; the CF itself is obtained simply by substituting $i\theta $ for $\theta $ in the MGF (68). For the purposes of this paper, however, the PDF of ${\Sigma}_{n}^{2}$ is not required.

Solution (61), in contrast to solution (58), incorporates the Bernoulli processes that generate particle decay. If, in a simulation of Brownian motion to be implemented for n time steps from Equation (61), the Bernoulli variate ${\epsilon}_{k}=0$ at time step $k\le n$ , it follows from Equation (66) that all the variates ${\omega}_{j}=0$ for $j=1,2,\cdots ,k$ , and therefore ${\Sigma}_{k}^{2}=0$ , and $X\left(t=k\text{d}t\right)=N\left(0,0\right)=0$ . (Note: $N\left(0,0\right)=0$ signifies that $\underset{\sigma \to 0}{\mathrm{lim}}\left({\sigma}^{-1}\mathrm{exp}\left(-{x}^{2}/2{\sigma}^{2}\right)\right)=0$ of the Gaussian PDF.) Thus, the random walk of the decayed particle has been realistically terminated at the randomly selected ${k}^{\text{th}}$ time step. With regard to notation, the subscript on ${\Sigma}_{n}^{2}$ can be taken to represent the full length (i.e. number of time steps) of a simulated random walk, and not (as before) a pre-determined number of sequential steps survived by the particle. In the limit of an infinitely long random walk, it follows from Equation (72) that the (infinitely) many particle trajectories arising from the introduction of a new particle after decay of each previous one, yield a root-mean-square (RMS) displacement equal to the mathematical diffusion length $\sqrt{2D/\lambda}$ , in accord with the ensemble averages (53) and (57).

An empirical test of the ensemble-averaged root-mean-square (RMS) displacement (53) can be made using the computer-generated Brownian trajectories of Figure 2; the relevant data are recorded in Table 1. Columns 1 and 2 list

Table 1. Test of Theoretical Ensemble-Averaged RMS Displacement Using Simulated Trajectories of Figure 2.

chronologically the decay events and corresponding times of decay. Columns 3 and 4 show the net displacement (to 4 significant figures) and square of displacement (truncated to 2 significant figures) between the time of creation of the particle at the origin ${x}_{0}=0$ up to the time step just prior to its decay. As shown, the mean of the 11 RMS values is very close to the theoretical prediction calculated from Equation (53) with the parameters used in the simulations of Figure 2: $D=0.5$ , $\text{d}t=1$ , ${p}_{s}=0.99$ , in which case $\lambda =\left(1-{p}_{s}\right)/\text{d}t=0.01$ .

Now that the stochastically correct solution (61) to the Langevin equation has been derived and shown to justify PDF (59), it is informative to compare the latter with PDF (60) derived from the Fokker-Planck Equation (18). Figure 3 and Figure 4 respectively show plots of the Langevin and Fokker-Planck PDFs as a function of displacement for different random walk durations. The principal feature to notice is that, as the duration of the random walk increases, the PDF derived from the Langevin solution approaches a steady-state Gaussian distribution of mean 0, maximum ${\left(4\text{\pi}D/\lambda \right)}^{-1/2}$ , and variance $2D/\lambda $ , whereas the distribution derived from the Fokker-Planck solution vanishes, i.e. approaches a Gaussian distribution of mean 0, maximum 0, and infinite width. In light of the previous discussion, these differences are understandable as follows:

・ The Fokker-Planck Equation (18), although it includes the particle decay rate, describes one long continuous process, i.e. the evolution of the probability of a particular diffusing particle to be found at an arbitrary location x at time t, provided the particle survived to time t. As t increases, the survival probability of that particular particle decreases exponentially as ${\text{e}}^{-\lambda t}$ , and the mean- square displacement of its Brownian motion spreads as a power law $2Dt$ . However, only after an infinitely long time is the probability for the particle to reach any location precisely zero.

・ In contrast to the preceding, the Langevin Equation (40) describes a potentially infinite number of independent Brownian trajectories, disconnected one from the other by events of particle decay. As t increases, the ensemble of independent trajectories come from particles that have survived for different lengths of time. In the limit of an infinitely long time, the probability density (59) does not vanish, but describes, instead, the ensemble-averaged statistics characterized by a Gaussian distribution with mean-square displacement $2D/\lambda $ .

Figure 3. Spatial variation of the probability density function (59) obtained from solution of the Langevin stochastic Equation (40) for time units t = 1 (red), 2 (gold), 5 (green), 10 (blue), infinite (violet). Process parameters are initial location ${x}_{0}=0$ , diffusivity $D=0.5$ , decay rate $\lambda =0.1$ .

Figure 4. Spatial variation of the probability density function (60) obtained from solution of the Fokker-Planck Equation (18) for the same time units, process parameters, and color code as in Figure 3.

Figure 5. Temporal variation of the probability density function (59) of the Langevin Equation (red) and (60) of the Fokker-Planck Equation (black) for spatial coordinates x = 0 (dash), 1 (solid), 2 (dot-dash). Process parameters are diffusivity $D=0.5$ , decay rate $\lambda =0.1$ .

Further perspective on the differences between the two approaches (Langevin vs Fokker-Planck) is given by Figure 5, which shows plots of PDFs (59) and (60) as a function of time for several different locations. The probability for a particle to remain at the origin $\left({x}_{0}=0\right)$ decreases monotonically in both cases, but asymptotically approaches 0 in the Fokker-Planck PDF and ${\left(4\text{\pi}D/\lambda \right)}^{-1/2}$ in the Langevin PDF. The probability for a particle to remain at a location away from the origin is 0 to begin with, due to the imposed boundary condition (19) which applies to both PDFs, rises to a maximum, and subsequently decreases to the preceding asymptotic limits.

It is to be emphasized that the Fokker-Planck and Langevin descriptions are both valid. It is not that one is right and the other wrong; rather, the two analytical methods provide different perspectives on the process of Brownian motion with decay. The Fokker-Planck equation describes a continuous process; it gives a statistical description of the displacement of a decaying particle for as long as the particle might survive in the course of an infinite time span. The Langevin equation describes a sequence of discontinuous processes; it gives a statistical description of the displacement of an ensemble of particles up to the instant each actually fails to survive. It is understandable, therefore, why the theoretical mean-square displacements obtained by the two approaches are different. A significant point of this paper, however, is that for stable particles the two approaches lead to identical results because both would then be describing a single continuous process of infinite duration.

3. First-Passage Times

The horizontal dashed blue lines in Figure 1 and Figure 2 visually address an important question that arises in the measurement of radioactive atoms and other unstable particles undergoing Brownian motion: At what time will the particle, starting from a given point (e.g. ${x}_{0}=0$ ), first reach some other specified location? This is the problem of first-passage time (FPT) or exit time. If, upon reaching the designated location, the particle is removed from the system, the location is said to be an absorbing boundary. At an absorbing boundary, the probability density vanishes. Other kinds of boundaries can be transmitting, i.e. have no effect on the particle’s subsequent motion, or reflecting, in which case the particle current density, but not probability density, vanishes. This paper is concerned with absorbing boundaries, since these conditions characterize a variety of measurement protocols employing passive diffusion of radioactive atoms in gas, aerosol, and dust, as well as unstable molecules in a fluid medium..

In Figure 1 and Figure 2, two absorbing boundaries were arbitrarily set at ${x}_{b}=\pm 15$ units from the origin. As shown in Figure 1, the unstable particle (black) first reached boundary ${x}_{b}=+15$ at approximately ${t}_{b}=125$ before decaying at ${t}_{d2}=250$ . Note, however, that a prior decay occurred close to ${t}_{d1}=33$ . Therefore the first particle did not reach the boundary before decaying, and the second particle, which began a new process, reached the boundary after a diffusion time of ${t}_{b}-{t}_{d1}=92$ . In the simulation shown in Figure 2, there are four decays, the last occurring at time step ${t}_{d4}=99$ , before the fifth unstable particle reached boundary ${x}_{b}=+15$ for the first time at time step ${t}_{b}=158$ , and then decayed shortly afterward at time step ${t}_{d5}=163$ . The FPT for particle 5 is therefore ${t}_{b}-{t}_{d4}=59$ . In a FPT problem with unstable particles and absorbing boundaries, the process of Brownian motion is terminated by decay or by reaching a boundary; there is no second-passage time. Table 2, which will be used shortly to test an important theoretical result, summarizes all the first-passage times of the trajectories in Figure 2.

FPT problems have been studied in great detail for non-decaying diffusers; see, for example [23] . In this paper, however, the problem of FPT is solved for unstable diffusers. Since the FPT is a random variable, for which the symbol T will again be used, solving the FPT problem ordinarily means calculating the expectation value ${T}_{\text{FPT}}$

${T}_{\text{FPT}}\equiv \langle T\rangle ={\displaystyle {\int}_{0}^{\infty}t{p}_{T}\left(t\right)\text{d}t}={\displaystyle {\int}_{0}^{\infty}\left(1-{F}_{T}\left(t\right)\right)\text{d}t}$ (74)

in which ${p}_{T}\left(t\right)$ is the probability density function (PDF) of T and ${F}_{T}\left(t\right)$ is the cumulative probability function (CPF)

${F}_{T}\left(t\right)={\displaystyle {\int}_{0}^{t}{p}_{T}\left({t}^{\prime}\right)\text{d}{t}^{\prime}}$ . (75)

A proof of the second equality in relation (74) for random variables such as T defined on the non-negative real numbers is given in Appendix 2 of [32] .

In general, two mathematically different approaches have been used to find ${T}_{\text{FPT}}$ for a stable particle. One approach involves calculation of ${F}_{T}\left(t\right)$ , starting with the solution to the Fokker-Planck equation. The other approach yields ${T}_{\text{FPT}}$ directly as the solution to a Poisson equation. These two methods are generalized in the following two sub-sections so as to apply to Brownian motion with decay. In the example analyzed, the particle is absorbed upon reaching either of the symmetrically located boundaries ${x}_{b}=\pm \mathcal{l}$ before decaying. With the methods developed below the boundary conditions can be easily modified to apply to other systems.

3.1. Method 1: Calculation of T_{FPT} by Image Functions

The problem is to find a PDF ${p}_{x}\left(x,t\right)$ of the diffusing particle with initial condition (19) for ${x}_{0}=0$

${p}_{x}\left(x,0\right)=\delta \left(x\right)$ (76)

that satisfies boundary conditions ${p}_{x}\left(\mathcal{l},t\right)={p}_{x}\left(-\mathcal{l},t\right)=0$ . Since the Fokker- Planck Equation (18) is linear, the desired PDF can be constructed from solution (20) by the method of images [35] , and leads to the infinite sum

$\begin{array}{l}{p}_{X}\left(x,t\right)=\underset{N\to \infty}{\mathrm{lim}}{\displaystyle \sum _{n=-N}^{N}{p}_{X}\left(x,t,n\right)}\\ =\underset{N\to \infty}{\mathrm{lim}}{\displaystyle \sum _{n=-N}^{N}\frac{1}{\sqrt{4\text{\pi}Dt}}\left(\mathrm{exp}\left(-\frac{{\left(x+4n\mathcal{l}\right)}^{2}}{4Dt}\right)-\mathrm{exp}\left(-\frac{{\left(x+\left(4n-2\right)\mathcal{l}\right)}^{2}}{4Dt}\right)\right){\text{e}}^{-\lambda t}}\end{array}$ (77)

in which the second line of (77) defines the function ${p}_{X}\left(x,t,n\right)$ .

A sense of the structure of solution (77) is provided by Figures 6-10. The figures depict the probability density of an unstable particle with diffusivity $D=1.0\text{\hspace{0.17em}}{\text{cm}}^{2}\cdot {\text{s}}^{-1}$ and decay rate $\lambda =2.0\times {10}^{-5}\text{\hspace{0.17em}}{\text{s}}^{-1}$ undergoing Brownian motion between absorbing barriers at $\mathcal{l}=1\text{\hspace{0.17em}}\text{m}$ to the left and right of the origin. In each figure the dashed black plot is the Gaussian PDF ${p}_{X}^{\left(\text{FP}\right)}\left(x,t\right)$ , Equation (60), of the unconstrained particle at stated time t; the solid red plot is the specified linear superposition of functions ${p}_{X}\left(x,t,n\right)$ contributing to solution (77); the vertical dashed blue lines mark the boundaries, and the horizontal solid blue line marks the physical space along the axis within which the particle is required to be confined.

Figure 6 shows the PDFs ${p}_{X}^{\left(\text{FP}\right)}\left(x,t\right)$ and ${p}_{X}\left(x,t,0\right)$ at time $t=4000\text{\hspace{0.17em}}\text{s}$ . At this time, there is a significant probability that the unconstrained particle (black), described by ${p}_{X}^{\left(\text{FP}\right)}\left(x,t\right)$ , could be found outside the designated physical region. However, the $N=0$ component ${p}_{X}\left(x,t,0\right)$ , which is a superposition of ${p}_{X}^{\left(\text{FP}\right)}\left(x,t\right)$ and a negative Gaussian image centered on $2\mathcal{l}$ , satisfies the right- side boundary condition ${p}_{x}\left(\mathcal{l},t\right)=0$ , although it does not satisfy the boundary condition on the left side. Figure 7 shows that the $N=1$ component with superposition of ${p}_{X}\left(x,t,-1\right)$ , ${p}_{X}\left(x,t,0\right)$ , ${p}_{X}\left(x,t,+1\right)$ , which include Gaussian images centered on $-4\mathcal{l}$ , $-2\mathcal{l}$ , $4\mathcal{l}$ , $6\mathcal{l}$ , satisfies both boundary conditions at $t=4000\text{\hspace{0.17em}}\text{s}$ . The figure depicts the entire waveform of the superposition, but only the portion of the probability density (red) between the two boundaries describes the particle in the physically allowed region of Brownian motion.

As time increases, the particle can diffuse to greater distances from the origin, and the probability density functions contributing to PDF (77) spread. The $N=1$ solution that satisfied the boundary conditions at $t=4000\text{\hspace{0.17em}}\text{s}$ , no longer satisfies these conditions at the later time $t=20000\text{\hspace{0.17em}}\text{s}$ , as shown in Figure 8 and Figure 9. Figure 8 shows that the $N=0$ solution satisfies just the right boundary condition; Figure 9 shows that the $N=1$ solution still does not satisfy the

Figure 6. Spatial variation (red) of the $N=0$ component to ${p}_{X}\left(x,t\right)$ , Equation (77), with two absorbing boundaries (dashed blue lines at $\pm 1$ ) and parameters $t=4000$ , $D=1.0\times {10}^{-4}$ , $\lambda =2.0\times {10}^{-5}$ . Superposed is the probability density (dashed black) of the unconstrained particle initially at ${x}_{0}=0$ . The solid blue line marks the physical region between intended boundaries. The particle is confined only on the right side: ${p}_{X}\left(+1,t\right)=0$ . Numerical values are in MKS units.

Figure 7. Spatial variation of the superposition of $N=0,\pm 1$ components in solution (77) for ${p}_{X}\left(x,t\right)$ at $t=4000$ . Parameters and color codes are the same as in Figure 6. The particle is confined on both sides: ${p}_{X}\left(\pm 1,t\right)=0$ .

Figure 8. Spatial variation of the $N=0$ component to ${p}_{X}\left(x,t\right)$ , Equation (77), at $t=20\uff0c000$ . Other parameters and color codes are the same as in Figure 6. The particle is again confined only on the right side.

Figure 9. Spatial variation of the $N=0,\pm 1$ components to ${p}_{X}\left(x,t\right)$ , Equation (77), at $t=20\uff0c000$ . Parameters and color codes are the same as in Figure 6. The particle is still confined only on the right side.

left boundary condition. In Figure 10 it is seen that truncation of PDF (77) at $N=2$ , which includes image functions centered on $-8\mathcal{l}$ , $-6\mathcal{l}$ , $-4\mathcal{l}$ , $-2\mathcal{l}$ , $2\mathcal{l}$ , $4\mathcal{l}$ , $6\mathcal{l}$ , $8\mathcal{l}$ , $10\mathcal{l}$ , does confine the particle between the two boundaries

Figure 10. Spatial variation of the $N=0,\pm 1,\pm 2$ components to ${p}_{X}\left(x,t\right)$ , Equation (77), at $t=20\uff0c000$ . Parameters and color codes are the same as in Figure 6. The particle is now confined on both sides.

as required. However, since time extends over an infinite range $\left(\infty \ge t\ge 0\right)$ , the necessity for an infinite number of image functions in (77) becomes apparent.

The integral of PDF (77) over the physically allowed region defines the survival function

$\begin{array}{l}S\left(t\right)={\displaystyle {\int}_{-\mathcal{l}}^{\mathcal{l}}{p}_{X}\left(x,t\right)\text{d}x}\\ ={\displaystyle \underset{n=-\infty}{\overset{\infty}{\sum}}\frac{1}{2}\left[\text{erf}\left(\frac{\left(4n+1\right)\mathcal{l}}{\sqrt{4Dt}}\right)+\text{erf}\left(\frac{\left(4n-3\right)\mathcal{l}}{\sqrt{4Dt}}\right)-2\text{erf}\left(\frac{\left(4n-1\right)\mathcal{l}}{\sqrt{4Dt}}\right)\right]{\text{e}}^{-\lambda t}}\end{array}$ (78)

which is the probability that the particle has not reached either boundary at time t. The probability $S\left(t\right)$ that the particle remains in the interval $\left[-\mathcal{l},\mathcal{l}\right]$ after time t means that $T>t$ , or, in terms of the definition (75) of cumulative probability, one has

$S\left(t\right)=\mathrm{Pr}\left(T>t\right)=1-{F}_{T}\left(t\right)$ . (79)

It then follows from relation (74) that the expectation value ${T}_{\text{FPT}}$ is given by the integral of the survival function over time

${T}_{\text{FPT}}={\displaystyle {\int}_{0}^{\infty}S\left(t\right)\text{d}t}=\frac{{\left({\text{e}}^{\mathcal{l}/\zeta}-1\right)}^{2}}{\lambda \left({\text{e}}^{2\mathcal{l}/\zeta}+1\right)}$ (80)

where $\zeta $ , the characteristic diffusion length, is defined by relation (46).

Figure 2 provides an empirical test of the theoretical expression (80) with the relevant data recorded in Table 2. The first two columns display chronologically

Table 2. Empirical Test of Theoretical FPT Using Simulated Brownian Trajectories of Figure 2.

the instances of particle decay and the time step at which decay occurred. The third column shows the time step at which a particle first reached or exceeded the right ( $\mathcal{l}=+15$ ) or left ( $\mathcal{l}=-15$ ) absorbing boundary. A dash signifies that the particle decayed before reaching either boundary. The fourth column is the difference between columns 3 and 2, which gives the time interval measured from the point at which the particle was placed at the origin $x=0$ . As shown, the mean of the five FPT values is very close to the theoretical prediction 76.4 calculated from Eq. (80) with the parameters used in the simulations shown in Figure 2: $D=0.5$ , $\lambda =\left(1-{p}_{s}\right)/\text{d}t=0.01$ in which ${p}_{s}=0.99$ and $\text{d}t=1$ . For comparison, the theoretical FPT for a particle that does not decay is derived in the following section (Equation (92)) and reduces to ${\mathcal{l}}^{2}/2D=225$ for the case at hand.

Although the PDF ${p}_{T}\left(t\right)$ of the random variable T was not required to derive the mean first-passage time (80), it can be obtained, if needed, by differentiating the survival function

${p}_{T}\left(t\right)=\text{d}{F}_{T}\left(t\right)/\text{d}t=-\text{d}S\left(t\right)/\text{d}t$ . (81)

The resulting expression, however, is complicated and not needed in this paper.

3.2. Method 2: Calculation of T_{FPT} from a Screened Poisson Equation

The method employed in Section 3.1 yielded the CPF ${F}_{T}\left(t\right)$ (or, equivalently, the survival function $S\left(t\right)$ ) from which all statistical properties of the FPT can be calculated. However, if all one wanted was ${T}_{\text{FPT}}$ , it can be derived directly from a master equation by a method that has been employed for stable particles, such as the diffusion of molecules of biological significance [36] . This method must be generalized to account for the finite lifetime of decaying particles.

Consider an unstable particle in one-dimensional Brownian motion with time steps of $\delta t$ , on a lattice with coordinate spacing $\delta x$ and absorbing boundaries at ${x}_{b}=\pm \mathcal{l}$ . The probability to step either left or right is equal to 1/2, and the particle does not jump over any lattice points. If $T\left(x\right)$ is the time to reach either boundary from point x, it must satisfy the following discrete equation

$T\left(x\right)=\delta t+\frac{1}{2}\left(T\left(x+\delta x\right)+T\left(x-\delta x\right)\right)\left(1-\lambda \delta t\right)$ (82)

because: a) point x can be reached only from points $x+\delta x$ and $x-\delta x$ , b) the probability of a transition from these two points is the same (1/2), and c) the transition to x during interval $\delta t$ can be made only if the particle has survived the transition. Recall from Equation (27) that the survival probably per step is ${p}_{s}=1-\lambda \delta t$ .

Rearrangement of Equation (82) to the form

$\begin{array}{l}\frac{1}{2}\left[\left(T\left(x+\delta x\right)-T\left(x\right)\right)-\left(T\left(x\right)-T\left(x-\delta x\right)\right)\right]\\ =-\delta t+\lambda \delta t\frac{1}{2}\left[T\left(x+\delta x\right)+T\left(x-\delta x\right)\right]\end{array}$ (83)

which is expressible as a second derivative

$\frac{{\text{d}}^{2}T\left(x\right)}{\text{d}{x}^{2}}{\left(\delta x\right)}^{2}=-2\delta t\left(1-\lambda T\left(x\right)\right)$ , (84)

leads in the limit $\delta x\to 0$ and $\delta t\to 0$ to a screened Poisson equation, encountered in the physics of ionized gases [37] [38] ,

$D\frac{{\text{d}}^{2}T\left(x\right)}{\text{d}{x}^{2}}-\lambda T\left(x\right)=-1$ (85)

in which the ratio ${\left(\delta x\right)}^{2}/2\delta t$ is again taken to be the diffusion coefficient D as defined in Equation (32). In the absence of the term proportional to $\lambda $ , Equation (85) takes the standard form of Poisson’s equation [39] .

Although a derivation will not be given here, one can also arrive in several steps at Equation (85) by integrating the backward Fokker-Planck equation

$\frac{\partial p\left(x,t|{x}_{0},{t}_{0}\right)}{\partial {t}_{0}}=-D\frac{{\partial}^{2}p\left(x,t|{x}_{0},{t}_{0}\right)}{\partial {x}_{0}{}^{2}}+\lambda p\left(x,t|{x}_{0},{t}_{0}\right)$ (86)

for which Equation (18) is the associated forward Fokker-Planck equation. It is important to note, however, that the coordinate x appearing in Equation (85) and in all expressions involving $T\left(x\right)$ refers to the initial point from which the particle diffuses to a boundary, and therefore actually corresponds to coordinate ${x}_{0}$ in Equation (86). Where there is no confusion, it is standard notation to use the variable x rather than ${x}_{0}$ as the argument of the FPT Equation (85) and solution.

The solution to Equation (85) with implementation of boundary conditions $T\left(\pm \mathcal{l}\right)=0$ is

$T\left(x\right)=\frac{1}{\lambda}\left(1-\frac{\mathrm{cosh}\left(x/\zeta \right)}{\mathrm{cosh}\left(\mathcal{l}/\zeta \right)}\right)$ . (87)

Comparison with solution (80) of Section 3.1 can be made by setting the initial location $x=0$ , whereupon Equation (87) reduces to

$T\left(0\right)=\frac{1}{\lambda}\left(1-\frac{1}{\mathrm{cosh}\left(\mathcal{l}/\zeta \right)}\right)=\frac{{\left({\text{e}}^{\mathcal{l}/\zeta}-1\right)}^{2}}{\lambda \left({\text{e}}^{2\mathcal{l}/\zeta}+1\right)}={T}_{\text{FPT}}$ . (88)

Methods 1 and 2, although very different, lead to the same final result, as they must for consistency.

3.3. Critical Role of Particle Decay

The Fokker-Planck Equation (18) incorporates particle instability by means of a decay term $-\lambda p\left(x,t|{x}_{0},{t}_{0}\right)$ that affects the resulting transition probability density (20) only through a global exponential factor ${\text{e}}^{-\lambda t}$ . Thus, letting $\lambda $ approach 0 smoothly generates the probability density of the unconstrained stable particle. Likewise, the vanishing of $\lambda $ in the survival function (78) smoothly generates the survival function of a stable particle confined between absorbing boundaries. This ostensible continuity between decay and stability gives a misleading impression of the effect of particle decay on the FPT problem.

To illustrate the radical change in outcome that can be engendered by the decay process, consider again by method 2 the FPT problem of a decaying particle in Brownian motion between two non-symmetrically placed absorbing boundaries $b>a$ . The solution to Equation (85) with boundary conditions $T\left(a\right)=T\left(b\right)=0$ is

${T}_{ab}\left(x\right)=\frac{1}{\lambda}\left[\left(\frac{{\text{e}}^{-a/\zeta}-{\text{e}}^{-b/\zeta}}{{\text{e}}^{\left(a-b\right)/\zeta}-{\text{e}}^{\left(b-a\right)/\zeta}}\right){\text{e}}^{x/\zeta}+\left(\frac{{\text{e}}^{b/\zeta}-{\text{e}}^{a/\zeta}}{{\text{e}}^{\left(a-b\right)/\zeta}-{\text{e}}^{\left(b-a\right)/\zeta}}\right){\text{e}}^{-x/\zeta}+1\right]$ . (89)

In the limit $a\to -\infty $ , there is only a single boundary b on the right, and the FPT (89) reduces in this limit to

${T}_{b}\left(x\right)=\frac{1}{\lambda}\left(1-{\text{e}}^{-\sqrt{\lambda /D}\left(b-x\right)}\right)$ (90)

in which the $\lambda $ dependence is explicitly shown. Note that ${T}_{b}\left(x\right)$ (90) is a finite quantity although the decaying particle can be located at any point in the infinite range to the left of b. In the limit that $b\to \infty $ as well, the particle is free to walk the entire x-axis, and the FPT reduces to the statistical lifetime (4) ${\lambda}^{-1}$ as expected.

For a stable particle $\left(\lambda =0\right)$ , however, the limit of the FPT (90) is infinite no matter how close to b the particle is located initially. This is evident from the leading term in the series expansion of ${T}_{b}\left(x\right)$

${T}_{b}\left(x\right)\sim \frac{b-x}{\sqrt{D\lambda}}-\frac{{\left(b-x\right)}^{2}}{2D}+\frac{{\left(b-x\right)}^{3}}{6}\sqrt{\frac{\lambda}{{D}^{3}}}+\cdots $ . (91)

The first term of (91), in which the denominator is the characteristic diffusion velocity [Ref [15] ], is singular as $\lambda \to 0$ . Physically, even if the particle is close to the boundary, there is an infinite number of random paths that the particle can take to arrive at point b. In the case of two finite boundaries, ${T}_{ab}\left(x\right)$ remains finite as $\lambda \to 0$

$\underset{\lambda \to 0}{\mathrm{lim}}{T}_{ab}\left(x\right)\to \frac{\left(b-x\right)\left(x-a\right)}{2D}$ . (92)

For a particle initially located at $x=0$ , the left boundary $-a>0$ , and therefore $\underset{\lambda \to 0}{\mathrm{lim}}{T}_{ab}\left(0\right)\to \left|ab\right|/2D$ .

Further insight is gained by examining the probability density of the FPT in the case of a single absorbing boundary at $b>0$ and a particle initially located at 0. ${p}_{T}\left(t\right)$ is obtained by method 1 of the preceding section as applied to PDF (77) truncated to include only the component ${p}_{X}\left(x,t,0\right)$ , since a Gaussian centered on 0 and a negative Gaussian image centered on $2b$ suffice to maintain ${p}_{X}\left(b,t\right)=0$ for all time t. The calculation yields

${p}_{T}\left(t\right)=\left[\frac{b\mathrm{exp}\left(-{b}^{2}/4Dt\right)}{\sqrt{4\text{\pi}D{t}^{3}}}+\lambda \text{erf}\left(\frac{b}{\sqrt{4Dt}}\right)\right]{\text{e}}^{-\lambda t}$ . (93)

The first term in brackets, known as the Smirnov density [40] or an inverse Gaussian density [41] , is independent of $\lambda $ and characterized by a heavy tail, i.e. an asymptotic power law ${t}^{-3/2}$ . It is this density that leads to an infinite first moment ${T}_{\text{FPT}}$ for a stable particle

${\langle T\rangle |}_{\lambda =0}={\displaystyle {\int}_{0}^{\infty}\frac{bt\mathrm{exp}\left(-{b}^{2}/4Dt\right)}{\sqrt{4\text{\pi}D{t}^{3}}}\text{d}t}=\infty $ . (94)

However, the Smirnov density multiplied by an exponential ${\text{e}}^{-\lambda t}$ yields a finite first moment

$\int}_{0}^{\infty}\frac{bt\mathrm{exp}\left(-{b}^{2}/4Dt\right){\text{e}}^{-\lambda t}}{\sqrt{4\text{\pi}D{t}^{3}}}\text{d}t}=\frac{b{\text{e}}^{-b/\zeta}}{\sqrt{4D\lambda$ (95)

because the exponential function decreases faster than a power law over the infinite extent of the tail, as shown in Figure 11. The insert in the figure extends the time axis far into the region of the tail. The contribution to the FPT of the second term in (93) is

${\int}_{0}^{\infty}\lambda t\text{erf}\left(\frac{b}{\sqrt{4Dt}}\right){\text{e}}^{-\lambda t}\text{d}t}=\frac{1}{\lambda}\left[1+\left(\frac{b+2\zeta}{2\zeta}\right){\text{e}}^{-b/\zeta}\right]$ ,

which, together with expression (95), sum to ${\lambda}^{-1}\left(1-{\text{e}}^{-b/\zeta}\right)$ , in agreement with Equation (90) for ${T}_{b}\left(0\right)$ obtained directly from the screened Poisson Equation (85).

Figure 11. Temporal variation of the probability density ${p}_{T}\left(t\right)$ of first-passage time to a single absorbing boundary at ${x}_{b}=10$ . Parameters are diffusion coefficient $D=0.5$ and decay rate $\lambda =$ (a) 0 (red), (b) $5\times {10}^{-3}$ (blue), (c) $10\times {10}^{-3}$ (black). The insert shows the faster decrease of ${p}_{T}\left(t\right)$ compared to the asymptotic power law of an inverse Gaussian (red).

4. Validity of the Stochastic Model

This paper is concerned principally with effects of particle instability (decay) on Brownian motion. Therefore, to have included a frictional force in the Langevin Equation (LE) or a drift term in the Fokker-Planck Equation (FPE) would have unnecessarily complicated the problem, increased the length of the analysis, and detracted from the primary objective.

Nevertheless, the physical cause of fluctuations in Brownian motion is in fact ascribable to random impacts on the observed particle by collisions with other particles of the medium. In the diffusion of radon gas in air, for example, a particular radon atom (or a particular dust particle to which is attached a radioactive polonium ion from a radon decay) is buffeted by impacts from ambient oxygen and nitrogen molecules. In keeping with the fluctuation-dissipation theorem [42] , impacts by particles of the medium are responsible not only for Brownian motion (i.e. displacement fluctuations) but also for the friction or viscosity of the medium. Under what circumstances, therefore, is it legitimate to ignore the dissipation term in the LE?

^{1}Mobility
${\mu}_{p}$ is the proportionality factor in
$V={\mu}_{p}F$ relating velocity V and dissipative force F.

Since neglect of friction in the analysis of Brownian motion is the essence of the approach taken by Einstein [Ref. [1] ], the question has been examined thoroughly. In brief, neglect of friction is justified provided the time interval
$\tau \equiv m{\mu}_{p}$ , in which m is the particle mass and
${\mu}_{p}$ is the particle mobility^{1}, is large enough that the displacement
$x\left(t+\tau \right)$ is independent of the displacement
$x\left(t\right)$ , but small compared to the time interval
$\Delta \text{\hspace{0.05em}}t$ between observations [43] . This criterion leads to an inequality [Ref. [8] , pp. 66-67]

$\Delta t>mD/{k}_{\text{B}}{T}_{\text{e}}$ (96)

in which D is the particle diffusion coefficient, ${k}_{\text{B}}$ is Boltzmann’s constant, and ${T}_{\text{e}}$ is the equilibrium temperature. For a radon-222 atom diffusing at room temperature ${T}_{\text{e}}\approx 300\text{\hspace{0.17em}}\text{K}$ , Equation (96) becomes $\Delta t>1\text{\hspace{0.17em}}\text{ns}$ , which is readily satisfied in most experiments or measurement protocols.

Another issue that also has arisen in the past is the validity of using Fick’s law (14) to model diffusion. Serber [44] has shown in the context of neutron diffusion that Fick’s law should be applicable if the net neutron flux across a surface is small compared to the flux in either direction. As shown in Appendix 2, application of this argument to radioactive particles of statistical lifetime ${\lambda}^{-1}$ leads to an inequality analogous to that of (96)

${\lambda}^{-1}\gg 4mD/3{k}_{\text{B}}{T}_{\text{e}}$ . (97)

Using radon-222 at room temperature as an example, one has ${\lambda}_{\text{Rn-222}}^{-1}\approx 4.8\times {10}^{5}\text{\hspace{0.17em}}\text{s}\gg 1.3\text{\hspace{0.17em}}\text{ns}$ , which is very well satisfied. In fact, relation (97) is well satisfied even for the short-lived progeny of radon-222, such as polonium- 214 for which ${\lambda}_{\text{Po-214}}^{-1}\approx 3.4\times {10}^{-4}\text{\hspace{0.17em}}\text{s}\gg 1.3\text{\hspace{0.17em}}\text{ns}$ One can conclude, therefore, that the stochastic model in this paper is valid for the Brownian motion of radioactive atoms in air probed at time intervals in excess of a few nanoseconds. Moreover, Fick’s law as applied to radon diffusion has been confirmed experimentally [45] .

5. Conclusions

Motivated by new experimental methods to measure radioactive atoms and ions diffusing in gas, on dust, or in liquids, this paper analyzed the Brownian motion of decaying particles, a process that has received relatively little prior attention. In particular, equations to be interpreted as the Fokker-Planck and Langevin equations of an unstable particle were derived and solved. Also, the equations of time of first passage of unstable particles to absorbing boundaries were derived and solved.

The phenomenon of particle decay introduces into Brownian motion an additional time parameter (the decay rate $\lambda $ or statistical lifetime ${\lambda}^{-1}$ ) that leads to marked differences in the analysis of unstable, compared to stable, diffusing particles. Whereas Brownian motion of a stable particle is a continuous Markov process that can in principle be followed for an infinite length of time, the process terminates abruptly at the decay of the unstable particle. A mathematical consequence of this discontinuity is reflected in the different analytical content of the Fokker-Planck Equation (FPE), which yields a transition probability density $p\left(x,t|{x}_{0},{t}_{0}\right)$ , and the Langevin Equation (LE), which yields the distribution function and trajectories of the corresponding process variable $X\left(t\right)$ . For stable particles, the FPE and LE both describe a continuous Wiener process and provide entirely equivalent information. For decaying particles, however, the former (FPE) describes the probability of displacement of a single particle throughout the course of its existence, which terminates with 100% probability only after an infinite time interval. In contrast, the latter (LE) gives an ensemble statistical description of the discontinuous trajectories of numerous sequential particles that have undergone Brownian motion up to the instant of their actual, randomly occurring decays. The two approaches are both valid, but provide different, and differently interpreted, statistical results relating to means, variances, and higher moments.

In regard to the statistical description of the stochastic process, there is a critical difference in mathematical structure of the LE and its solution for an unstable particle compared to the LE and solution of a stable particle. Assuming that Brownian motion of a stable particle is modeled by a frictionless Wiener process with constant diffusivity D, the resulting solution is a normal distribution of zero mean and variance $2Dt$ . However, the LE for the decaying particle entails a mixture of Wiener processes for displacement and Bernoulli processes for decay and leads to a solution in the form of a normal distribution of zero mean but with a variance that is itself a random variable. This random variable, ${\Sigma}_{n}^{2}$ , although not among any known to the author, depends on the number of time steps n and survival probability ${p}_{s}$ per step and is completely characterized by its moment-generating function. It is shown in Appendix 1 that in the limit of an infinite number of time steps, ${\Sigma}_{\infty}^{2}$ behaves increasingly like an exponential distribution as ${p}_{s}$ approaches unity.

Another fundamental distinction in the Brownian motion of unstable, compared to stable, particles concerns the first-passage time (FPT) to absorbing boundaries. The FPT of a stable particle is infinite irrespective of how close the initial position is to a single absorbing boundary. Mathematically, this is a consequence of the heavy tail―or asymptotic power law behavior―of the associated probability density (Smirnov density). In contrast, the FPT of an unstable particle is finite because the exponential decrease in time of the probability density is faster than that of a power law. Indeed, even if the unstable particle is free to walk the entire real axis, the FPT to reach $-\infty $ or $+\infty $ is finite and equal to the statistical lifetime ${\lambda}^{-1}$ of the particle.

Cite this paper

Silverman, M.P. (2017) Brownian Motion of Decaying Particles: Transition Probability, Computer Simulation, and First-Passage Times. Journal of Modern Physics, 8, 1809-1849. https://doi.org/10.4236/jmp.2017.811108

References

- 1. Einstein, A. (1905) Annalen der Physik, 17, 549-560. Reproduced in Stachel, J. (1989) The Swiss Years Writings: 1900-1909, 2, 223-236.
- 2. Bachelier, L. (1900) Annales scientifique de l’Ecole Normale Superieure, Series 3, 17, 21-86.
- 3. Gardiner, C. (2009) Stochastic Methods: A Handbook for the Natural and Physical Sciences. Springer-Verlag, Berlin.
- 4. Dattagupta, S. (2014) Diffusion: Formalism and Applications. CRC Press, New York, NY.
- 5. Mazo, R.M. (2002) Brownian Motion: Fluctuations, Dynamics, and Applications. Oxford University Press, Oxford, UK.
- 6. Crank, J. (1979) The Mathematics of Diffusion. 2nd Ed., Oxford University Press, Oxford, UK.
- 7. Ghez, R. (2001) Diffusion Phenomena. Kluwer, New York, NY. https://doi.org/10.1007/978-1-4757-3361-7
- 8. Gillespie, D.T. and Seitaridou, E. (2013) Simple Brownian Diffusion. Oxford University Press, Oxford, UK.
- 9. McCauley, J.L. (2013) Stochastic Calculus and Differential Equations for Physics and Finance. Cambridge University Press, Cambridge, UK. https://doi.org/10.1017/CBO9781139019460
- 10. Lamarsh, J.R. and Baratta, A.J. (2001) Introduction to Nuclear Engineering. Pearson, New York, NY.
- 11. Byrne, J. (2011) Neutrons, Nuclei, & Matter: An Exploration of the Physics of Slow Neutrons. Dover, New York, 276-296, 318-323.
- 12. Brenner, H. and Edwards, D. (1993) Macrotransport Processes. Butterworth-Heine-mann, Boston, MA, USA, 407-470.
- 13. Harrison, L.G. (1993) Kinetic Theory of Living Pattern. Cambridge University Press, Cambridge, UK, 173-190. https://doi.org/10.1017/CBO9780511529726
- 14. Murray, J.D. (2002) Mathematical Biology I: An Introduction. Springer, New York, 437-459.
- 15. Silverman, M.P. (2016) World Journal of Nuclear Science and Technology, 6, 232-260. https://doi.org/10.4236/wjnst.2016.64024
- 16. Uhlenbeck, G.E. and Ornstein, L.S. (1930) Physical Review, 36, 823-841. https://doi.org/10.1103/PhysRev.36.823
- 17. Merzbacher, E. (1970) Quantum Mechanics. 2nd Ed, Wiley, New York, 475-486.
- 18. Silverman, M.P. (2014) A Certain Uncertainty: Nature’s Random Ways. Cambridge University Press, Cambridge, UK, 122-128. https://doi.org/10.1017/CBO9781139507370
- 19. Silverman, M.P., Strange, W., Silverman, C.R. and Lipscombe, T.C. (1999) Physics Letters, A262, 265-273. https://doi.org/10.1016/S0375-9601(99)00668-4
- 20. Silverman, M.P. (2016) Search for Anomalies in the Decay of Radioactive Mn-54. Europhysics Letters, 114, 62001(1-6).
- 21. Wiersema, U.F. (2008) Brownian Motion Calculus. Wiley, Chichester, UK, 103-125.
- 22. Mantegna, R.N. and Stanley, H.E. (2000) An Introduction to Econophysics: Correlations and Complexity in Finance. Cambridge University Press, Cambridge, UK, 4, 17.
- 23. Redner, S. (2001) A Guide to First-Passage Processes. Cambridge University Press, Cambridge, UK. https://doi.org/10.1017/CBO9780511606014
- 24. Risken, H. (1996) The Fokker-Planck Equation: Methods of Solution and Applications. 2nd Ed., Springer-Verlag, Heidelberg, 29-31. https://doi.org/10.1007/978-3-642-61544-3_4
- 25. Wallace, P.R. (1984) Mathematical Analysis of Physical Problems. Dover, New York, 355-383.
- 26. Silverman, M.P. (2000) Probing the Atom: Interactions of Coupled States, Fast Beams, and Loose Electrons. Princeton University Press, Princeton, NJ, USA, 117-138.
- 27. Silverman, M.P., Strange, W. and Lipscombe, T.C. (2004) American Journal of Physics, 72, 1068-1081. https://doi.org/10.1119/1.1738426
- 28. Gibbs, J.W. (1902) Elementary Principles in Statistical Mechanics, Developed with Especial Reference to the Rational Foundation of Thermodynamics. Reproduced by Dover, New York, 1960.
- 29. Panagiotopoulos, A.Z. (1995) Gibbs Ensemble Techniques. In: Baus, M., et al., Eds., Observation, Prediction, and Simulation of Phase Transitions in Complex Fluids, NATO ASI Series C, Kluwer, Dordrecht, The Netherlands, Vol. 460, 463-501. https://doi.org/10.1007/978-94-011-0065-6_11
- 30. Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H. and Teller, E. (1953) Journal of Chemical Physics, 21, 1087-1092. https://doi.org/10.1063/1.1699114
- 31. Hastings, W.K. (1970) Biometrika, 57, 97-109. https://doi.org/10.1093/biomet/57.1.97
- 32. Silverman, M.P. (2017) World Journal of Nuclear Science and Technology, 7, 252-273. https://doi.org/10.4236/wjnst.2017.74020
- 33. Feller, W. (1950) An Introduction to Probability Theory and Its Applications. Vol 1. Wiley, New York, 304-305, 411-413.
- 34. Mood, A.M., Graybill, F.A. and Boes, D.C. (1974) Introduction to the Theory of Statistics. 3rd Ed., McGraw-Hill, New York, 72-81.
- 35. Carslaw, H. and Jaeger, J. (1959) Conduction of Heat in Solids. 2nd Ed., Oxford University Press, Oxford, UK, 273-281.
- 36. Berg, H.C. (1993) Random Walks in Biology. Princeton University Press, Princeton, 42-47.
- 37. Cap, F. (1970) Einführung in die Plasmaphysik I: Theoretische Grundlagen. [Introduction to Plasma Physics I: Theoretical Foundations.] Akademie-Verlag, Berlin, 16-18.
- 38. Bellan, P.M. (2006) Fundamentals of Plasma Physics. Cambridge University Press, Cambridge, UK, 569-579. https://doi.org/10.1017/CBO9780511807183
- 39. Jackson, J.D. (1999) Classical Electrodynamics. 3rd Ed., Wiley, New York, 34-35. https://doi.org/10.1119/1.19136
- 40. Hughes, B.D. (1995) Random Walks and Random Environments: Random Walks. Oxford University Press, Oxford, UK, Vol. 1, 205-206.
- 41. Folks, J.L. and Chihikara, R.S. (1978) Journal of the Royal Statistical Society: Series B, 40, 263-289.
- 42. Callen, H.B. and Welton, T.A. (1951) Physical Review, 83, 34-40. https://doi.org/10.1103/PhysRev.83.34
- 43. Coffey, W.T. and Kalmykov, Y.P. (2012) The Langevin Equation. 3rd Ed., World Scientific, Singapore, 16-17. https://doi.org/10.1142/8195
- 44. Serber, R. (1992) The Los Alamos Primer. University of California Press, Los Angeles, CA, USA, 65-75.
- 45. Sasaki, T., Gunji, Y. and Iida, T. (2007) Journal of Nuclear Science and Technology, 44, 1330-1336. https://doi.org/10.1080/18811248.2007.9711379

Appendix 1: Random Variable ${\Sigma}_{n}^{2}$

Random variable ${\Sigma}_{n}^{2}$ is defined in Equation (62) as a sum of products of Bernoulli random variables ${\epsilon}_{j}\left(j=1,\cdots ,n\right)$ for a discrete Brownian motion of n time steps. The moment generating function (MGF) is defined by the expectation value

$\begin{array}{c}{g}_{{\Sigma}_{n}^{2}}\left(\theta \right)\equiv \langle \mathrm{exp}\left({\Sigma}_{n}^{2}\theta \right)\rangle ={\displaystyle \underset{\left\{\epsilon \right\}}{\sum}\left({\displaystyle \underset{j=1}{\overset{n}{\prod}}{\text{e}}^{\theta \left({\epsilon}_{j}\cdots {\epsilon}_{n}\right)}}{\displaystyle \underset{j=1}{\overset{n}{\prod}}{\pi}_{j}}\right)}\\ ={\displaystyle \underset{\left\{\epsilon \right\}}{\sum}\left({\text{e}}^{\theta \left({\epsilon}_{1}\cdots {\epsilon}_{n}\right)}{\text{e}}^{\theta \left({\epsilon}_{2}\cdots {\epsilon}_{n}\right)}{\text{e}}^{\theta \left({\epsilon}_{3}\cdots {\epsilon}_{n}\right)}\cdots {\text{e}}^{\theta {\epsilon}_{n}}\right){\pi}_{1}{\pi}_{2}\cdots {\pi}_{n}}\end{array}$ (98)

in which ${\pi}_{j}$ is the probability of the ${j}^{\text{th}}$ Bernoulli outcome

${\pi}_{j}=\{\begin{array}{l}p\text{if}{\epsilon}_{j}=1\\ q=1-p\text{if}{\epsilon}_{j}=0\end{array}$ (99)

and the sum represented by $\left\{\epsilon \right\}$ is over all possible outcomes of the Bernoulli variables.

An illustration of the simple case $n=3$ can give a sense of the content of Equation (98). All possible outcomes with corresponding expectations are given in Table 3.

Summing the results of Table 3 and substituting $q=1-p$ lead to the expression for ${g}_{{\Sigma}_{3}^{2}}\left(\theta \right)$

${g}_{{\Sigma}_{3}^{2}}\left(\theta \right)={p}^{3}{\text{e}}^{3\theta}+{p}^{2}\left(1-p\right){\text{e}}^{2\theta}+p\left(1-p\right){\text{e}}^{\theta}+\left(1-p\right)$ . (100)

Proceeding in the same way for $n=4$ leads to ${g}_{{\Sigma}_{4}^{2}}\left(\theta \right)$

${g}_{{\Sigma}_{4}^{2}}\left(\theta \right)={p}^{4}{\text{e}}^{4\theta}+{p}^{3}\left(1-p\right){\text{e}}^{3\theta}+{p}^{2}\left(1-p\right){\text{e}}^{2\theta}+p\left(1-p\right){\text{e}}^{\theta}+\left(1-p\right)$ . (101)

The above pattern holds for arbitrary n, whereupon Equation (98) reduces to the form

${g}_{{\Sigma}_{n}^{2}}\left(\theta \right)={p}^{n}{\text{e}}^{n\theta}+\left(1-p\right){\displaystyle \underset{k=1}{\overset{n-1}{\sum}}\left({p}^{k}{\text{e}}^{k\theta}\right)}+\left(1-p\right)$ , (102)

which, upon summing the geometric series in Equation (102) yields the compact

Table 3. Outcomes of 3 Bernoulli trials.

expression

${g}_{{\Sigma}_{n}^{2}}\left(\theta \right)=\frac{\left(1-p\right)+{p}^{n+1}{\text{e}}^{n\theta}\left(1-{\text{e}}^{\theta}\right)}{1-p{\text{e}}^{\theta}}$ . (103)

In the analyses of Sections 2.3 and 2.4 the survival probability p was eventually replaced by the expression (27), leading to the limit ${p}_{s}^{n}={\left(1-\lambda t/n\right)}^{n}\underset{n\to \infty}{\to}{\text{e}}^{-\lambda t}$ . However, if in Equation (103) p is taken to be a fixed finite parameter satisfying the requirement of a probability $1\ge p\ge 0$ , then ${p}^{n}\to 0$ for $n\to \infty $ . In that limit, MGF (103) reduces to the approximate form

$\underset{n\to \infty}{\mathrm{lim}}{g}_{{\Sigma}_{n}^{2}}\left(\theta \right)=\frac{\left(1-p\right)}{1-p{\text{e}}^{\theta}}$ (104)

which can be approximated further

$\underset{\begin{array}{l}n\to \infty \\ \theta \ll 1\end{array}}{\mathrm{lim}}{g}_{{\Sigma}_{n}^{2}}\left(\theta \right)=\frac{\left(1-p\right)}{1-p\left(1+\theta \right)}=\frac{1}{1-\left(p/q\right)\theta}$ (105)

by expansion to first order of the exponential about $\theta =0$ . (Note: $\theta $ is always set to 0 after differentiation of the MGF.) Expression (105) is the MGF of the exponential distribution $E\left(p/q\right)$ [Ref. [34] , pp. 540-541]; the form of the PDF is given by Equation (54). A characteristic feature of the exponential distribution is that the standard deviation $\sigma $ equals the mean $\mu $ , which is equal to the parameter of the distribution. In the present case the parameter is $p/q=p/\left(1-p\right)$ .

A test of this relation is shown in Figure 12, which plots $\mu \left(p\right)$ (solid red)

Figure 12. Plots of the mean $\mu \left(p\right)=\langle {\Sigma}_{n}^{2}\rangle $ (solid red) and standard deviation $\sigma \left(p\right)=\sqrt{\langle {\left({\Sigma}_{n}^{2}\right)}^{2}\rangle -{\langle {\Sigma}_{n}^{2}\rangle}^{2}}$ (dashed blue) as a function of survival probability p for $n=1000$ , calculated from the exact expressions (70) for the first and second moments. In the limit of $n\gg 1$ , the ratio $\mu \left(p\right)/\sigma \left(p\right)$ (dotted black) is equal to $\sqrt{p}$ (dashed green).

and $\sigma \left(p\right)$ (dashed blue), given by the exact expressions (70) derivable from the MGF (103), over the full range of p for $n=1000$ . The equality is very close at the scale of the figure, but the higher-resolution insert shows that $\sigma \left(p\right)$ exceeds $\mu \left(p\right)$ by a small amount that decreases as p approaches 1. The theoretical ratio $\mu \left(p\right)/\sigma \left(p\right)$ (dotted black), equals $\sqrt{p}$ (dashed green) in the limit $n\gg 1$ . As a sum of correlated products of Bernoulli random variables, the compound random variable ${\Sigma}_{n}^{2}$ does not lend itself to an obvious physical interpretation. Nevertheless, its approximate exponential distribution in the asymptotic limit $p\to 1$ (or $q=\left(1-p\right)\to 0$ ) might be understood in the following way. A negative exponential distribution characterizes the time intervals between events in a Poisson process. The Poisson distribution itself, however, evolves from a binomial distribution under the circumstances that the number of samples $n\gg 1$ and the probability of an event (e.g. the event of particle decay in the present context) $q\ll 1$ such that $nq$ approaches a constant mean number of events. Given that a Bernoulli random variable is a special case of binomial random variable, the asymptotic limit at which ${\Sigma}_{n}^{2}$ is well-represented by a negative exponential distribution is precisely the condition for a Poisson process.

Appendix 2: Validity of Fick’s Law for a Decaying Particle

The point at issue is whether diffusion theory based on Fick’s law (14) is valid for a particle with a finite statistical lifetime. Serber [Ref. [44] ], one of the theoretical physicists of the Manhattan Project, raised this question in the context of the diffusion of neutrons. Although a free neutron has a half-life of about 10 minutes, the problem faced by Serber did not involve neutron decay, but the loss of neutrons through the surface of the host material. According to Serber, ordinary diffusion theory is valid only when the size of the confining region is large compared to the mean free path of the diffusing particles. Nevertheless, since the process of neutron diffusion terminates the moment a neutron escapes the material, the theoretical problem faced by Serber is identical to the problem of termination of Brownian motion by radioactive decay.

In terms of notation used in this paper, the criterion derived by Serber for the validity of Fick’s law is the following

$\frac{1}{3}\mathcal{l}\left|\frac{\text{d}j\left(x\right)}{\text{d}x}\right|\ll \frac{1}{2}j\left(x\right)$ (106)

^{2}Jeans, J. (1954) The Dynamical Theory of Gases. 4^{th} Ed., Dover Publications, New York, 307-310.

^{3}Reif, F. (1965) Fundamentals of Statistical and Thermal Physics. McGraw-Hill, Boston, 483-486.

in which
$j\left(x\right)$ is the particle flux to the right or to the left along the x-axis, and
$\mathcal{l}$ is the mean free path between collisions with particles of the medium. The diffusion coefficient D, derivable from the elementary kinetic theory of gases to an approximation adequate for the purpose of this appendix, is^{2}^{,}^{3}

$D=\frac{1}{3}\mathcal{l}{v}_{\text{rms}}$ (107)

in which the root-mean-square speed of particles of mass m in equilibrium at temperature ${T}_{\text{e}}$ is

${v}_{\text{rms}}=\sqrt{3{k}_{\text{B}}{T}_{\text{e}}/m}$ (108)

where ${k}_{\text{B}}$ is Boltzmann’s constant. The self-diffusion coefficient (107) follows in the simple case of a dilute gas by calculation of the net flux $j\left(x\right)={v}_{\text{rms}}\left[n\left(x-\mathcal{l}\right)-n\left(x+\mathcal{l}\right)\right]/6$ of particles originating from a distance $\mathcal{l}$ right and left of a plane perpendicular to the flow. $n\left(x\right)$ is the number density of particles; the factor 1/6 arises because on average 1/3 of the particles move along the x-axis, and, of these, 1/2 move in the positive (negative) direction. Expansion in a Taylor series to first order in $\mathcal{l}$ then leads to

${j}_{x}\approx -\frac{1}{3}{v}_{\text{rms}}\mathcal{l}\left(\text{d}n/\text{d}x\right)\equiv -D\left(\text{d}n/\text{d}x\right)$ (109)

which defines the diffusion coefficient D.

From the macroscopic theory of diffusion of radioactive particles [Ref. [15] ], one can show that

$\left|\frac{1}{j\left(x\right)}\frac{\text{d}j\left(x\right)}{\text{d}x}\right|=\sqrt{\frac{D}{\lambda}}\equiv \zeta $ (110)

where $\zeta $ is the characteristic diffusion length. Combining relations (106) through (110) leads to the succinct expression for the validity of Fick’s law applied to radioactive particles

$2{v}_{\text{d}}\ll {v}_{\text{rms}}$ (111)

in which

${v}_{\text{d}}\equiv \sqrt{D\lambda}$ (112)

is the characteristic diffusion velocity. Substitution of relations (112) and (108) into (111) yields ${\lambda}^{-1}\gg 4mD/3{k}_{\text{B}}T$ , which is the inequality (97).