Brownian Motion of Radioactive Particles : Derivation and Monte Carlo Test of Spatial and Temporal Distributions

Stochastic processes such as diffusion can be analyzed by means of a partial differential equation of the Fokker-Planck type (FPE), which yields a transition probability density, or by a stochastic differential equation of the Langevin type (LE), which yields the time evolution of a statistical process variable. Provided the stochastic process is continuous and certain boundary conditions are met, the two approaches yield equivalent information. However, Brownian motion of radioactively decaying particles is not a continuous process because the Brownian trajectories abruptly terminate when the particle decays. Recent analysis of the Brownian motion of decaying particles by both approaches has led to different mean-square displacements. In this paper, we demonstrate the complete equivalence of the two approaches by 1) showing quantitatively and operationally how the probability densities and statistical moments predicted by the FPE and LE relate to one another, 2) verifying that both approaches lead to identical statistical moments at all orders, and 3) confirming that the analytical solution to the FPE accurately describes the Brownian trajectories obtained by Monte Carlo simulations based on the LE. The analysis in this paper addresses both the spatial distribution of the particles (i.e. the question of displacement as a function of diffusion time) and the temporal distribution (i.e. the question of first-passage time to fixed absorbing boundaries).


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 222 Rn and thoron 220 Rn, or are associated with power reactor releases such as iodine 129 I and 131 I, cesium 131 Cs and strontium 90 Sr, or are produced in nuclear weapons manufacture and testing such as tritium 3 H, among many other possibilities [1].
Although diffusion as a physical process has been investigated theoretically since Einstein's seminal work on Brownian motion in 1905 [2], 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, [3] [4] [5] 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 [6] leads to the transition probability density ( ) 0 0 , , p t t x x for finding a particle at location x at time t, given that it was at 0 x at time 0 t . 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 0 0 x = at 0 0 t = 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 [6] takes the form of an update stochastic differential equation [7] [8] based on a Wiener process [9] [10] whose graphical solution shows the sequential displacements ( ) x t of the particle as a function of time and whose analytical solution yields the statistical distribution of the displacement variable X. 1 The symbol t n , which governs particle displacement, represents a random sample from the unit normal distribu- where the subscript t and superscript explicitly denote the temporal range with respect to which t n is associated. In other words, two in which the probability of radioactive decay within the differentially small time interval dt is given by dt λ , the transformation law first recognized by Rutherford and Soddy and used by von Schweidler to derive the law of exponential decay [11]. Thus, s p 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 [6] to be Note that the expression in (5) is not a function, but a statistical distribution (Gaussian) whose variance 2 2 d k D t Σ after k particle displacements is expressed in terms of another random variable 2 k Σ that was derived and explained in [6] (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 [12] [13]. It was shown theoretically in [6], 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 X (which is 0 in the absence of advection) and mean-square displacement 2 X . FPE solution (2) and LE solution (5) (for 1D Brownian motion in the continuum limit) led to expressions It 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. 2 The symbol ( ) 2 , N µ σ stands for a normal (Gaussian) distribution of mean μ and variance σ 2 . 3 The symbol ( ) , B n p stands for a binomial distribution of n trials with probability of success p. A Bernoulli distribution is the special case 1 n = . x X X ∆ = − . 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 2 k Σ 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.

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 [14]. Such processes are relevant to the quantitative measurement of radioactivity by methods that sample a radioactive source by means of diffusion [15] [16] [17], to studies of the dispersal of radioactive contaminants in the atmosphere [18], and to determinations of radiation exposure to critical organs such as the lungs [19].
In [6] the survival function ( ) S t 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 [6] in the case of two absorbing boundaries at b is the probability that the particle, initially located at 0 0 x = , 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 ( ) x is 0 for all values of time t. Equation (8) was solved by a method of images [20], 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 0 t ∞ ≥ ≥ .
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 differen-M. P. Silverman, A. Mudvari tiation of the survival function [21] ( ) ( ) 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. [6] 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 ( ) T p t was not given. The density ( ) T p t 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.

