Dynamics Underlying the Gaussian Distribution of the Classical Harmonic Oscillator in Zero-Point Radiation

In the past decades, Random Electrodynamics (also called Stochastic Electrodynamics) has been used to study the classical harmonic oscillator immersed in the classical electromagnetic zero-point radiation. Random Electrodynamics (RED) predicts an identical probability distribution for the harmonic oscillator compared to the quantum mechanical prediction for the ground state. Moreover, the Heisenberg minimum uncertainty relation is also recovered with RED. To understand the dynamics that gives rise to this probability distribution, we perform an RED simulation and follow the motion of the oscillator. This simulation provides insight in the relation between the striking different double-peak probability distribution of the classical harmonic oscillator and the Gaussian probability distribution of the RED harmonic oscillator. A main objective for RED research is to establish to what extent the results of quantum mechanics can be obtained. The present simulation method can be applied to other physical systems, and it may assist in evaluating the validity range of RED.


I. INTRODUCTION
According to quantum electrodynamics, the vacuum is not a tranquil place. A background electromagnetic field called the zero-point field, or vacuum field, is always present, independent of any external electromagnetic source [1]. This has inspired an interesting modification to classical mechanics, called random electrodynamics (RED) [2]. As a variation of classical electrodynamics, RED adds a classical zero-point electromagnetic field with random phases to the original theory. In this paper, the name "RED vacuum field" will be used for this classical electromagnetic field. The RED vacuum field has some properties identical to the QED vacuum field 1 such as the Lorentz-invariant spectral energy density [1,2]. With the aid of this background field, RED is able to reproduce a number of results that were originally thought to be pure quantum effects [1,2,[4][5][6][7]. In particular, Boyer has shown that the moments x n of a harmonic oscillator immersed in the RED vacuum field, called the RED harmonic oscillator, are identical to those of the quantum harmonic oscillator in the ground state [8]. As a consequence, the Heisenberg minimum uncertainty relation is satisfied for an RED harmonic oscillator, and the probability distributions of the RED harmonic oscillator is also the same as that of the ground state quantum harmonic oscillator.
Although classical mechanics and RED are both theories that gives trajectories of particles, the position probability distributions of the harmonic oscillator in both * email: u910134@alumni.nthu.edu.tw † email: hbatelaan2@unlnotes.unl.edu 1 However, the RED vacuum field does not fully correspond to the vacuum field derived in quantum electrodynamics, as discussed in [1,3].
theories are dramatically different. In classical mechanics, it is most likely to find the particle at the two turning points of the trajectory. Consequently, the probability distribution has two peaks as shown in Fig. 1. On the other hand, the probability distribution for an RED harmonic oscillator is a Gaussian distribution. What happens to the dynamics of the classical harmonic oscillator in the presence of the RED vacuum field such that its probability distributions would become Gaussian? Why do the widths of these distributions satisfy Heisenberg's minimum uncertainty relation? In this paper, we perform an RED simulation that answers these questions. Most of the results obtained with the RED theory reproduce quantum mechanical results analytically. Whereas Cole has a series of numerical studies on the hydrogen atom [9][10][11][12], numerical studies in this field are rare. The major challenge for RED simulation is to incorporate the RED vacuum field which has an unbounded spectrum. A representative selection of vacuum field modes is thus the key to successful RED simulations. In this study, our approach for the selection of modes has been extensively tested for convergence, and we document the details of our simulation so it can be used by others.
The organization of this paper is the following. First, in Sec. II Boyer's results on the RED harmonic oscillator are briefly reviewed. Based on these results, the probability distribution for the RED harmonic oscillator is derived. Second, in Sec. III a detailed numerical method for simulating the RED vacuum field and the RED harmonic oscillator is presented. Third, in Sec. IV the trajectory of the RED harmonic oscillator is solved numerically, and the constructed probability distribution is compared to the analytical probability distribution. Two sampling methods are used in constructing the probability distribution from the simulated trajectories. The first approach is "sequential sampling", which is suitable for studying the relation between the dynamics of Top: Without any external force, the classic harmonic oscillator that is initially displaced from equilibrium performs a simple harmonic oscillation with constant oscillation amplitude. The resulting probability distribution has peaks at the two turning points. Bottom: Without any external force except for the RED vacuum field, the RED harmonic oscillator undergoes a motion that results in a Gaussian probability distribution. This motion and its relation to its classical counterpart is investigated in detail in this paper.
the RED harmonic oscillator and its probability distribution. The second approach is "ensemble sampling", which lends itself well for parallel computing and is convenient for statistical interpretations. Lastly, in Sec. VI we discuss some potential applications of RED simulation in studying other quantum phenomena, including the electron doubleslit diffraction. Numerical studies have an advantage over the analytical solutions in that they can be adopted to a range of physical systems, and we hope that the current method of simulation may assist in assessing RED's validity range.

A. Brief Review of Boyer's Work
In his 1975 papers [2,8], Boyer calculated the statistical features of an RED harmonic oscillator, and the Heisenberg minimum uncertainty relation is shown to be satisfied for such an oscillator. The RED vacuum field used in Boyer's work arises from the homogeneous solution of Maxwell's equations, which is chosen to be zero in classical electrodynamics [2]. In unbounded space, the RED vacuum field has an integral form 2 , where ω = c|k|, andθ(k, λ) is the random phase uniformly distributed in [0, 2π]. The integral is to be taken over all k-space. The two unit vectors, ε(k, 1) and ε(k, 2), describe a polarization basis in a plane that is perpendicular to the wave vector k, Furthermore, the polarization basis vectors are chosen to be mutually orthogonal, To investigate the dynamics of the RED harmonic oscillator, Boyer used the dipole approximation, to remove the spatial dependence of the RED vacuum field, Eq. (1), from the equation of motion. Therefore, the equation of motion for an RED harmonic oscillator used in Boyer's analysis is where m is the mass, e is the charge, ω 0 is the natural frequency, and Γ is the radiation damping parameter. The x-component of the RED vacuum field in Eq. (5) is and the steady-state solution is obtained as where C(k, λ) ≡ −ω 2 + ω 2 0 − iΓω 3 . Additionally, using the condition of sharp resonance, Boyer further calculated the standard deviation of position and momentum from Eq. (7) by averaging over the random phaseθ [8], The above result satisfies the Heisenberg minimum uncertainty relation, From an energy argument, Boyer showed that this uncertainty relation can be derived from a delicate balance between the gain of energy from RED vacuum field driving and the loss of energy through radiation damping [2].

B. Position Probability Distribution
Given the knowledge of the moments x n θ , the Fourier coefficients Fθ(k) of the position probability distribution 3 3 An RED probability distribution can be constructed from either data collected at a sequence of time, or data drawn from an RED Pθ(x) can be determined through the relation Using Eq. (7) and the relation from Boyer's paper [8] the moments x n θ can be evaluated, where m is a natural number. Consequently, only evenpower terms are contributing in Eq. (12), and the Fourier coefficients Fθ(k) are determined, Therefore, although not explicitly given, it is implied by Boyer's work [8] that the position probability distribution of the RED harmonic oscillator is which is identical to the position probability distribution of the quantum harmonic oscillator in the ground state 4 .
ensemble. In the RED ensemble, each member is characterized by the random phaseθ of the RED vacuum field. In this paper, the probability distribution constructed from the ensemble sampling method is denoted as Pθ(x). 4 This result is consistent with the phase space probability distribution given in Marshall's work [13].

A. RED Vacuum Field in Bounded Space
A field confined in a space of volume V with zero value boundary condition has a discrete spectrum, and a summation over the wave vectors k is required. On the other hand, the field in unbounded space is not subject to any boundary condition, so every wave vector k is allowed. Thus, the field in unbounded space takes an integral form [2], while the field in bounded space takes a summation form [1,3]. In a simulation, it is convenient to write the RED vacuum field in its summation form, where ω = c|k|,θ k,λ is the random phase uniformly distributed in [0, 2π], and V is the volume of the bounded space. A derivation of the summation form of the RED vacuum field in bounded space is given in Appendix A 2.
Since the range of the allowed wave vectors k is over all k-space and the RED vacuum field diverges for infinitely large |k|, we choose to sample only the wave vectors k whose frequencies is within the finite range [ω 0 − ∆/2, ω 0 + ∆/2]. Such sampling is valid as long as the chosen frequency range ∆ completely covers the characteristic resonance width Γω 2 0 of the harmonic oscillator, On the other hand, the distribution of the allowed wave vectors k depends on the specific shape of the bounded space. In a cubic space of volume V , the allowed wave vectors k are uniformly distributed at cubic grid points, and the corresponding RED vacuum field is The sampling density is uniform and has a simple relation with the space volume V , Nevertheless, such uniform cubic sampling is not convenient for describing a frequency spectrum. In order to sample only the wave vectors k in the resonance region, spherical coordinates are used. In addition, for the sampling to be uniform, each sampled wave vector k must occupy the same size of volume element ∆ 3 k = k 2 sin θ∆k∆θ∆φ. To satisfy both conditions, we use a set of specifically chosen numbers (κ ijn , ϑ ijn , ϕ ijn ) to sample the wave vectors k. Namely, for i = 1 . . . N κ , j = 1 . . . N ϑ , and n = 1 . . . N ϕ , where R (0) is a random number uniformly distributed in [0, 2π], and the stepsizes are constant, Such a set of numbers (κ ijn , ϑ ijn , ϕ ijn ) is then used for assigning the spherical coordinates to each sampled wave vector k, Consequently, each sampled wave vector k ijn will be in the resonance region and occupy the same size of volume element, Under the uniform spherical sampling method (described by Eqs. (22), (26), and (27)), the expression for the RED vacuum field, Eq. (18), becomes where k = k ijn . It is worth noticing that when the total number of wave vectors N k becomes very large, both uniform spherical and cubic sampling effectively sample all the wave vectors k in k-space. Therefore, in the limit of large sampling number N k → ∞, the two sampling methods become equivalent 5 , and Eq. (21) can be used for both sampling methods to calculate the volume factor V in Eqs. (20) and (28), where V k = 4π 3 In a computer simulation, the summation indices in Eq. (28) can be rewritten as where the multiple sums indicate a numerical nested loop, and the wave vector k = k ijn is chosen according to the uniform spherical sampling method. In practice, N ϑ and N ϕ need to be sufficiently large so that the wave vector k at a fixed frequency may be sampled isotropically. Also, a large N κ is required for representative samplings in frequency. As a result, N k = N κ N ϑ N ϕ has to be very large to achieve convergence when using uniform spherical sampling. Therefore, sampling k for many angles at each frequency makes RED simulation computationally intensive. To improve the efficiency of our simulation, we sample k for only one random angle (θ i , φ i ) at each frequency (i.e. N κ = N ω , N ϑ = 1, and N ϕ = 1). Namely, N k = N ω , and for i = 1 . . . N ω , where and i .
the limit of large sampling number (i.e. N k → ∞) is equivalent to the limit of unbounded space (i.e. V → ∞). At this limit, the volume element becomes differential (denoted as d 3 k) and is free from any specific shape associated with the space boundary. Therefore, all sampling methods for the allowed wave vectors k become equivalent, and the summation form approaches the integral form. This is consistent with the fact that no volume factor V is involved in the integral form of the RED vacuum field in unbounded space, as shown in Eq. (1).
The stepsize ∆κ is specified in Eq. (23), ϑ = R (1) is a random number uniformly distributed in [−1, 1], and ϕ = R (2) is another random number uniformly distributed in [0, 2π]. As the number of sampled frequencies N ω becomes sufficiently large, the random angles (θ, φ) will be uniformly distributed on a unit sphere, and the angular distribution of the wave vectors k is isotropic.
In the limit ∆ ω 0 1, the above sampling method (described by Eqs. (32), (33), and (34)) and the uniform spherical sampling method both approach a uniform sampling on a spherical surface at the radius r k = ω 0 c . In this limit, the addition of the condition in Eq. (19) leads to the possible choices of the frequency range ∆, Within this range (Eq. (35)), the expression for the RED vacuum field in Eq. (18) becomes where k = k i . For large N k = N ω , the volume factor V is calculated using Eq. (29). Finally, for a complete specification of the simulated RED vacuum field, Eq. (36), the polarizations ε k,λ need to be chosen. From Eq. (29), we notice that large N k gives large V . Since for large V the RED vacuum field is not affected by the space boundary, there is no preferential polarization direction. Therefore, the polarizations are isotropically distributed. The construction for isotropically distributed polarizations is discussed in detail in Appendix B. Here we give the result that satisfies the property of isotropy and the properties of polarization (described by Eqs. (2) and (3)), where χ is a random number uniformly distributed in [0, 2π]. With the wave vectors k (described by Eqs. (32), (33), and (34)) and the polarizations ε k,λ (described by Eq. (37)), the endpoints of the sampled RED vacuum field vector are plotted on a unit sphere as shown in Fig. 2, which illustrates the isotropy of the distribution. In summary, the RED vacuum field, Eq. (36), can be specified by a set of four numbers (κ i , ϑ i , ϕ i , χ i ), chosen according to Eqs. (32), (33), (34), and Eq. (37). The only assumption used in determining these numbers is Eq. (35), which is equivalent to the sharp resonance condition, Eq. (8), used in Boyer's analysis.

B. Simplified Equation of Motion
An RED harmonic oscillator is a charged classical harmonic oscillator immersed in the RED vacuum field. The complete equation of motion is where m is the mass, e is the charge, and ω 0 is the natural frequency. The radiation damping parameter . Only the motion in the x-direction is considered in this paper, so the position of the particle is r = (x, 0, 0), and the velocity of the particle Three terms are involved in the equation of motion -the spring force, radiation damping, and the RED vacuum field, In unbounded space, the RED vacuum field is given by Eq. (1). In a bounded space of volume V , the RED vacuum field, including both the electric part E vac as in Eq. (18) and the magnetic part B vac , is where ω = c|k|, andθ k,λ is the random phase uniformly distributed in [0, 2π]. Throughout this section the RED vacuum field as described by Eq. (42) will be used. Before a simulation of the RED harmonic oscillator is performed, we may increase the efficiency of our computation by making some approximations to the equation of motion, Eq. (38). Namely, the terms F damp , F vac can be approximated by using the dipole approximation k · r 1 (Eq. (4)), and the sharp resonance condition Γω 0 1 (Eq. (8)). Denoting the order of magnitude for length, time, and electric field as x 0 , t 0 , and E 0 respectively, the order of magnitude of the two terms of F vac in Eq. (41) are The sharp resonance condition Eq. (8) makes the harmonic oscillator responsive to only frequencies that are close to ω 0 , thus the time scale of the particle's motion can thus be evaluated, Furthermore, the dipole approximation from Eq. (4) implies which can be used to approximate F vac with its dominating term, The dipole approximation also removes the spatial de- 6 The order of magnitude for dimension A is denoted as [A].
x in Eq. (48) is well known to give runaway solutions. Boyer's approach is to find the solution through Fourier transformation [8]. Following this approach, the steady-state solution for Eq. (48) may be obtained, (49) could be evaluated numerically, this would not allow for straightforward generalization to physical systems where the equation of motion is modified, which is one of the purposes of this work. Hence, we follow the approach in [14][15][16] to avoid runaway solutions. According to classical electrodynamics, the equation of motion for an electron with radiation damping is which can be approximated with using the assumption mγ ... x F ext . Such an assumption is necessary for a point-particle description of the electron to be valid [14,16]. With the reduced equation (Eq. (51)), the radiation damping mγ ...
x may be estimated as and then iterated back to the original equation (Eq. (50)), This approximate equation of motion is free of runaway solutions. Applying Eq. (53) to the RED harmonic oscillator, Eq. (48), we obtain The order of magnitude of each force term is To evaluate the magnitude of the terms in F damp , a random walk model may be used. For a fixed time t = t 0 , a particular value of E     vac (t 0 ) and x(t 0 ). In a two-dimensional random walk model [17], the root-mean-squared distance D rms is given by where N s is the number of steps taken, and ∆s is a typical stepsize; for D The order of magnitude for F damp is evaluated accordingly, Again, using the sharp resonance condition Γω 0 1 (Eq. (8)), F damp is approximated with its dominating term, (29) and (30)), the value of x 0 can be estimated as x 0 3/π h/2mω 0 , which is consistent with Boyer's calculation for the standard deviation of position in Eq. (9). and the equation of motion, Eq. (54), is further simplified As an additional note, given the estimation of E 0 and x 0 in Eq. (57), the three force terms in Eq. (60) have the following relation, The implication is that the spring force F spring dominates, while the RED vacuum field F vac and the radiation damping F damp are perturbations. The driving force F vac and the damping force F damp will establish a balance, keeping the amplitude of the motion of the RED harmonic oscillator in the vicinity of x 0 . Finally, after the main equations for simulation are completely defined (Eqs. (47) and (60)), the range of the integration time τ int (i.e. how long the simulation runs) needs to be chosen. Upon careful inspection, two important time scales can be identified from the analytical solution of Eq. (60), where A is a coefficient determined by the initial conditions, and The first term in Eq. (62) is the transient part of the solution. The characteristic time for the transient motion is The position should be evaluated for τ trans τ int , if we are only interested in the steady-state motion. The second term in Eq. (62) is the steady-state part of the solution. Because the steady-state solution is a finite discrete sum of periodic functions, it would have a nonphysical repetition time τ rep . The choice of τ int should satisfy τ int ≤ τ rep to avoid repetitive solutions. A detailed discussion about τ rep can be found in Appendix C.
Here we give the choice of τ int , where ∆ω is the gap between adjacent frequencies. The frequency gap ∆ω can be estimated using Eqs. (23), (33), and (34), To summarize, Eq. (60) is the equation of motion that will be used for simulating the RED harmonic oscillator. The RED vacuum field in Eq. (60) is simulated by Eq. (47), and the specifications of its modes (k, λ), polarizations ε k,λ , and other relevant variables can be found in Sec. III A. The simplification from Eq. (38) to Eq. (60) relies on only two conditions -namely, the dipole approximation Eq. (4) and the sharp resonance condition Eq. (8), which are also the only two conditions used in Boyer's analysis [8]. The parameters e, m, and ω 0 should be chosen to satisfy these two conditions. Lastly, the integration time τ int of the simulation is to be chosen within the range τ trans τ int ≤ τ rep , where τ tran and τ rep are given in Eqs. (66) and (67) respectively.

