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

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.


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 timedependent perturbation theory in quantum mechanics [17], that the probability of particle decay within a short time interval t ∆ takes the form ( ) where λ is the intrinsic decay rate, it is readily deducible [18] that the proba- bility density of survival to time t satisfies the differential equation ( ) ( ) with normalized solution ( ) e t p t λ λ − = . (3) From Equation (3) it follows that the mean lifetime τ is the reciprocal of λ ( ) and that the half-life, i.e. duration within which the survival probability is 50%, is The premise that the probability (1), with λ 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 t 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] ( ) ( ) ( ) ( ) ( ) ( )  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 ( ) one can express the corresponding random variable as and the sum of two independent Gaussian random variables as 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) ( ) 0 0 , , p x t x t of finding the particle at location x at time t, given that the particle was initially at 0 x at time 0 t t < .The partial differential equation corresponding to the stochastic differen- tial Equation ( 6) is the forward Fokker-Planck equation [Ref [3], pp.117-122]

A x t p x t x t B x t p x t x t p x t x t t
x x 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) ( ) , 2 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.

Fokker-Planck Equation and Transition Probability
Consider a quantity Application of the divergence theorem then leads from the integral relation (12) to the differential equation for the conservation law of a disintegrating quantity.
Upon dividing Equation ( 13) by the initial number of particles 0 1 N  and relating the current density to the gradient of the particle density by Fick's law with diffusion coefficient D (here taken to be a constant), one obtains an equation of the form x is interpretable as the probability density for diffusion of a single decaying particle.However, it is demonstrable that the solution ( ) ,  w t x under the special initial condition at 0 t ( ) ( ) is identical to the transition probability density ( ) , which is the conditional probability for an unstable particle to have reached position x at time t given that it was at 0 x at 0 t t < .The equivalence is readily established by substitution of condition ( 16) into the Chapman-Kolmogorov equation [24] ( ) w t p t t w t x x x x x x x (17) In the case of one-dimensional Brownian motion treated in this paper, Equation (15) then takes the form which, together with the initial delta-function condition ( 16) ) 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 ( From the form of relation (20) 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 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 th k power of displacement of a particle in time interval t starting from po- 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 meansquare displacement.Employing solution (20) as the Gaussian PDF ( ) , w x t in the expectation value (21), one obtains the mean location 1 0 m x = and variance about the mean of a non-decaying diffusing particle.However, calculation of expectation values (23) with the TPF ( ) 0 , ,0 p x t x leads to the initial location 0 x x = of the diffusing particle and its mean-square distance from the point of origin ( ) 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.