1D Brownian Motion
As discussed in [6], 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 e t λ − , the probability that the particle has decayed becomes 100% only at t = ∞ .
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 [22] 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 0 0 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 S w for the Brownian particle to diffuse to x and survive-hence the subscript "S".
Multiplying FPE solution (2) by dt λ 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 D w in Equation (11) signifies that the particle has decayed at x at some time between t and d t t + . Thus, the probability density that a radioactive particle decayed at x irrespective of the moment of decay within the finite interval [ ] 0,t is obtained by integration leading to ( ) where D ζ λ ≡ (14) is the characteristic diffusion length and is the characteristic diffusion velocity (See Appendix 1 of Ref. [15].).
Probability densities (10) and (13)  1 exp 4 e 4π is normalized to unity ( ) In the limit t → ∞ , Equation (13) and therefore Equation (17) reduce to the asymptotic form which is the probability density of a Laplace distribution [23].
The general forms of the densities     graphical presentation. For example, given the designated value of λ, the probability of a single particle decay in time step d 1 t = 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 ( (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 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. [6] (Equation (73)) from the LE 4 .
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.

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" [24]. This problem dates back to the end of the 19 th 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 4 Regrettably, due to a typographical error in Equation (73) of Ref. [6], the printed 4 th moment was a factor 3 smaller than expression (26) provided it exists 5 , where θ serves only as an expansion variable. A Taylor series expansion in θ of relation (27) leads to a superposition of the statistical moments (28) from which follows the operational expression by which individual moments can be calculated for non-negative integer k. The term 0 k = 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 k th power of the random variable 2 n Σ 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 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. 5 Some 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).

M. P. Silverman, A. Mudvari
In Appendix 1 of Ref. [6] the MGF of the variable 2 n Σ that occurs in the distribution (5) derived from the LE was shown to take the form where S 1 d p t λ = − 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 dt t n = in which t is the total diffusion time. Substitution of MGF (33) into Equation (29), followed by taking the continuum limit ( ) n → ∞ , leads to a closed form analytical expression for the moments of spatial displacement as a function of time where the upper incomplete gamma function ( ) 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 t → ∞ , the first term in Equation (34) whose sum Dashed curves show the theoretical asymptotic ( ) 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 ( ) X g θ directly from the probability densities derived from the FPE as follows The integrals on the right side reduce to  (50) and (33) are not the same because they pertain to two different random variables, X and 2 k Σ , 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. [6], the two approaches to analyzing the diffusion of radioactively decaying particles yield statistically equivalent information.

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 , w r t for the particle to reach a radial distance r from the origin at time t, and 2) calculation of the density ( ) 3D , w r t 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 ( )

3S
, w r t , which is defined by the differential rela- which is equivalent to a chi-square distribution 2 k χ of k = 3 degrees of freedom.
(The equivalence is established by the change of variable The integral in Equation (55) is a generalized incomplete gamma function [28], defined by ( ) The total density for undergoing a radial displacement r by Brownian motion in time t is then in which is known as the lower incomplete gamma function [27]. The complete k th moment of r is then The lowest five moments (including the completeness relation 0 k = ) 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 ζ of Equation (14). Plots (on a 10 log scale) of all the moments in Table 1 M. P. Silverman, A. Mudvari 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 t r (red curve), (b) root-mean-square displacement      Table 1.
An alternative approach, in analogy to Equation (50), is to determine the statistical moments by differentiation of the 3D moment-generating function

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 b x = ± is derivable from the FPE by the method of images [6] ( ) Although Equation (66)  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 ab-M. P. Silverman, A. Mudvari sorbing boundaries for the first time at t. In other words, ( ) ( ) 0 T p t 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 ( ) ( ) 0 T p t (blue curve) as a function of t. A characteristic feature of density (68), which represents a special case of a Levy process [30], is its heavy tail, as evidenced by the 3 2 t − 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 ( ) ( ) 0 T p t with that of the corresponding Smirnov density [31] (dashed cyan curve) to which Equation (68) reduces for 0 The second function expressed by (69) is the original survival function ( ) S t multiplied by the intrinsic nuclear decay constant λ. In other words, ( ) ( ) 1 T p t 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 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  The mean time T 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 ζ is defined by relation (14). Two other ways of deriving relation (70) were shown in [6]; these included (a) integration of the survival function (8) 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 ( ) 1 ζ   . In the opposite limit ( ) 1 ζ   , the mean first-passage time of a decaying particle asymptotically approaches the mean particle lifetime 1 λ − . This behavior is illustrated in Figure 11, which shows a plot of ( ) log T as a function of boundary location  for 1) decaying particles (red curve) with 0.01 λ = and 2) stable particles (blue curve); in both cases the diffusivity is The mean time 0 T 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  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 d v 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)

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 Fer- Ulam [32], and modernized in response to the development of programmable computers capable of generating large sets of pseudo-random numbers 6 , this approach now comprises the most widely used class of numerical methods to solve statistical problems in physics, chemistry, finance, and other fields [33].
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 0 0 x = to one or the other of two symmetrically located absorbing boundaries at b x = ± . 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, 6 10 p N = 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 1 t ∆ = 6 Pseudo-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.
(arbitrary units) for up to a maximum duration max T K t = ∆ with 5000 K = . For the parameters used in these simulations, the probability of a particle surviving to max T was infinitesimal as shown in Appendix 2. At each time step, the displacement along each Cartesian axis was determined by three independent samples ( ) , , , , x y z at some time step within a specified interval k t k t = ∆ for 0 K k ≥ ≥ were then tallied and recorded in histograms (plots of frequency vs displacement) comprising b N bins (statistical categories) of width x ∆ (same value for y ∆ and z ∆ ). b N and x ∆ were selected to facilitate analysis and the graphical display of data, and could vary for different diffusion times k t .
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 0.5 D = and 0.01 λ = 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 (13) (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 1000 t = , is indistinguishable at the scale shown from the asymptotic density for t = ∞ , as displayed by the red curve calculated from Equation (19). Figure 13 displays a histogram of the asymptotic 3D radial displacement r 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 10 b x = ± , (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 0.5 D = and 0.01 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 ( ) , , , X Y Z R 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.

Conclusions
The Fokker-Planck equation (FPE) and the Langevin equation (LE) provide dif-M. P. Silverman, A. Mudvari ferent 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 [34], the nuclear transformation of particles is not a continuous process, and the equivalence of the two    Ref. [6] 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 [35]. 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 [6] will serve as useful guides in the construction, solution, and interpretation of stochastic differential equations for

Relation between First and Second Moments of the First-Passage Time
In Appendix 2 of Reference [16] 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 ( ) F t , defined in terms of the probability density function By definition, the survival function ( ) S t 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) The numerical values associated with integrals (97) and (98) are generally referred to as P-values, in which ( ) ( ) Consider, as an example, histogram D of Figure 12 for Brownian motion of a radioactive particle over a time period of 1000 t = units, which is effectively in the asymptotic limit t → ∞ . The ratios k k n m for 37 K = sufficiently occupied bins, with k m calculated from Equations (90) and (89) with ( ) , w x t 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 1.001 0.091 ± . The statistics for the ensemble ratio test were obs 37.05 R = , 2 0.67 S = , leading to ( ) obs Pr 0.47 R R ≥ = . The statistical threshold below which to reject the hypothesis that the theoretical profile fits the histogram is ordinarily 0.05.

M. P. Silverman, A. Mudvari
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. [26], 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 here 7 .
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 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 P of the 37 bins was 0.594, which falls within the range P P σ ± . 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).