1. Introduction: Brownian Motion with Decay
1.1. Fokker-Planck and Langevin Equations
The diffusion of radioactive particles, particularly in the gaseous and aerosol state, plays an important role in the dispersal, detection, and monitoring of a variety of radioactive isotopes affecting public health and safety. Examples include nuclides that arise naturally in the environment such as radon 222Rn and thoron 220Rn, or are associated with power reactor releases such as iodine 129I and 131I, cesium 131Cs and strontium 90Sr, or are produced in nuclear weapons manufacture and testing such as tritium 3H, among many other possibilities  .
Although diffusion as a physical process has been investigated theoretically since Einstein’s seminal work on Brownian motion in 1905  , the preponderance of studies has focused on stable matter that did not undergo transformations. A few exceptions known to the authors pertain to the diffusion or percolation of radon as modeled by steady-state solutions to a macroscopic transport equation based on Fick’s law. See, for example,    for examples of steady-state models of radon diffusion in air, water, and percolation through construction materials. The primary objective of these kinds of models was to obtain the equilibrium concentrations of radon and its progeny in indoor environments or mines.
In a recent publication  one of the authors (MPS) developed the theory of diffusion of unstable particles from a comprehensive micro-statistical perspective by deriving and solving the associated time-dependent Fokker-Planck and Langevin equations for Brownian motion with decay. The three-dimensional (3D) Fokker-Planck equation (FPE) for this process, which takes the form
leads to the transition probability density for finding a particle at location at time t, given that it was at at time . The constant parameters of the equation are the diffusion coefficient (or diffusivity) D and the intrinsic nuclear decay rate λ. The solution to Equation (1) in the case of one-dimensional (1D) Brownian motion of a particle located at at was shown to be
and is readily generalized for independent Brownian motion in three dimensions (as will be discussed in Section 2.3).
In contrast to the FPE, the Langevin equation (LE) expresses the time-variation of a process random variable such as displacement or velocity, rather than the transition probability density. In the case of 1D Brownian motion of a decaying particle, the LE derived in  takes the form of an update stochastic differential equation   based on a Wiener process  
whose graphical solution shows the sequential displacements of the particle as a function of time and whose analytical solution yields the statistical distribution of the displacement variable X.1 The symbol , which governs particle displacement, represents a random sample from the unit normal distribution2 where the subscript t and superscript explicitly denote the temporal range with respect to which is associated. In other words, two samples and corresponding to distributions and are independent for . The symbol , which governs particle decay, represents a Bernoulli random variable3 defined by
in which the probability of radioactive decay within the differentially small time interval dt is given by , the transformation law first recognized by Rutherford and Soddy and used by von Schweidler to derive the law of exponential decay  . Thus, defined in Equation (4), is the probability of surviving the interval dt. The analytical solution to Equations ((3) and (4)) for a decaying particle undergoing k discrete one-dimensional Brownian displacements in time t at intervals of dt was shown in  to be
Note that the expression in (5) is not a function, but a statistical distribution (Gaussian) whose variance after k particle displacements is expressed in terms of another random variable that was derived and explained in  (and will be used in Section 2.2).
For a system of non-decaying particles, the statistical content of the FPE is ordinarily equivalent to that of the LE. It is a matter of choice which method to apply, and there is a standard formulaic procedure for generating the LE from the FPE and vice versa   . It was shown theoretically in  , however, that for decaying particles FPE (1) and LE (3) can yield different, but complementary, information. For example, a statistical measure of the dispersal of a radioactive gas by diffusion is given primarily by the first two moments, i.e. the mean displacement (which is 0 in the absence of advection) and mean-square displacement . FPE solution (2) and LE solution (5) (for 1D Brownian motion in the continuum limit) led to expressions
where the mean-square displacement about the mean is defined by .
Now two theoretical approaches to the same physical problem of Brownian motion of radioactive particles may yield complementary information, but they cannot yield inconsistent information and both be correct. In this paper we establish the complete equivalence of the FPE and LE approaches to the Brownian motion of radioactive particles by
1) deriving probability density functions from the FPE that account for the distribution of particle displacements obtained from the LE;
2) clarifying operationally the difference in predicted variances (6) and (7);
2) demonstrating by means of the random variable that the two approaches lead to the same statistical moments at all orders; and
4) testing the distributions derived from the FPE by a Monte Carlo method based on the LE.
1.2. Survival Function and First-Passage Time
Solutions to the FPE and LE address the general question of how far a diffusing particle has gone in a given time. The converse question of how much time a particle takes to diffuse a given distance is part of stochastic theory that analyzes first-passage processes  . Such processes are relevant to the quantitative measurement of radioactivity by methods that sample a radioactive source by means of diffusion    , to studies of the dispersal of radioactive contaminants in the atmosphere  , and to determinations of radiation exposure to critical organs such as the lungs  .
In  the survival function and corresponding mean first-passage time―i.e. the mean time for a particle to first reach one or more designated boundaries―was calculated for 1D Brownian motion of unstable particles to reach either of two absorbing boundaries. An absorbing boundary is one that removes a particle from the system. In mathematical terms, the probability density of the particle at the boundary is zero (in contrast to a reflecting boundary where the current density of the particle is zero). The survival function derived in  in the case of two absorbing boundaries at ,
is the probability that the particle, initially located at , has not reached either boundary at time t. The summation over an infinite number of terms in (8) arises from the requirement that the probability density of finding the particle at is 0 for all values of time t. Equation (8) was solved by a method of images  , and the number of image functions necessary to maintain the boundary conditions increases as the diffusion time increases. An infinite number is therefore required for the range .
In the analysis of first-passage processes of non-decaying particles, the probability density that a particle has not reached a boundary is obtained by differentiation of the survival function 
Correspondingly, the random variable T, for which expression (9) is the probability density (as indicated explicitly by the subscript on p), is referred to as the first-passage time (FPT). However, it is important to note that in the case of decaying particles, T takes on a broader meaning because the process of decay, as well as the process of diffusion, may be the cause of a particle failing to reach a boundary. Nevertheless, unless otherwise noted, for convenience of reference we will maintain the convention of referring to T as the first-passage time and to the density (9) as the FPT probability density.
In Ref.  the mean FPT for various boundary configurations was calculated by two methods that avoided the need to use the FPT density (9), and therefore the explicit form of was not given. The density will be derived in the present paper because it is useful, not only for practical reasons (e.g. to calculate FPT moments higher than the first), but also for conceptual reasons bearing on the consistency of the FPE and LE approaches to Brownian motion of radioactive particles.
In the survival function (8), the decay parameter λ is present only in the exponential factor. One might therefore surmise erroneously that the effect of particle decay would simply be to shorten the mean FPT relative to what it would be for a stable particle. The effects, however, are more profound than that. In the present paper, Equation (9) is evaluated explicitly and shown to take a form substantially different from the corresponding FPT density for non-decaying particles. To clarify these matters we
1) calculate the FPT probability density from the FPE and interpret its mathematical structure, and
2) test the expressions derived from the FPE (relations (8) and (9)) by a Monte Carlo simulation based on the LE applied to Brownian motion of radioactive particles confined between absorbing boundaries.
This paper is organized as follows. In Section 2 spatial distributions (1D and 3D) obtained from the solution to the FPE are derived, the mean-square displacements of the particles are calculated, an explanation is given for the two statistical expressions (6) and (7), and the equivalence of the FPE and LE analyses of Brownian motion of radioactive particles is demonstrated by showing that both equations give rise to the same statistical moments at all orders. In Section 3 the distribution of first-passage times is obtained from the solution to the FPE, and the probability density is shown to comprise two components, one relating to removal of a particle from the system by absorption at the boundaries, and the other relating to removal of a particle by radioactive decay. In Section 4 a Monte Carlo simulation based on the stochastic LE is used to generate distributions of Brownian trajectories under specified conditions to test the density functions derived from the FPE in Sections 2 and 3. Conclusions regarding the equivalence of the two approaches are summarized in Section 5.
2. Spatial Distributions
2.1. 1D Brownian Motion
As discussed in  , an examination of the forms of the two dynamical equations, (1) and (3), provides a qualitative insight into the seemingly discordant variances given by relations (6) and (7). FPE (1) is a continuous partial differential equation whose solution characterizes the continuous time-dependent probability of displacement of a specific particle throughout the full range of time . Since the probability of a radioactive particle surviving to time t is , the probability that the particle has decayed becomes 100% only at . In contrast, LE (3) is a discontinuous stochastic differential equation that describes a discontinuous set of Brownian trajectories with interruptions at each instant of actual particle decay, at which time the transformed particle is removed from the system and replaced by a new particle at the origin. Thus, the LE describes a type of Gibbsian ensemble  of Brownian trajectories, rather than a single continuous Brownian trajectory. In brief, the two equations describe Brownian motion from two different perspectives. In the case of a stable particle, there is no decay and therefore no removal and replacement; both the FPE and LE then describe the same continuously evolving Brownian motion of an arbitrary particle.
There is, however, a mathematical connection between the different perspectives of Brownian motion afforded by the FPE and LE, which relates the mean-square displacements (6) and (7) consistently. Consider first the interpretation of solution (2) to the 1D FPE, which for emphasis is re-expressed in the following way (with initial conditions taken to be , for simplicity)
The factor in square brackets is the probability density for the particle to reach location x at time t. The subsequent exponential decay factor is the independent probability that the particle has survived to time t. The product of the two factors in relation (10) is therefore the probability density for the Brownian particle to diffuse to x and survive―hence the subscript “S”.
Multiplying FPE solution (2) by leads to a differential expression with the same initial conditions as in (10)
in which the bracketed factor is the survival probability density (10), and the differential factor that follows is the probability of subsequent radioactive decay within interval dt. The subscript D on the probability density in Equation (11) signifies that the particle has decayed at x at some time between t and . Thus, the probability density that a radioactive particle decayed at x irrespective of the moment of decay within the finite interval is obtained by integration
is the characteristic diffusion length and
is the characteristic diffusion velocity (See Appendix 1 of Ref.  .).
Probability densities (10) and (13) are not individually normalized to unity because each density covers an incomplete set of possible outcomes. For example, the integral of density (10) over all available space
yields the probability of survival to time t. The sum of the two densities, however,
is normalized to unity
because it covers a complete set of mutually exclusive, uncorrelated Brownian processes that originate at : 1) survival at , 2) decay at x at some instant .
In the limit ,Equation (13) and therefore Equation (17) reduce to the asymptotic form
which is the probability density of a Laplace distribution  .
The general forms of the densities , , are shown as functions of displacement in Figure 1 for time units and diffusion and decay parameters , . is of Gaussian form; is approximately exponential with reflection symmetry across the plane . The total density takes on an intermediate form depending on the diffusion time. Figures 2-4 respectively show the individual variations in , , as functions of diffusion time. In particular, one sees how the total density evolves from Gaussian form to the mirror-reflected exponential form of a Laplace distribution as time increases from 10 units to infinity. The units of space and time in the figures are arbitrary; the numerical values were chosen to facilitate the subsequent Monte Carlo computer simulations and
Figure 1. Variation of probability densities (Survive), (Decay), (Total) as a function of displacement for t = 50 time units and parameters , .
Figure 2. Variation of survival density as a function of displacement at times (arbitrary units): (a) 10, (b) 50, (c) 100, (d) ¥. Parameters: , . The figures depict the spreading of a Gaussian distribution in time.
Figure 3. Variation of decay density as a function of displacement at times (arbitrary units): (a) 10, (b) 50, (c) 100, (d) ¥. Parameters: , . Figure (d) depicts a Laplace distribution.
Figure 4. Variation of total density as a function of displacement at times (arbitrary units): (a) 10, (b) 50, (c) 100, (d) ¥. Parameters: , .
graphical presentation. For example, given the designated value of λ, the probability of a single particle decay in time step is predicted from Equation (4) to be 1%. This is a convenient probability to use in the Monte Carlo simulations of Brownian motion presented in Section 4.
The first and second moments of the distributions (10), (13), and (17) are respectively calculated to be
from which it is seen from Equation (23) that the total variance of the Brownian displacements, derived from the solution to the 1D FPE (1), is exactly the same result (7) for obtained from the analytical solution of the 1D LE (3). The demonstration of this equality provides an operational interpretation of the difference between the two expressions for variance, Equations ((6) and (7)), which are displayed as functions of time in Figure 5. The plot designated “Survive” is the mean-square displacement (21) of particles that have not decayed by time t. As time increases, the displacement of surviving particles reaches a maximum and then decreases exponentially. Correspondingly, the plot designated “Decay”, is the mean-square displacement (22) of particles that have decayed at or before time t. The sum of the two plots generates the dashed plot designed “Total” given by Equation (23).
An experimental procedure to test the foregoing predictions would be to follow the Brownian diffusion of a sample of radioactive particles initially concentrated
Figure 5. Variation with time of the mean-square displacements , , . Parameters: , .
at the origin (in approximation to a delta function distribution) and record (a) the location of each particle as a function of time, as well as (b) the time at which any particle has decayed. An experiment of this kind is difficult to execute physically, but can be simulated on a fast computer by a Monte Carlo method. In general, the Monte Carlo method utilizes random numbers of a specified probability distribution to simulate the dynamics of a stochastic process by calculating and aggregating the outcomes of a large number of trials for a prescribed set of initial conditions. Commercially available symbolic mathematical software for science, engineering, and statistics usually includes random number generating algorithms for common probability distributions (e.g. binomial, Gaussian, Poisson, exponential) in their accompanying libraries of functions. The results of a Monte Carlo simulation of 1 million radioactively decaying particles undergoing 1D and 3D Brownian motion are discussed in Section 4.
It is also of interest to calculate the fourth moments of the displacement, since the variances of the mean-square distributions depend on them. Straightforward integration leads to
It is again to be noted that the total fourth moment (26), derived above directly from solution of the FPE, is the result obtained previously in Ref.  (Equation (73)) from the LE4.
Agreement of the FPE and LE predictions of the second and fourth moments demonstrates consistency of the two approaches at a practical level, since analysis of diffusion by Brownian motion rarely requires higher moments. However, for complete theoretical equivalence, which is essential if the two approaches to treating diffusion of radioactive matter are to be equally trusted, it is required to show that the FPE and LE yield equivalent statistical moments at all orders. This demonstration is given in the next section.
2.2. Equivalence of the Fokker-Planck and Langevin Analyses
The general problem of whether a distribution function is uniquely determined by the totality of its moments is referred to in theoretical statistics as the “problem of moments”  . This problem dates back to the end of the 19th century and has generated a copious mathematical literature whose content overall is beyond the scope of this paper. Of relevance to the present work, however, is the version of the problem associated with the name of Hamburger, who treated the case of random variables defined on the entire real axis. In brief, if the random variable has a finite moment-generating function (MGF), then the totality of the moments uniquely determines the distribution of that random variable  . This criterion is applicable here. In this section we show that the MGF derived from the Langevin equation (LE) for Brownian motion of radioactive particles leads to the identical sequence of moments derived from the probability density function of the corresponding Fokker-Planck equation (FPE). Since the sequence of moments is unique, the complete equivalence of the two approaches is established.
The MGF of a random variable Z is defined by the expectation value 
provided it exists5, where θ serves only as an expansion variable. A Taylor series expansion in θ of relation (27) leads to a superposition of the statistical moments
from which follows the operational expression by which individual moments can be calculated
for non-negative integer k. The term represents the completeness relation .
From the form of Equation (5) for the displacement variable X derived from the LE, it follows that the even statistical moments of X are given by the product of factors
, where the first factor is an expectation of the kth power of the random variable that characterizes the decay process, and the second factor is an expectation of the unit normal (Gaussian) distribution
that characterizes the spatial displacement of the particle during its survival. The function is the standard gamma function defined by the integral
One sees from Equation (31) that all odd statistical moments of spatial displacement vanish. This is a consequence of the left-right symmetry of Brownian motion in the absence of advection.
In Appendix 1 of Ref.  the MGF of the variable that occurs in the distribution (5) derived from the LE was shown to take the form
where is the survival probability per step defined in Equation (4), and the subscript n signifies a Brownian trajectory of n discrete steps each of duration in which t is the total diffusion time. Substitution of MGF (33) into Equation (29), followed by taking the continuum limit , leads to a closed form analytical expression for the moments of spatial displacement as a function of time
where the upper incomplete gamma function is defined by the integral 
An explicit listing of the five lowest non-vanishing moments shows an interesting series pattern arising primarily from the incomplete gamma function
In brief, apart from a numerical coefficient, the moments as a function of t take the form:
In the limit , the first term in Equation (34) vanishes exponentially, and the incomplete gamma function (35) likewise vanishes since the upper and lower limits to the integral become identical. Equation (34) then reduces asymptotically to
Figure 6 shows plots (solid curves) on a scale of the variation in moments (36) to (40) as a function of diffusion time. Dashed curves mark the asymptotic limits calculated from Equation (43).
We consider next the moments deducible from the Fokker-Planck equation by direct integration using the density function , Equation (17), which is the sum of the survival density , Equation (10), and the decay density , Equation (13). Evaluating the integrals, one obtains
Figure 6. Variation of moments (on scale) as a function of time t (solid curves) for powers k = (a) 2, (b) 4, (c) 6, (d) 8, (e) 10. Parameters: , . Dashed curves show the theoretical asymptotic limits.
is identical to Equation (34) derived from the LE moment-generating function.
Although the preceding demonstration is sufficient, it is useful to note that one can apply Equation (27) to obtain the moment-generating function directly from the probability densities derived from the FPE as follows
The integrals on the right side reduce to
which sum to the compact expression
in terms of the characteristic diffusion length z. Applying Equation (29) to the MGF (50) generates the spatial moments given by Equations (34) and (46). One final remark to avoid any confusion: Although the LE and FPE lead to identical statistical moments for all non-negative integers k, the moment-generating functions (50) and (33) are not the same because they pertain to two different random variables, X and , related by Equation (5).
We have therefore established that, despite the structural differences in form of the Fokker-Planck and Langevin equations derived in Ref.  , the two approaches to analyzing the diffusion of radioactively decaying particles yield statistically equivalent information.
2.3. 3D Brownian Motion
Since Brownian motions along the x, y, z axes are all uncorrelated, and since diffusion in the atmosphere without advection is isotropic, it is a straightforward matter to deduce the analogous 3D probability density of the radial displacement , by extension of the analysis used to obtain the 1D density function (17). The procedure involves two parts: 1) calculation of the survival density for the particle to reach a radial distance r from the origin at time t, and 2) calculation of the density for the particle to decay upon reaching a radial distance r from the origin at some moment in the time interval between 0 and t. The two sets of events are mutually exclusive and complementary (i.e. complete the sample space), whereupon the total density
is normalized to unity.
Consider first the density , which is defined by the differential relation
where numerical subscripts explicitly signify the dimensionality of a Brownian process. In the left side of relation (52) it is understood that the survival probability occurs as a factor only once (not three times) since the same time axis applies to motion along all three spatial dimensions. It then follows from Equation (10) that
which is equivalent to a chi-square distribution of k = 3 degrees of freedom. (The equivalence is established by the change of variable in the definition of the chi-square distribution where .)
Next, in analogy to the 1D relation (11), the density is obtained by integration of the differential relation
to yield the density
The integral in Equation (55) is a generalized incomplete gamma function  , defined by
The total density for undergoing a radial displacement r by Brownian motion in time t is then
and satisfies the normalization requirement for completeness
In the limit of infinite time, all particles will have decayed, and reduces to the asymptotic form
The function is the probability density of a gamma distribution of index 2  . Gamma distributions occur widely in statistical and nuclear physics; special cases include the chi-square distribution and the distribution of squares of Fourier amplitudes in a time series of radioactive decay  .
Figure 7 shows the spatial variation of densities (blue), (red), and (solid black) for diffusion time . As t increases, fewer particles survive and the peak of decreases toward the baseline. Correspondingly, the probability of decay increases and the peak of increases. In the limit of infinite time, the total density approaches the
Figure 7. Variation of 3D densities (a) , (b) , (c) (solid curves) as a function of radial displacement for diffusion time , and (d) (dashed curve). Parameters: , .
asymptotic density (dashed black) given by Equation (59). Although the peak of is lower than that of for finite t, its heavy tail falls off more slowly, and the area under both densities is 1.
The spatial moments of the foregoing 3D densities can be evaluated in closed form as given below (where ):
is known as the lower incomplete gamma function  . The complete kth moment of r is then
The lowest five moments (including the completeness relation ) calculated from Equation (63) are given in Table 1. Note that the complete moments (column 4) are expressible in terms of a single length, the characteristic diffusion length z of Equation (14). Plots (on a scale) of all the moments in Table 1 as a function of time look very much like the plots in Figure 6 and need not be shown. Instead, Figure 8 shows the variation in time of functions of the first two moments most often used to characterize the dispersal of a diffusing sample: (a) mean radial displacement (red curve), (b) root-mean-square displacement (blue curve), and (c) standard deviation about the mean radial displacement (green curve). The diffusion and decay parameters are
Table 1. Lowest radial moments of 3D Brownian trajectories.
Figure 8. Variation of 3D statistical measures of dispersion as a function of diffusion time: (a) mean , (b) root mean square , (c) standard deviation about the mean . To facilitate comparison at the same scale, the dimension of all three statistics is the first power of length. Parameters: , .
, . The asymptotic values are given by the relations in column 5 of Table 1.
An alternative approach, in analogy to Equation (50), is to determine the statistical moments by differentiation of the 3D moment-generating function (MGF) defined by the integral
which reduces to the expression
Although Equation (65) is complicated in appearance compared to relations (60)-(62), there are potential advantages to having the MGF at one’s disposal. In regard to the analytical or numerical calculation of moments, symbolic mathematical software can generally perform derivative operations more quickly and efficiently than integrations. And secondly, the MGF can be useful, as shown in the preceding sections, for deriving properties and establishing equivalence of statistical distributions.
3. Distribution of First-Passage Time (FPT)
The probability density for a decaying particle to diffuse in time t from the origin 0 to a location x between two absorbing boundaries at is derivable from the FPE by the method of images 
Although Equation (66) is not a closed form, since the summation includes an infinite number of terms, in practice the number of terms required for an accurate calculation can be relatively few depending on the values of the parameters and the duration t of the Brownian motion. For example, for the parameters , used in the preceding section, the density vanishes for at both absorbing boundaries for just the set of images (i.e. summation indices ), as shown in Figure 9. As the diffusion time increases, the peak of the density function drops toward the baseline (dashed black) and by (not shown in the figure) the left boundary is no longer maintained, whereupon one must sum the set of images (i.e. summation indices ), in Equation (66).
The survival function given explicitly in Equation (8) was obtained by integration of the density (66) over the range . As such,
Figure 9. Probability density (with ) of a decaying particle confined by absorbing boundaries to the interval . The diffusion times (arbitrary units) are: (a) 5, (b) 10, (c) 20, (d) 50, (e) 100. Parameters: , .
represents the probability that the diffusing particle has not reached boundary points or by time t―hence the usual term “survival function”. The probability density for termination of a Brownian trajectory is given by Equation (9) and can be expressed as a sum of two terms
The first function, expressed by (68), is interpretable as the actual FPT probability density, i.e. the density function for a particle to reach one of the two absorbing boundaries for the first time at t. In other words, describes a particle that has survived radioactive decay to reach a designated boundary, at which point the particle is removed from the system because the boundary is absorbing. We designate density (68) by the label “Absorption” in Figure 10, which shows the variation of (blue curve) as a function of t. A characteristic feature of density (68), which represents a special case of a Levy process  , is its heavy tail, as evidenced by the power-law dependence of the prefactor. However, the faster fall-off of the exponential decay factor dominates the slower fall-off of the power law, as is seen in Figure 10 by comparing the plot of with that of the corresponding Smirnov density  (dashed cyan curve) to which Equation (68) reduces for .
The second function expressed by (69) is the original survival function multiplied by the intrinsic nuclear decay constant λ. In other words, describes a particle that has decayed at time t before reaching a designated boundary. Bear in mind that it is probability (or, in practical terms, relative frequency) that is measurable, rather than probability density which is a mathematical function to be summed or integrated. Thus, to interpret the two expressions comprising density (67) one must actually examine the differential expression
Figure 10. Distribution of first-passage times of (a) particles that are absorbed at the boundaries; (b) particles that decay before reaching the boundary; (c) particles that are removed from the system because of either decay or absorption; (d) stable particles that reach the boundary [Smirnov]. Parameters: , . In the case of unstable particles the exponential decay function leads to a faster falloff in time than the power-law heavy tail.
, whereupon the factor is recognized from Equation (4) as the probability of radioactive decay. The plot of density (69) (red curve) in Figure 10 is therefore identified by the label “Decay”.
The sum of the two contributions in Equation (67), which together represent the loss of the particle from the system by either absorption or decay, is designated “Total” and shown (black curve) in Figure 10 as a function of time t. Thus, the process of radioactive decay, which introduces into the FPT diffusion problem an additional time parameter (the mean lifetime ), is consequential in several ways. First, it alters the form of the total density from that of the Smirnov density (for infinitely long-lived particles) in the regions of the time domain below the peak of the density and in the tail of the density. Second, this modification of form significantly affects the associated statistical moments. For example, the mean FPT for a non-decaying particle to reach a single absorbing boundary is infinite, whereas this moment always remains finite in the case of Brownian motion with decay  .
The mean time for a decaying particle to be removed from the system by reaching either of two symmetrically disposed absorbing boundaries at or by radioactive decay is the first moment of the total density (67)
where the characteristic diffusion length z is defined by relation (14). Two other ways of deriving relation (70) were shown in  ; these included (a) integration of the survival function and (b) solution of a screened Poisson differential equation.
In the limit of vanishing decay rate, the mean time (70) reduces to
which is the same mean time obtained for decaying particles in the limit . In the opposite limit , the mean first-passage time of a decaying particle asymptotically approaches the mean particle lifetime . This behavior is illustrated in Figure 11, which shows a plot of as a function of boundary location for 1) decaying particles (red curve) with and 2) stable particles (blue curve); in both cases the diffusivity is .
The mean time for just the process of particle absorption at the boundaries (blue curve in Figure 10)―in other words, for what strictly corresponds to the mean first-passage time―is the first moment of density defined in Equation (68). However, alone does not span a complete sample space of events because a particle can decay at time t without reaching a boundary; the mathematical consequence of incompleteness is that is not normalized to unity. The correctly normalized first moment is given by
Figure 11. Mean time (on scale) for (a) first passage to either of two absorbing boundaries or decay prior to first passage compared with (b) corresponding first passage of a stable particle. In the asymptotic limit the mean time of the unstable particle approaches the particle mean lifetime. Parameters: , .
where the characteristic diffusion velocity is defined by relation (15) and the normalization constant is
As in the case of the spatial distributions analyzed in the previous section, it is useful to have the second moment of the temporal distribution
Comparison of relations (74), (70), and (72) leads to the identity
which can be derived directly from the survival function (8), based on calculation of statistical moments by means of the cumulative distribution function rather than the probability density function; see Appendix 2 of Ref  for use of the cumulative distribution function to calculate statistical moments. Derivation of Equation (75) from the survival function is given in Appendix 1 of the present paper.
4. Monte Carlo Simulation of Brownian Motion of Radioactive Particles
The first application of a statistical method to study the dynamics of a physical system subject to numerous random interactions is attributable to Enrico Fermi’s unpublished investigations of neutron diffusion in the 1930s. Designated the Monte Carlo method in the 1940s by Los Alamos scientists Metropolis and Ulam  , and modernized in response to the development of programmable computers capable of generating large sets of pseudo-random numbers6, this approach now comprises the most widely used class of numerical methods to solve statistical problems in physics, chemistry, finance, and other fields  . Although there are numerous ways to implement a Monte Carlo algorithm, the basic core of the method is to: 1) define a domain of possible inputs, 2) generate these inputs randomly from an appropriate probability distribution, 3) substitute the inputs into the equations that determine the dynamics of the system being investigated, 4) extract from the aggregated results the statistical information of interest (such as probability densities, cumulative distribution functions, statistical moments, transition probabilities, and other statistical quantities).
To study the diffusion of a radioactive gas, we used a Monte Carlo method to create a large set of Brownian trajectories arising from the Langevin stochastic differential Equation (3) which incorporates two independent random number generators (RNG): 1) a Gaussian RNG to account for spatial displacement at each time step, and 2) a Bernoulli RNG to account for the possibility of decay at each step. The simulations comprised two parts: Part 1 was to determine the distribution of displacements as a function of time; Part 2 was to determine the distribution of first-passage times from the point of origin to one or the other of two symmetrically located absorbing boundaries at . Because of radioactive decay, it is possible that a particle will disintegrate before reaching any specified location. Under the conditions of the simulation, the Brownian trajectories in the first part were always finite since they terminated at a decay. The events in the second part fell into one of two categories: a) those for which first passage to a boundary preceded decay, and b) those for which decay eliminated the possibility of first passage to a boundary.
In implementation of Part 1, particles were sequentially created at the origin (a new particle being generated upon decay of a preceding one) and tracked in time as they underwent 3D Brownian motion in time steps (arbitrary units) for up to a maximum duration with . For the parameters used in these simulations, the probability of a particle surviving to was infinitesimal as shown in Appendix 2. At each time step, the displacement along each Cartesian axis was determined by three independent samples of a Gaussian RNG in accordance with Equation (3). The numbers of particles reaching locations at some time step within a specified interval for were then tallied and recorded in histograms (plots of frequency vs displacement) comprising bins (statistical categories) of width (same value for and ). and were selected to facilitate analysis and the graphical display of data, and could vary for different diffusion times .
Simulations were performed for different values of the diffusion time, diffusivity D, and radioactive decay rate λ. Histograms of 1D Brownian motion displayed in Figure 12 for values and were generated by projection of the 3D trajectories onto the x, y, z axes. The frequency scale of each histogram is in units of for comparison of the empirical histograms with the theoretical probability density (red curves) evaluated at the specified times. The progression of histogram panels A through D in Figure 12, corresponding to the increase in diffusion time from 10 to 1000 units, shows a nearly perfect visual match to the structure of the probability density function as it transforms from a Gaussian shape to the mirror-reflected exponential shape of a Laplace distribution. Histogram D, generated for , is indistinguishable at the scale shown from the asymptotic density for , as displayed by the red curve calculated from Equation (19). Figure 13 displays a histogram of the asymptotic 3D radial displacement superposed by the theoretical curve (red) calculated from Equation (63). Again, agreement between theory and simulation is in excellent visual accord. A comparison of 1D and 3D mean and root-mean-square displacements obtained by Monte-Carlo simulation and stochastic theory is given in Table 2 for specified diffusion times.
The three panels of Figure 14 display histograms, superposed respectively by the corresponding theoretical density functions (67)-(69), for (A) first-passage time to absorbing boundaries at , (B) time of particle loss by decay prior to reaching a boundary, and (C) total time of particle loss by either mechanism; simulation parameters are and . The Monte-Carlo simulation is again in excellent accord with theory and illustrates that the probability densities derived from the FPE account for the statistics of the process variable determined numerically from the LE.
Quantitative statistical tests, employing the ratio of observed to predicted frequencies as described in Appendix 3, showed for all histograms that deviations between computer-simulated outcomes and theoretical predictions could be attributable to pure chance.
The Fokker-Planck equation (FPE) and the Langevin equation (LE) provide different mathematical approaches to the analysis of Brownian motion of radioactive particles. Although the two equations differ structurally―the FPE is a multivariable partial differential equation, the LE takes the form of a stochastic differential equation based on Wiener and Bernoulli processes―we have demonstrated analytically that they are fully equivalent. The demonstration of complete equivalence was based on showing that the probability density functions derived from the FPE and the moment-generating function derived from the LE led to the identical complete sequence of statistical moments for 1D and 3D Brownian motion. Although such equivalence is expected for a continuous stochastic process that meets appropriate boundary conditions  , the nuclear transformation of particles is not a continuous process, and the equivalence of the two
Figure 12. Monte Carlo simulated distributions of radioactive particles with and undergoing 1D Brownian motion for diffusion times (in arbitrary units): (a) 10, (b) 50, (c) 100, (d) 1000. The horizontal axis is partitioned into bins of respective widths , , , . Histograms are normalized to unit area; the vertical axis records frequency in units of for comparison with the theoretical probability density (red curve) Equation (17). The envelope to Histogram (d) is density Equation (19) for .
Figure 13. Monte Carlo simulated radial distribution of radioactive particles with bins and bin width , undergoing 3D Brownian motion for time (arbitrary units). Parameters: , . The envelope (red curve) is the theoretical density, Equation (59) for .
Table 2. Comparison of theoretical and monte-carlo simulated Brownian statistical moments for D = 0.5, λ = 0.01.
Figure 14. Monte Carlo simulated distributions of the times for (a) first passage to absorbing boundaries located at (arbitrary units), (b) decay prior to first passage, and (c) combined processes of particle removal. Sample sizes for first passage and decay are respectively , . Parameters: , . Red curves are the respective theoretical densities , , and .
methods was not a priori evident or demonstrated at the time of publication of Ref.  in which the Fokker-Planck and Langevin equations for radioactive systems were initially constructed and investigated.
Besides the analytical demonstration, we have tested the two approaches by means of a Monte Carlo method that generated histograms of Brownian trajectories from the LE whose empirical profiles were matched nearly perfectly by the probability density functions derived from the FPE. Statistical tests applied to the sets of sufficiently occupied bins showed that any deviations between theory and observation could be attributable to pure chance.
In the course of these demonstrations, we have derived and reported in this paper explicit mathematical expressions for the probability density functions, spatial moments, and moment-generating functions of 1D and 3D Brownian trajectories of radioactive particles. We have also derived an explicit expression for the probability density function for the distribution of first-passage times of a radioactively decaying particle to diffuse to specified absorbing boundaries. These functions, which should prove useful in the measurement and analysis of radioactive particles dispersed in the environment or investigated in the laboratory, differ markedly from corresponding statistical functions for Brownian motion of non-decaying particles.
In statistical terminology, radioactive decay as described by the Rutherford-Soddy law, whereby the intrinsic decay rate is constant and independent of interactions external to the nucleus, is one of the simplest examples of a birth-and-death process of which there exists an extensive literature  . The excellent agreement between theory and Monte-Carlo simulated results reported in this paper confirms the validity of treating radioactive decay as a continuous exponential process in the Fokker-Planck equation and as a discrete Bernoulli process in the associated Langevin equation. However, looking beyond the physics relating to measurement and analysis of radioactivity, we believe that the approaches taken in this paper and the antecedent reference  will serve as useful guides in the construction, solution, and interpretation of stochastic differential equations for more complex birth-and-death processes.
Relation between First and Second Moments of the First-Passage Time
In Appendix 2 of Reference  it is shown that the first and second moments of a random variable T defined over the non-negative domain take the form
in which the cumulative distribution function , defined in terms of the probability density function , is
By definition, the survival function is the probability that the diffusing particle has not reached a boundary―or, in other words, the probability that the random variable T, the first-passage time (FPT), exceeds t. Therefore
Substitution of Equation (79) into Equations (67), (76), (77), with recognition that Equation (69) has the form
which is equivalent to Equation (75).
Probability of Survival: The Geometric Distribution
Let p be the probability of surviving a single time step and the corresponding probability of decay. If survival or decay at each time step is independent of the outcome at any other time step, then the probability of surviving time steps and decaying on the kth is given by the geometric probability law
The first few moments of the distribution are
Mean number of steps:
Mean square number of steps:
Given probabilities of survival and decay to be and respectively at each time step, the probability of surviving 5,000 time steps is
and the number of surviving particles out of a sample of size of 106 is then or effectively 0. The mean number of time steps to a decay event is . The standard deviation in the number of time steps to a decay event is .
Statistical Tests of the Monte-Carlo Histograms
Operationally, a histogram is a plot of event frequency as a function of event class. Mathematically, a histogram is a multinomial distribution of outcomes (classes or bins) with set of frequencies and probabilities governed by the discrete probability function (Ref.  , 10-12, 154-156)
subject to the constraint
For a radioactive particle undergoing Brownian motion, the probability that the particle is located between coordinates and that defines the kth bin of constant width is obtained from Equation (17)
The mean frequency of the kth bin is then predicted to be
where n is the total sample size (i.e. number of particles). For a multinomial distribution, the variance of is given by
For large sample size and numerous bins the transformation of observed frequencies
is by the Central Limit Theorem (Ref  , 28-32) an approximate realization of a unit Gaussian distribution. Dividing both sides of Equation by and rearranging the expression leads to the ratio of the observed frequency to the corresponding theoretical mean
which is distributed as a Gaussian of mean 1 and variance . Equation (93) follows from the identity (Ref  , 18)
Summing Equation (93) over all K classes generates the ensemble random variable R
It then follows from relations (93) and (95) that the probability of obtaining under repeated trials a value greater or equal to the observed value for the kth bin, or an ensemble value R greater or equal to the observed value is given by
The numerical values associated with integrals (97) and (98) are generally referred to as P-values, in which for some specified cumulative distribution function .
Consider, as an example, histogram D of Figure 12 for Brownian motion of a radioactive particle over a time period of units, which is effectively in the asymptotic limit . The ratios for sufficiently occupied bins, with calculated from Equations (90) and (89) with given by Equation (19), all varied by small amounts about unity. The mean ratio for bins with more than 5 events (the usual statistical cutoff) was . The statistics for the ensemble ratio test were , , leading to . The statistical threshold below which to reject the hypothesis that the theoretical profile fits the histogram is ordinarily 0.05.
Up to this stage, the preceding test is similar in procedure to that of a standard Pearson chi-square test, except that it involves the ratio of observed to predicted frequencies, whereas a chi-square test involves the square of frequency differences. The chi-square test, although widely used, is based on a chain of assumptions that do not necessarily apply to Brownian motion (e.g. that the probability of an outcome is approximately Poissonian; see Ref.  , 61-71) and has long been known to do poorly―i.e. to reject a prevailing hypothesis that is in fact true―under conditions of large sample size such as pertains here7.
Furthermore, to discern whether an empirical histogram is a good representation of the overall profile of a proposed probability density function, it is of utility to go beyond the ensemble P-value and to calculate the P-values (97) individually for all sufficiently occupied bins. In the above example for histogram D of Figure 12, only 2 of the 37 calculated P-values fell below the 5% threshold. However, the P-value is itself a realization of a random variable governed by a uniform distribution over the range [0, 1] with mean and standard deviation (Ref.  , 34-38). Thus, one would expect 5% of 37―i.e. approximately 2 outcomes―to fall below the threshold on the basis of chance alone. This is what was observed. Moreover, the observed mean of the 37 bins was 0.594, which falls within the range . Finally, it was also observed from a scatter plot that the 37 P-values were distributed randomly between 0 and 1, as is expected for a uniformly distributed random variable .
Statistical tests for goodness-of-fit can never prove that a hypothesized theory is the “true” explanation of some set of empirically derived results. What the preceding statistical tests show, however, is that they provide no statistical grounds for rejecting the conclusion, reached by rigorous analysis, that the Monte-Carlo simulations of the Brownian motion of radioactive particles, based on the Langevin Equation (3), are accounted for by the theoretical probabilities derived from the Fokker-Planck Equation (1).
1It is standard notation in statistical physics to use an upper-case letter to represent a random variable and a corresponding lower-case letter to represent a realization (i.e. sample) of that variable.
2The symbol stands for a normal (Gaussian) distribution of mean μ and variance σ2.
3The symbol stands for a binomial distribution of n trials with probability of success p. A Bernoulli distribution is the special case .
4Regrettably, due to a typographical error in Equation (73) of Ref.  , the printed 4th moment was a factor 3 smaller than expression (26) above. A generalization of this calculation is given in Section 2.2.
5Some distributions, such as the Cauchy distribution (“Lorentzian line shape” in physics) do not have a MGF. In such cases, one can calculate statistical moments from the characteristic function (CF) obtained by substituting iθ for θ in Equation (27).
6Pseudo-random numbers generated by a computer may satisfactorily pass a battery of statistical tests for randomness, but are not truly random because the generating algorithm produces the same output sequence of numbers whenever initiated by the same seed value.
7J. Berkson, Journal of the American Statistical Association 33 (1938) 526-536.