Mean-Square Displacement
Consider a one-dimensional Gaussian random walk on a lattice with time step t δ and displacement j X δ during the jth time step given by ( ) where the lattice spacing 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 t δ is , and the probability that the process ends at a particular time step is t λδ .Thus, although the distance traveled in the jth time step is de- termined by a unit normal distribution ( ) , the probability that the displacement is made at all is determined by a Bernoulli random variable ( ) A Bernoulli distribution is the special case 1 n = of the binomial distribution ( ) , B n p of n trials with probability of success p ; the corresponding discrete probability function is 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 n X ( ) ( ) The expectation value of the mean square displacement at time t n t δ which can be re-expressed in the form Upon defining the diffusion constant D in the standard way and taking the limit 0 t δ → , 0 X δ → , with requirement that D remain finite, 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 x t in the Wiener term of the forward Fokker- Planck Equation (10).Under the conditions that ( ) is constant and there is no drift, ( ) , 0 A x t = , Equation (10) then reduces to the standard one- dimensional diffusion equation in which D is the physically measurable diffusion constant first introduced by Einstein in his theory of Brownian motion.

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 ( ) ( ) 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 t n is a random sample from the unit normal distribution ( ) , where the subscript t and superscript d t t + explicitly denote the temporal range with respect to which t n is asso- ciated.Thus, two samples n , corresponding to distributions ( ) ( ) , are independent for 2 1 d t t t − > .Given an initial value ( ) , a sequence of points by iterative use of Equation ( 35) up to time d t n 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 d 0 t → .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 x t and ( ) , B x t 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 does not involve the spatial derivative of ( ) 0 0 , , p x t x t .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 ( ) Integration of Equation ( 36) yields with boundary and initial conditions The error function is defined by the integral 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 ( ) , which determines whether or not the particle survives the interval dt with a survival probability s p given by Equation (27).
The appropriate Langevin equation would then take the update form where the Bernoulli variate t ε is defined in relation (27).Note that an occur- rence of 0 t ε = terminates the entire process by setting ( ) just before the moment of decay, the particle will have reached location ( ) 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 0 0 x = are shown explicitly below: The mean displacement ( ) d 0 x n t = follows immediately from the properties of the unit Gaussian ( ) n N + = = .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 dt t n = and eventually takes the limit n → ∞ and d 0 t → such that t is a fixed quantity.Equation ( 41) then yields where the Bernoulli variates 1 t ε = 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 ( ) summing the terms in Equation ( 42), and taking the forementioned limit, one arrives at 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 ( ) , s B n p .By comparison, one can think of the outcome (44) as resulting from a sum of n compound stochastic processes of the form ( ) ( ) It is of interest to examine the two limiting cases of Equation ( 44): 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 rootmean-square distance ( ) x t reached by a particle very likely to decay during the diffusion time is equivalent, within a factor 2 , to the characteristic diffusion length [15] 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 x t appearing in the Fokker-Planck Equation (10).In this paper, therefore, the mathematical diffusion length will be defined Consider next the mean-square displacement as obtained from Equation ( 44) with n discrete time steps of duration t δ ( ) ( ) ( ) where use was made of the summation formula for a geometric series 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 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 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 ( ) expresses the probability of a process with n successes in sequence followed by a single failure that terminates the process.The distribution is normalized The outcome of a large number of Brownian random walks of a decaying particle, each starting from the same initial condition 0 0 x = , leads to a distribu- tion of time steps n given by Equation (50).One can ask, therefore, for the ensemble-average of the mean-square displacement (49) where, again, the term 1 t λδ  was dropped.The result (53) of the ensemble average is a root-mean square displacement 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 ( ) The distribution is normalized The ensemble average of the mean-square displacement (49) is then 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].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)).
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 1000 n = for the same values of D and s p 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

Langevin Equation: Analytical Solution
It is possible to solve stochastic Equation (40) analytically in closed form for the displacement ( ) X t and PDF ( ) , X p x t 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 1 ε = , then ( ) X t must also be expressible in closed form as a normal random varia- ble 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 d n t t → ) must be ( ) whereupon the corresponding PDF, generalized on the basis of space and time translational symmetry to include initial conditions ( ) 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 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 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 ε , in departure from the conventional notation to use an upper-case letter.
The question then becomes: What kind of random variable is (which always exists), or its probability function (for discrete outcomes) or probability density (for continuous outcomes) [34].The variable θ in the ar- gument 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 2 n Σ in Appendix 1, the MGF was used progressively, start- ing from the relation ( ) and If the set of random variables { } j ω 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 ω variables, e.g.
 , could contain some identical factors and be highly correlated.Nevertheless, although less easily done, the MGF of 2 n Σ can be shown to be thereby characterizing 2 n Σ uniquely.