IV. RESULTS
In Sec. II B, it was shown that the RED harmonic oscillator probability distribution is Gaussian, contrasting to the double-peak probability distribution of the classical harmonic oscillator. We set out to investigate how the dynamics of the RED harmonic oscillator can give rise to such a Gaussian distribution. In this section, the result of the simulation for an RED harmonic oscillator is presented, and the relation between the trajectory and the probability distribution is discussed.
In constructing the probability distribution from the particle trajectory, two sampling methods are used, namely sequential sampling and ensemble sampling. In a sequential sampling, the position or velocity is recorded in a time sequence from one trajectory, while in an ensemble sampling, the same is recored only at the end of the simulation from an ensemble of trajectories. The recorded positions or velocities are collected in histogram and then converted to a probability distribution, which is compared to the analytical result, Eq. (17). Whereas the sequential sampling illustrates the relationship between the buildup of the probability distribution and the dynamics of the trajectory, the ensemble sampling is convenient for statistical interpretation. In addition, the ensemble sampling is suitable for parallel computing, which improves the computation efficiency. is compared to the trajectory of the RED harmonic oscillator (black solid line). Bottom: A magnified section of the trajectory shows that there is not a fixed phase or amplitude relation between the particle trajectory (black solid line) and the instantaneous driving field (red solid line). The modulation time for the field (red solid line) is also shown to be longer than that for the motion of the harmonic oscillator (black solid line).

A. Sequential Sampling and the Relation to its Probability Distribution
By solving Eq. (60) numerically, the steady-state trajectory for the RED harmonic oscillator is obtained as shown in Fig. 3. For a comparison, the RED vacuum field is also included. The differential equation (Eq. (60)) is solved using the adaptive 5th order Cash-Karp Runge-Kutta method [18], and the integration stepsize is set as one twentieth of the natural period, 1 20 2π ω 0 . The charge e, mass m, and natural frequency ω 0 are chosen as q e , 10 −4 m e , and 10 16 (rad/s) respectively, where q e is the electron charge and m e is the electron mass. The frequency range ∆ of the RED vacuum field is chosen to be much broader than the resonance width, 220 × Γω 2 0 . The choice of 10 −4 m e is made to bring the modulation time and the natural period of the harmonic oscillator closer to each other. In other words, the differential equation (Eq. (60)) covers time scales at two extremes, and the choice of mass 10 −4 m e brings these scales closer to each other while keeping the integration time manageable without losing the physical characteristics of the problem. Some features of the trajectory, shown in Fig. 3, are worth mentioning. First, there is not a fixed phase or amplitude relation between the particle trajectory and the instantaneous driving field. Second, the rate of amplitude modulation in the particle trajectory is slower than that in the driving field. While these results may seem a little surprising at first glance, the steady-state solution in Eq. (62) provides some insights into both of these behaviors.
To understand the features of the trajectory, we study the temporal response of the harmonic oscillator as revealed in the Green formulation of the steady-state solution [19], (69) The solution indicates that the driving field E (x) vac (t ) at any given time t has an effect on the particle for a finite . In other words, the particle trajectory x(t) at a particular time t is affected by the driving field E (x) vac (t ) from all previous moments t ≤ t. This explains why the particle trajectory does not seem to have a fixed phase and amplitude relation with the instantaneous driving field. Another implication of Eq. (69) is that it always takes a finite period of time, 1 Γω 2 0 , for the particle to dissipate the energy gained from the driving field. Thus, even if the field already changes its amplitude, it would take a while before the particle can respond. This is why the amplitude modulation in the particle trajectory is slower compared to that in the driving field 9 . Using a sequential sampling method, a probability distribution is constructed from the simulated trajectory of a RED harmonic oscillator as shown in Fig. 4. The RED probability distribution turns out to be a Gaussian distribution, which is the same as the ground state probability distribution of a quantum harmonic oscillator. To understand the relation between the trajectory and the Gaussian probability distribution, we investigate the dynamics of the RED harmonic oscillator at two time scales. First, at small time scale, the particle oscillates as a classical harmonic oscillator. The oscillation amplitude is constant, and the period is T = 2π/ω 0 . Such an oscillation makes a classical double-peak probability distribution. Second, at large time scale, the oscillation amplitude modulates. As a result, different parts of the trajectory have double-peak probability distributions associated with different oscillation amplitudes, which add to make the final probability distribution a Gaussian distribution as shown in Fig. 5. To further understand how amplitude modulation turns the final probability distribution into a Gaussian distribution, we again inspect the steady-state solution in Eq. (62), The solution can be rewritten in terms of the complex positionx(t), x(t) = Re (x(t)) , Since the frequencies are sampled symmetrically around ω 0 , the amplitude modulation is made evident as the complex positionx(t) is factorized into a modulation term A(t) and an oscillation term e −iω0t , The amplitude A(t) is a sum of many complex components, which rotate in the complex plane at different rates, ω − ω 0 . At any given time, the amplitude A(t) is determined by the configuration of these components, as shown in Fig. 6. As time elapses, the configuration of the components evolves and the amplitude A(t) may change dramatically. However, when the elapsed time ∆t is much smaller than the coherence time τ coh , which is given by the shortest rotating period, the change of amplitude A(t) is negligible, Therefore, since an oscillation amplitude A(t) is approximately constant within one coherence time τ coh 10 , the 10 The coherence time τ coh defined in this work has the same meaning as the temporal width of the first-order correlation function. We can take the autocorrelation of the simulated trajectory and show that its temporal width is the same as the coherence time computed in this work.
amplitude modulation as previously discussed can be characterized by the coherence time τ coh . In other words, the coherence time τ coh is the modulation time.
Now, if we take a representative sampling of the oscillation amplitudes (i.e. each sampled amplitude is several coherence times τ coh apart), we obtain an amplitude distribution as shown in Fig. 7. From the amplitude distribution, it is clear that the occurrence of very large or very small amplitudes is rare. This is because it requires all the components of the amplitude A(t) to be either completely aligned or completely misaligned to obtain the extreme values of A(t). For most of the time, the configuration is in partial alignment, and the amplitude A(t) is of the average size as calculated previously (i.e. x 0 ). Using the amplitude distribution in Fig. 7, a Gaussian distribution can be reconstructed by adding up all the classical double-peak distributions associated with these sampled oscillation amplitudes. This reconstructed Gaussian probability distribution resembles the probability distribution obtained directly from the simulated data as shown in Fig. 7. It is worth pointing out that it is impossible to reconstruct such a Gaussian distribution with any amplitude distribution other than the one given by the particle trajectory (Fig. 7). This result shows that the Gaussian probability distribution of the RED harmonic oscillator is a consequence of its intrinsic amplitude modulation.