For example, moments 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 (Note: ( ) ) Thus, the random walk of the decayed particle has been realistically terminated at the randomly selected th k time step.With regard to notation, the subscript on 2 n Σ 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 2D λ , 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  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 ( ) , and variance 2D λ , whereas the dis- tribution 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 e t λ − , and the meansquare displacement of its Brownian motion spreads as a power law 2Dt .Journal of Modern Physics 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 λ .18) for the same time units, process parameters, and color code as in Figure 3.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 ( ) decreases monotonically in both cases, but asymptotically approaches 0 in the Fokker-Planck PDF and ( ) 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.

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.0 0 x = ), first reach some other speci- fied 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  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 expec- T p t is the probability density function (PDF) of T and ( ) 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 FPT T for a stable particle.One approach involves calculation of ( ) with the solution to the Fokker-Planck equation.The other approach yields FPT T directly as the solution to a Poisson equation.These two methods are ge-Journal of Modern Physics neralized 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 b x = ± before decaying.With the methods developed below the boundary conditions can be easily modified to apply to other systems.

Method 1: Calculation of TFPT by Image Functions
The problem is to find a PDF ( ) x p x t of the diffusing particle with initial dition (19) that satisfies boundary conditions . 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 in which the second line of (77) defines the function  .At this time, there is a significant probability that the unconstrained particle (black), described by ( ) ( ) x t , could be found outside the designated physical region.
However, the 0 x t and a negative Gaussian image centered on 2 , satisfies the right- side boundary condition ( ) , although it does not satisfy the boundary condition on the left side.Figure 7 shows that the 1 N = component with su- perposition of ( ) , , 1 which is the probability that the particle has not reached either boundary at time t.The probability ( ) S t that the particle remains in the interval [ ] where ζ , 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 The resulting expression, however, is complicated and not needed in this paper.

Method 2: Calculation of TFPT from a Screened Poisson Equation
The method employed in Section 3.1 yielded the CPF ( ) T F t (or, equivalently, the survival function ( ) S t ) from which all statistical properties of the FPT can be calculated.However, if all one wanted was FPT T , 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 t δ , on a lattice with coordinate spacing x δ and absorbing boundaries at b x = ± .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 x is the time to reach ei- ther boundary from point x, it must satisfy the following discrete equation ) which is expressible as a second derivative leads in the limit 0 x δ → and 0 t δ → to a screened Poisson equation, en- countered in the physics of ionized gases [37] [38], ( ) ( ) in which the ratio ( ) δ δ is again taken to be the diffusion coefficient D as defined in Equation (32).In the absence of the term proportional to λ , Equa- tion (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 and in all expressions involving ( ) T x refers to the initial point from which the particle diffuses to a boundary, and therefore actually corresponds to coordinate 0 x in Equation (86).Where there is no confusion, it is standard notation to use the variable x rather than 0 x as the argument of the FPT Equation (85) and solution.
The solution to Equation (85) with implementation of boundary conditions Comparison with solution (80) of Section 3.1 can be made by setting the initial location 0 x = , whereupon Equation (87) reduces to Methods 1 and 2, although very different, lead to the same final result, as they must for consistency.

Critical Role of Particle Decay
The Fokker-Planck Equation ( 18 In the limit a → −∞ , there is only a single boundary b on the right, and the FPT (89) reduces in this limit to ( ) in which the λ dependence is explicitly shown.Note that ( ) b T x (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 → ∞ as well, the particle is free to walk the entire x-axis, and the FPT reduces to the statistical lifetime (4) as expected.
For a stable particle ( ) 0 λ = , 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 The first term of (91), in which the denominator is the characteristic diffusion velocity [Ref [15]], is singular as 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, The first term in brackets, known as the Smirnov density [40] or an inverse Gaussian density [41], is independent of λ and characterized by a heavy tail, i.e. an asymptotic power law 3 2 t − .It is this density that leads to an infinite first moment FPT T for a stable particle

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 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?
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 p m τ µ ≡ , in which m is the particle mass and p µ is the particle mobility 1 , is large enough that the displacement ( ) x t τ + is independent of the displacement ( ) x t , but small compared to the time interval t ∆ between observations [43].This criterion leads to an inequality [Ref.[8], pp.66-67]  in which D is the particle diffusion coefficient, B k is Boltzmann's constant, and e T is the equilibrium temperature.For a radon-222 atom diffusing at room temperature e 300 K T ≈ , Equation (96) becomes , 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

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 λ or statistical lifetime 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 ( ) 0 0 , , p x t x t , and the Langevin Equation (LE), which yields the distribution function and trajectories of the corresponding process variable ( ) X t .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, 2 n Σ , although not among any known to the author, depends on the number of time steps n and survival probability s p 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, 2 ∞ Σ behaves increasingly like an exponential distribution as s p 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 −∞ or +∞ is finite and equal to the statistical lifetime 1 λ − of the particle.
is a deterministic function, whereas the Wiener term is the source of fluctuations.

Figure 1 .
Figure 1.Computer simulation of Brownian motion with decay.Parameters: 500 n = time steps d 1 t = with diffusivity 0.5 D = and survival probability 0.99 s p =

Figure 2 .
Figure 2. Computer simulation of Brownian motion with decay, employing 1000 n = 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 inTable1 and Table 2, are in accord with predictions of Equations ((53) and (80)).

yield 10 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 5000 n = 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 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 nΣ
68) does not correspond to any of the tabulated random variables known to the author, all the moments of can be determined by differentiation ( ) 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)

2 nΣ
72)upon substitution of relation(27) for s p 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 and the 4 th mo- ment of the displacement