B. Ensemble Sampling and Heisenberg Uncertainty Relation
In many RED analyses [2, 4-6, 8, 20], the procedure of random phase averaging is an essential tool in obtaining the statistical properties of a physical system. To compare the RED numerical results with the analyses, it is advantageous to use ensemble sampling. In ensemble sampling, particles in each trajectory of the ensemble are prepared with identical initial conditions, but the RED vacuum fields differ in the initial random phasesθ k,λ . Here, different random initial phasesθ k,λ represent the physical realization of random phase averaging. At the end of the simulation, physical quantities such as position and momentum are recorded from an ensemble of trajectories as shown in Fig. 8. In the case of RED harmonic oscillator, the recorded positions and momentum are converted into probability distributions, and the Heisenberg minimum uncertainty relation is found to be satisfied. Boyer suggested that the minimum uncertainty relation is established by a delicate balance between the gain of energy from RED vacuum field driving and the loss of energy through radiation damping [2]. As we turn off the radiation damping in the simulation (Fig. 9), the minimum uncertainty relation no long holds, although the range of the particle's motion is still constrained by the harmonic potential. This result supports the balancing mechanism proposed by Boyer. Unlike sequential sampling, ensemble sampling has the advantage that the recorded data are fully uncorrelated. As a result, the integration time does not need to be much longer than one coherence time τ coh after the particle reaches its steady-state. However, this does not mean that the computation time will be shorter. In fact, because only one data point is recorded from each trajectory, a simulation with ensemble sampling usually takes much longer than with sequential sampling. For example, a typical RED simulation with sequential sampling (i.e. number of sampled frequencies N ω = 20k) takes 2.3 hours to finish, but a typical RED simulation with ensemble sampling (i.e. number of particles N p = 200k, number of sampled frequencies N ω = 0.5k) takes 61 hours, which is 27 times longer.
A resolution to this problem is to incorporate parallel computing in the ensemble sampling approach. The parallelization of the RED simulation is straightforward for ensemble sampling, since each trajectory are independent except for the random initial phasesθ k,λ . To reduce the amount of interprocessor communication and computation overhead, each processor is assigned an equal amount of work. A typical simulation with ensemble sampling approach (i.e. number of particles N p = 200k, number of sampled frequencies N ω = 0.5k) takes 61 hours to run on a single processor. After parallelization, it takes only 4.5 hours to run on 80 processors, which is a factor of 14 speed-up 11 , and the computation time scales inversely with the number of processors as shown in Fig. 10. In addition, the parallelized code is advantageous for testing the numerical convergence of the simulation. In Fig. 11, the convergence of the ensemble-averaged energy as a function of sampled frequency number is shown. 11 The parallelization of the RED simulation program is developed and benchmarked with assistance from University of Nebraksa Holland Computing Center. The program is written in Fortran and parallelized using Message Passing Interface (MPI) [18,21]. The compiler used in this work is from the GNU Compiler Collection (GCC) "openmpi-1.4.3/gcc-4.4.1". However, with the commercial compiler "openmpi-1.4.4/intel-11.1", we have found that the computation time can be further reduced by a factor of 6.

V. SUMMARY & CONCLUSIONS
The analytical probability distribution of an RED harmonic oscillator is obtained in Sec. II B. The numerical method for the RED simulation is documented in Sec. III. Agreement is found between the simulation and the analytical results. Namely, the probability distributions constructed from the simulated data follows the same Gaussian distributions as the ground state quantum harmonic oscillator (Figs. 4, and 8). As a consequence, the Heisenberg minimum uncertainty relation is satisfied by the widths of the RED position and momentum distributions. In the absence of radiation damping, this minimum uncertain relation no longer holds (Fig. 9), which supports Boyer's energy argument that there is a delicate balance between the vacuum field driving and the radiation damping. Overall, the RED harmonic oscillator is distinct from the classical harmonic oscillator in that its oscillation amplitude modulates at the time scale of coherence time τ coh (Fig. 12), and the resulting RED probability distribution can be thought of as a sum of classical double-peak distributions with a specific amplitude distribution (Fig. 7).

VI. DISCUSSION-APPLICATION OF RED SIMULATION TO OTHER PHYSICAL SYSTEMS
In quantum mechanics, the harmonic oscillator has excited, coherent, and squeezed states. A natural extension of our current work is to search for the RED correspondence of such states. Currently, we are investigating how a Gaussian pulse with different harmonics of ω 0 will affect the RED harmonic oscillator. Can the RED harmonic oscillator support a discrete absorption spectrum, and if so, how does it compare with the prediction from quantum mechanics? Such a study is interesting in the broader view of Milonni's comment that RED is unable to account for the discrete energy levels of interacting atoms [1], and Boyer's comment that at present the line spectra of atoms is still unexplained in RED [22].
The current method of RED simulation may also be applicable to study other quantum systems that are related to the harmonic oscillator, such as a charged particle in a uniform magnetic field and the anharmonic oscillator [4,20]. For the first example, classically, a particle in a uniform magnetic field performs cyclotron motion. Such a system can be viewed as a two-dimensional oscillator, having the natural frequency set by the Larmor frequency. On the other hand, a quantum mechanical calculation for the same system reveals Landau quantization. The quantum orbitals of cyclotron motion are discrete and degenerate. Such a system presents a challenge to RED. For the second example, a harmonic potential can be modified to include anharmonic terms of various strength. Heisenberg considered such a system a critical test in the early development of quantum mechanics [23,24]. We think that a study of the anharmonic oscillator is thus a natural extension of our current study and may serve as a test for RED.
Lastly, over the last decades there has been a sustained interest to explain the origin of electron spin and the mechanism behind the electron doubleslit diffraction with RED [25][26][27][28]. Several attempts were made to con-struct a dynamical model that accounts for electron spin. In 1982, de la Peña calculated the phase averaged mechanical angular momentum of a three-dimensional harmonic oscillator. The result deviates from the electron spin magnitude by a factor of 2 [25]. One year later, Sachidanandam derived the intrinsic spin one-half of a free electron in a uniform magnetic field [26]. Whereas Sachidanandam's calculation is based on the phase averaged canonical angular momentum, his result is consistent with Boyer's earlier work where Landau diamagnetism is derived via the phase averaged mechanical angular momentum of a free electron in a uniform magnetic field [4]. Although these results are intriguing, the most important aspect of spin, the spin quantization, has not been shown. If passed through a Stern-Gerlach magnet, will the electrons in the RED description split into two groups of trajectories 12 ? At this point, the dynamics becomes delicate and rather complex. To further investigate such a model of spin, an RED simulation may be helpful.
On the other hand, over the years claims have been made that RED may predict doubleslit electron diffraction [2,27,28]. In order to explain the experimentally observed electron doubleslit diffraction 13 [35][36][37], different FIG. 13. Schematic illustration of RED doubleslit diffraction. Several authors have proposed a RED (or SED) mechanism that explains the electron doubleslit diffraction [2,27,28]. A central idea provided by RED is that the modes of a vacuum field in one slit are determined by the presence of both slits. This idea allows the superposition principle to be merged with the concept of particle trajectories. mechanisms motivated by RED were proposed [27,28], while a detailed calculation of the diffracting vacuum field in the vicinity of a slit structure has been given [38]. In 1999, Kracklauer suggested that particles steered by the modulating waves of the RED vacuum field should display a diffraction pattern when passing through a slit, since the RED vacuum field itself is diffracted [27]. In recent years, another diffraction mechanism is proposed by Cavalleri et al. in relation to a postulated electron spin motion [28]. Despite of these efforts, Boyer points out in a recent review article that at present there is still not a concrete RED calculation on the doubleslit diffraction [22]. Boyer suggests that as the correlation function of the RED vacuum field near the slits is modified by the slit boundary, the motion of the electron near the slits should be influenced as well. Can the scattering of the RED vacuum field be the central mechanism behind the electron doubleslit diffraction (Fig. 13)? As Heisenberg's uncertainty relation is a central feature in all matter diffraction phenomena, any proposed mechanism for electron doubleslit diffraction must be able to account for Heisenberg's uncertainty relation. This is exactly what RED does for the harmonic oscillator. Thus, we hope that the current simulation method may help providing a detailed investigation for any proposed RED diffraction mechanism. In "unbounded" space, the modes are continuous and the field is expressed in terms of an integral. In "bounded" space, the modes are discrete and the field is expressed in terms of a summation. In both cases, the expression for the field amplitude needs to be obtained (see A 1 and A 2). The integral expression helps comparison with analytical calculations in previous papers [2,[4][5][6]8], while the summation expression is what we use in our numerical work.