Figure 4 .
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 Figure3.

and Figure 2 , 33 dt
two absorbing boundaries were arbitrarily set at 15 b x = ± units from the origin.As shown in Figure1, the unstable particle (black) first reached boundary 15 however, that a prior decay occurred close to 1 = .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 In the simulation shown in Figure2, there are four decays, the last occurring at time step 4 the structure of solution (77) is provided by Figures 6-10.The fig-ures depict the probability density of an unstable particle with diffusivity and right of the origin.In each figure the dashed black plot is the Gaussian PDF ( ) ( ) particle at stated time t; the solid red plot is the specified linear superposition of functions n 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.

X p x t + , which include Gaussian images centered on 4 −  , 2 −
Figure 9. Figure 8 shows that the 0 N = solution satisfies just the right boun-

Figure 6 .
Figure 9. Figure 8 shows that the 0 N = solution satisfies just the right boun- dary condition; Figure 9 shows that the 1 N = solution still does not satisfy the

Figure 8 .
Figure 8. Spatial variation of the 0 N = component to

Figure 9 .Figure 10 .
Figure 9. Spatial variation of the 0, 1 N = ± components to 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 ( .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 0 x = .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 :
Figure 2: 0.5 D = , 82)because: a) point x can be reached only from points x x δ + and x x δ − , b) the probability of a transition from these two points is the same (1/2), and c) the transition to x during interval t δ can be made only if the particle has survived the transition.Recall from Equation (27) that the survival probably per step is 1 of Equation (82) to the form

T
) incorporates particle instability by means of a decay term resulting transition probability density(20) only through a global exponential factor e t λ − .Thus, letting λ ap- proach 0 smoothly generates the probability density of the unconstrained stable particle.Likewise, the vanishing of λ 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 ( ) ( ) 0 gained by examining the probability density of the FPT in the case of a single absorbing boundary at 0 b > and a particle initially located at 0.( )T p t is obtained by method 1 of the preceding section as applied to PDF (77) truncated to include only the component ( ) centered on 0 and a negative Gaussian image centered on 2b suffice to main- tain = for all time t.The calculation yields Smirnov density multiplied by an exponential e t 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

T
which, together with expression (95), sum to() obtained directly from the screened Poisson Equation (85).

Figure 11 .
Figure 11.Temporal variation of the probability density

1
Mobility pµ is the proportionality factor in p V F µ = relating velocity V and dissipative force F.
radon-222 at room temperature as an example, one has very well satisfied.In fact, relation (97) is well satisfied even for the short-lived progeny of radon-222, such as polonium-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].

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

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