Unbounded Space
The homogeneous solution of Maxwell's equations in unbounded space is equivalent to the solution for a wave equation, whereÃ(k, λ) is the undetermined field amplitude for the mode (k, λ) and has the unit of electric field (V/m),k is defined as the unit vector of k, and the two vectors, ε(k, 1) and ε(k, 2), describe an orthonormal polarization basis in a plane that is perpendicular to the wave vector k.

Bounded Space
The solution of homogeneous Maxwell's equations in bounded space has the summation form, E(r, t) = 1 2 k,λ Ã k,λ e i(k·r−ωt) +Ã * k,λ e −i(k·r−ωt) ε k,λ , B(r, t) = 1 2c k,λ Ã k,λ e i(k·r−ωt) +Ã * k,λ e −i(k·r−ωt) ξ k,λ , where ξ k,λ =k × ε k,λ ,Ã k,λ is the undetermined field amplitude for the mode (k, λ) and has the unit of electric field (V/m),k is defined as the unit vector of k, and the two vectors, ε k,1 and ε k,2 , describe an orthonormal polarization basis in a plane that is perpendicular to the wave vector k. Using we can follow the same argument in A 1 and obtain the phase averaged energy density in bounded space whereÃ k,λ = A k,λ e iθ k,λ . Again, if we postulate that the vacuum energy ishω/2 for each mode (k, λ), then in a bounded space of volume V the vacuum energy density is (A18) 14 If the two modes are not identical (i.e. k = k or λ = λ), then e iθ k,λ and e iθ k ,λ are independent, which leads to the factorization e i(θ k,λ ±θ k ,λ ) θ = e iθ k,λ θ e ±iθ k ,λ θ = 0.

Appendix B: Isotropic Polarization Vectors
A wave vector chosen along the z-axis, has an orthonormal polarization basis in the xy-plane, In the same manner, we can rotateε k,1 andε k,2 with the rotation matrixR, and obtain an isotropically dis-  (B5)

Appendix C: Repetitive Time
A field composed of finite discrete frequencies, repeats itself at the least common multiple (LCM) of all the periods of its frequency components, where T k = 2π ω k . An example of a two-frequency beat wave is given in Fig. 14. 15 After the rotation, the uniformly distributed circle will span into a uniformly distributed spherical surface.
Given the frequency spectrum of E(t), one can draw a relation between the repetition time τ rep and the greatest common divider (GCD) of the frequencies, τ rep = 2π (ω 1 , ω 2 , . . . , ω N ) GCD . (C4) The derivation of the relation in Eq. (C4) is the following. First, we factorize the sum of all the frequencies into two terms, where ∆ω is the smallest frequency gap (∆ω ≤ ∆ω gcd ), suffices our purpose. The frequency gap as a function of κ can be estimated using Eq. (33), Applying the sharp resonance condition (Eq. (8)) to Eqs. (23) and (34), it can be further shown that ∆κ is much smaller than κ 0 and κ κ 0 ,