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

Abstract

Stochastic electrodynamics (SED) predicts a Gaussian probability distribution for a classical harmonic oscillator in the vacuum field. This probability distribution is identical to that of the ground state quantum harmonic oscillator. Thus, the Heisenberg minimum uncertainty relation is recovered in SED. To understand the dynamics that give rise to the uncertainty relation and the Gaussian probability distribution, we perform a numerical simulation and follow the motion of the oscillator. The dynamical information obtained through the simulation provides insight to the connection between the classic double-peak probability distribution and the Gaussian probability distribution. A main objective for SED 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 SED.

1. Introduction

According to quantum electrodynamics, the vacuum is not a tranquil place. A background electromagnetic field, called the electromagnetic vacuum field, is always present, independent of any external electromagnetic source. The first experimental evidence of the vacuum field dates back to 1947 when Lamb and his student Retherford found an unexpected shift in the hydrogen fine structure spectrum. The physical existence of the vacuum field has inspired an interesting modification to the classical mechanics, known as stochastic electrodynamics (SED). As a variation of classical electrodynamics, SED adds a background electromagnetic vacuum field to the classical mechanics. The vacuum field as formulated in SED has no adjustable parameters except that each field mode has a random initial phase and the field strength is set by the Planck constant. With the aid of this background field, SED is able to reproduce a number of results that were originally thought to be pure quantum effects.

Despite that the classical mechanics and SED are both theories that give trajectories of particles, the probability distributions of the harmonic oscillator in both theories are very different. In a study of the harmonic oscillator, Boyer showed that the moments of the position of an SED harmonic oscillator are identical to those of the ground state quantum harmonic oscillator. As a consequence, the Heisenberg minimum uncertainty relation is satisfied, and the probability distributions of an SED harmonic oscillator are a Gaussian, identical to that of the ground state quantum harmonic oscillator. While in classical mechanics it is most likely to find an oscillator at the two turning points of the trajectory, hence the double-peak probability distribution, the SED Gaussian probability distribution has the maximum in the center. How do the dynamics differ in the two classical theories so that the probability distributions become so different?

Although the analytical solution of an SED harmonic oscillator was given a long time ago, it is not straightforward to see the dynamical properties of a single particle from the complicated solution. Additionally, many results in SED such as probability distribution are obtained from ensemble phase averaging. Thus, they cannot be used for the interpretation of a single particle’s dynamical behavior in time. Some may consider applying the central limit theorem and treat the positions in the time sequence as independent random variables; however, for an SED system this cannot work, because the correlation between the motion at two points in time persists beyond many cycles of oscillation, unless the system is proved to be ergodic. As an analytical proof of the ergodicity is very difficult, we take a numerical approach to study the particle’s dynamical behavior. Besides what is already known from the analytical solution, our numerical studies construct the probability distribution from a single particle’s trajectory. We investigate the relation between such a probability distribution and the particle’s dynamical behavior. Ultimately, we want to know the underlying mechanism that turns the classic double-peak distribution into the Gaussian distribution.

While most works in the field of SED are analytical, numerical studies are rare. The advantage of numerical simulation is that it may be extended to other physical systems with relative ease and is flexible in testing different assumptions and approximations. For example, in most SED analyses the effect of the Lorentz force due to the vacuum field is often neglected as the first-order approximation. However, when the field gradient is nonzero, the Lorentz force from the magnetic field can work with the electric field to give a nonzero drift to a charged particle, known as the ponderomotive force. Therefore, the Lorentz force of the vacuum field may play a significant role in SED. In fact, it is shown in the literature that the magnetic part of the vacuum field is responsible for the self-ionization of the atoms and the acquisition of an energy of ten to the twentieth electronvolts in a few nanoseconds on the part of a free electron in vacuum. Using numerical simulation, one can easily examine SED beyond the first-order approximation — in this work we limit ourselves to physical parameters where we expect the Lorentz force effects to be minimal, as this work is compared to analytical results where this approximation is made — and the Lorentz force effects in SED can be investigated.

The major challenge for the numerical simulation is to properly account for the vacuum field modes. A representative sampling of the modes is thus the key for successful simulations. In this study, one of our goals is to use a simple physical system, namely, the simple harmonic oscillator, to benchmark our numerical method of vacuum mode selection so that it can be used to test the validity range of SED as discussed in the following.

Over the decades, SED has been criticized for several drawbacks. Authors like Cavalleri argued that SED can neither explain electron slit-diffraction nor derive the nonlinear Schrödinger equation; moreover, SED implies broad radiation and absorption spectra for rarefied gases. In the case of the quartic anharmonic oscillator, Pesquera and Claverie showed that SED disagrees with quantum mechanics. Additionally, some results claimed in the literature were shown to be wrong due to improper relativistic approximation. While these theoretical analyses are documented in detail, it may be useful to use numerical simulation as an independent check in order to establish the validity range of SED, which is one of the objectives of our work.

Meanwhile, a modified theory called stochastic electrodynamics with spin (SEDS) was recently proposed. As SEDS is an extension of SED with a model of electron spin motion, it is argued that the introduction of electron spin can eliminate several drawbacks of SED. Among all, it is claimed that SEDS allows the derivation of the complete and even the generalized Schrödinger equation. Also, it is claimed that SEDS can explain electron slit-diffraction and the sharp spectral lines of the rarefied gases. In view of the fact that SEDS is a modified theory from SED and may extend its validity range, it is probably important as a next step to apply numerical simulation to such a model as a "quasi-experiment" and test some of its claims. Models that include constraints can be incorporated in our numerical model with, for example, Lagrange multipliers.

The organization of this paper is the following. First, in Section 2 Boyer’s results on the SED harmonic oscillator are briefly reviewed. Based on these results, the probability distribution for the oscillator is derived. Second, in Section 3 details for simulating the vacuum field and the SED harmonic oscillator are documented. Third, in Section 4 the trajectory of the SED 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 method is "sequential sampling," which is suitable for studying the relation between the dynamics of the SED harmonic oscillator and its probability distribution. The second approach is "ensemble sampling," which lends itself well to parallel computing and is convenient for statistical interpretations. Lastly, in Section 6 we discuss some potential applications of the numerical simulation in studying other quantum phenomena. 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 also assist in assessing SED’s validity range.

2. Theory of Stochastic Electrodynamics

2.1. Brief Review of Boyer’s Work

In his 1975 papers, Boyer calculated the statistical features of an SED harmonic oscillator, and the Heisenberg minimum uncertainty relation is shown to be satisfied for such an oscillator. The vacuum field used in Boyer’s work arises from the homogeneous solution of Maxwell’s equations, which is assumed to be zero in classical electrodynamics. In an unbounded, free space, the vacuum field has an integral form: a sum over two polarizations of an integral over all wave vectors, of the polarization vector times a field amplitude times a phase factor, where the angular frequency is the speed of light times the magnitude of the wave vector and the phase is a random variable uniformly distributed between zero and two pi. The two polarization unit vectors are perpendicular to the wave vector and mutually orthogonal.

To investigate the dynamics of the SED harmonic oscillator, Boyer used the dipole approximation — the product of the wave vector and the position is much smaller than one — to remove the spatial dependence in the vacuum field. Therefore, the equation of motion for an SED harmonic oscillator used in Boyer’s analysis is: mass times acceleration equals minus mass times the square of the natural frequency times the displacement, plus mass times the radiation damping parameter times the third time derivative of the displacement, plus charge times the component of the vacuum field along the axis of motion. The radiation damping parameter here is two thirds of the squared charge divided by the mass and the cube of the speed of light, taken with the Coulomb constant.

Boyer further calculated the standard deviation of position and momentum from the steady-state solution by averaging over the random phase. The standard deviation of position is the square root of the reduced Planck constant divided by twice the mass and the natural frequency, and the standard deviation of momentum is the square root of the reduced Planck constant times the mass times the natural frequency, divided by two. Here the phase averaging represents the ensemble average over many realizations; in each realization, the random phase of the vacuum field is different. The above result satisfies the Heisenberg minimum uncertainty relation: the product of the two standard deviations equals one half of the reduced Planck constant. From an energy argument, Boyer showed that this uncertainty relation can also be derived from a delicate balance between the energy gain from the vacuum field and the energy loss through radiation damping.

2.2. Probability Distribution

Given the knowledge of the moments of the position, the Fourier coefficients of the probability distribution can be determined by Taylor expanding the Fourier kernel in powers of the position, term by term. Using the steady-state solution together with Boyer’s relations for the averages of phase factors — the average of a factor carrying the sum of two random phases vanishes, while the average of a factor carrying their difference is non-zero only when the two modes coincide — the moments can be evaluated, and the resulting probability distribution is a Gaussian whose width is the standard deviation given above. This is identical to the position probability distribution of the ground state quantum harmonic oscillator.

3. Methods of Numerical Simulation

(Omitted for length: Section 3 documents the numerical method in full — the representation of the vacuum field in bounded space as a discrete sum over modes with a field strength fixed by the Planck constant, the scheme for selecting a representative sample of vacuum modes and the frequency gap and bandwidth that scheme implies, the random-walk estimate of the oscillation amplitude, and the reduction of the third-derivative equation of motion to a form suitable for numerical integration. The complete text is at the source.)

4. Simulation Results

In Section 2, it was shown that the probability distribution for an SED harmonic oscillator is a Gaussian. In Section 3, we develop the methods for a numerical simulation to investigate the dynamics of the SED harmonic oscillator and how it gives rise to the Gaussian probability distribution. In this section, the results of the simulation are presented, and the relation between the trajectory and the probability distribution is discussed.

To construct the probability distribution from a particle’s trajectory, two sampling methods are used. The first method is sequential sampling and the second method is ensemble sampling. In sequential sampling the position or velocity is recorded in a time sequence from a single particle’s trajectory, while in ensemble sampling the same is recorded only at the end of the simulation from an ensemble of particle trajectories. The recorded positions or velocities are collected in histogram and then converted to a probability distribution for comparison to the analytical result. Whereas the sequential sampling illustrates the relationship between the buildup of probability distribution and the dynamics of particle trajectory, the ensemble sampling is convenient for statistical interpretation. In addition, the ensemble sampling is suitable for parallel computing, which can be used to improve the computation efficiency.

4.1. Particle Trajectory and the Probability Distribution

By solving the equation of motion numerically, the steady-state trajectory for the SED harmonic oscillator is obtained. For a comparison, the temporal evolution of the vacuum field is also included. The ordinary differential equation is solved using the adaptive fifth order Cash-Karp Runge-Kutta method, and the integration stepsize is set as small as one-twentieth of the natural period to avoid numerical aliasing. The charge is the electron charge; the mass is ten to the minus fourth times the electron mass; the natural frequency is ten to the sixteenth radians per second; and the vacuum field frequency range is 220 times the resonance width of the harmonic oscillator. The choice of mass is made to bring the modulation time and the natural period of the harmonic oscillator closer to each other. In other words, the equation of motion covers time scales at two extremes, and this choice of mass brings these two scales closer so that the integration time is manageable without losing the physical characteristics of the problem.

Here we would like to highlight some interesting features of the simulated trajectory. First, there appears to be no 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. To gain insights into these dynamical behaviors, we study the steady-state solution in the Green function form: the displacement at a given time is the charge divided by the mass and the shifted resonance frequency, times an integral over all earlier times of the vacuum field, damped by an exponential in the elapsed time at half the radiation-damping rate and modulated by the sine of the shifted resonance frequency times the elapsed time. The solution indicates that the effect of the driving field at any given time lasts for a time period equal to the inverse of the radiation-damping rate beyond that time. In other words, the particle motion at a given time is affected by the vacuum field from all the previous moments. As the vacuum field fluctuates in time, the fields at two points in time only become uncorrelated when the time separation is much longer than one coherence time. This property of the vacuum field reflects on the particle trajectory, and it explains why the particle trajectory has no fixed phase or amplitude relation with the instantaneous driving field. Another implication is that it takes that same characteristic time for the particle to dissipate the energy gained from the instantaneous driving field. Thus, even if the field already changes its amplitude, it would still take a while for the particle to follow. This explains why the amplitude modulation in the particle trajectory is slower compared to that in the driving field. In case of slow field modulation, when the field bandwidth is shorter than the resonance width of the harmonic oscillator, the modulation time of the field and the particle trajectory are the same.

The sequential sampling of a simulated trajectory gives the probability distributions of position and momentum. While Boyer’s result is obtained through ensemble phase averaging, the Gaussian probability distribution shown here is constructed from a single trajectory and is identical to the probability distribution of a ground state quantum harmonic oscillator.

To understand how the trajectory gives rise to a Gaussian probability distribution, we investigate the particle dynamics at two time scales. At short time scale, the particle oscillates in a harmonic motion. The oscillation amplitude is constant, and the period is two pi divided by the natural frequency. Such an oscillation makes a classical double-peak probability distribution. 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. To verify this idea, we attempt to reconstruct the Gaussian probability distribution from the double-peak probability distributions at different sections of the trajectory. We approach this problem by numerically sampling the oscillation amplitudes at a fixed time-step.

To determine the appropriate sampling time-step, we inspect the steady-state solution in its complex form and factorize it into an amplitude term and an oscillation term at the natural frequency. The amplitude term is a sum of complex components, each rotating in the complex plane at a rate given by the detuning of its mode from the natural frequency. At any given time, the configuration of these components determines the magnitude of the amplitude. As time elapses, the configuration evolves and the amplitude changes with time. When the elapsed time is much shorter than the shortest rotating period — two pi divided by the largest detuning — the change in the amplitude is negligible. We denote this shortest rotating period as the coherence time. The coherence time as defined here is equivalent to the temporal width of the first-order correlation function, the autocorrelation: as the autocorrelation of the simulated trajectory is the Fourier transform of the spectrum according to the Wiener-Khinchin theorem, it has a temporal width the same as the coherence time calculated here. For our problem at hand, it is clear that the sampling of oscillation amplitudes should use a time-step greater than the coherence time.

A representative sampling of the oscillation amplitudes, with each sampled amplitude separated by three coherence times, shows in histogram that the occurrence of large or small amplitudes is rare. Most of the sampled amplitudes have a medium value. This is because the occurrence of extreme values requires complete alignment or misalignment of the complex components in the amplitude. For most of the time, the complex components are in partial alignment and thus give a medium value of the amplitude. Interestingly, the averaged value of the amplitude is close to the oscillation amplitude predicted by the random walk model. Using the amplitude distribution obtained this way, a probability distribution can be constructed by adding up the double-peak probability distributions: the total distribution is the integral over amplitude of the double-peak distribution for that amplitude, weighted by the amplitude distribution. This constructed probability distribution is a Gaussian and is identical to the simulation result. The reconstruction of the Gaussian probability distribution indicates the transitioning from double-peak distribution to the Gaussian distribution due to the amplitude modulation driven by the vacuum field.

4.2. Phase Averaging and Ensemble Sampling

In many SED analyses, the procedure of random phase averaging is often used to obtain the statistical properties of the physical system. A proper comparison between numerical simulation and analysis should thus be based on ensemble sampling. In each realization of ensemble sampling, the particle is prepared with identical initial conditions, but the vacuum field differs in its initial random phase. The difference in the initial random phase corresponds to the different physical realizations in random phase averaging. At the end of the simulation, physical quantities such as position and momentum are recorded from an ensemble of trajectories.

The ensemble sampling of the simulation gives position and momentum probability distributions that satisfy the Heisenberg minimum uncertainty as predicted by Boyer’s analysis. In addition, Boyer proposed a mechanism for the minimum uncertainty using an energy-balance argument. Namely, he calculated the energy gain from the vacuum field and the energy loss through radiation damping, and he found that the delicate balance results in the minimum uncertainty relation. We confirm this balancing mechanism by turning off the radiation damping in the simulation and see that the minimum uncertainty relation no longer holds, although the range of the particle’s motion is still bounded by the harmonic potential.

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 very long compared to the coherence time. However, since only one data point is recorded from each trajectory, a simulation with ensemble sampling actually takes longer time than with sequential sampling. For example, a typical simulation run with sequential sampling takes 2.3 hours to finish for twenty thousand sampled frequencies, but with ensemble sampling it takes 61 hours for two hundred thousand particles and five hundred sampled frequencies. A remedy to this problem is to use parallel computing for the simulation. The parallelization scheme for our simulation with ensemble sampling is straightforward, since each trajectory is independent except for the random initial phases. To reduce the amount of interprocessor communication and computation overhead, each processor is assigned an equal amount of work. The parallelized program is benchmarked and shows an inverse relation between the computation time and the number of processors. As the computation speedup is defined as the single-processor computation time divided by the multiprocessor computation time, the inverse relation indicates ideal performance of linear speedup. As an additional note, the parallelized code is advantageous for testing the numerical convergence of the simulation.

5. Conclusions

The analytical probability distribution of an SED harmonic oscillator is obtained in Section 2. The details of our numerical methods including vacuum mode selection are documented in Section 3. Agreement is found between the simulation and the analytical results, as both sequential sampling and ensemble sampling give the same probability distribution as the analytical result. Numerical convergence is reached with a low number of sampled vacuum field modes, five hundred, which is an indication that our method of vacuum mode selection is effective in achieving a representative sampling.

As the probability distribution constructed from a single trajectory is a Gaussian and satisfies the Heisenberg minimum uncertainty relation, we investigate the relation between the Gaussian probability distribution and the particle’s dynamical properties. As a result, the amplitude modulation of the SED harmonic oscillator at the time scale of the coherence time is found to be the cause for the transitioning from the double-peak probability distribution to the Gaussian probability distribution.

6. Discussions: Application of 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 SED correspondence of such states. Currently, we are investigating how a Gaussian pulse with different harmonics of the natural frequency will affect the SED harmonic oscillator. Can the SED harmonic oscillator support a discrete excitation 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 SED is unable to account for the discrete energy levels of interacting atoms and also Boyer’s comment that at present the line spectra of atoms are still unexplained in SED.

The methods of our numerical simulation may 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. 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 SED. For the second example, a harmonic potential can be modified to include anharmonic terms of various strengths. Heisenberg considered such a system a critical test in the early development of quantum mechanics. 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 SED.

Lastly, over the last decades there has been a sustained interest to explain the origin of electron spin and the mechanism behind the electron double-slit diffraction with SED. Several attempts were made to construct 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. One year later, Sachidanandam derived the intrinsic spin one-half of a free electron in a uniform magnetic field. 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 an electron in a uniform magnetic field. 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 SED description split into two groups of trajectories? The electron Stern-Gerlach effect is an interesting but controversial topic in its own right: whereas Bohr and Pauli asserted that an electron beam cannot be separated by spin based on the concept of classical trajectories, Batelaan and colleagues and Dehmelt argue that one can do so with certain Stern-Gerlach-like devices. At this point, the dynamics become delicate and rather complex. To further investigate such a model of spin, a numerical simulation may be helpful.

On the other hand, over the years claims have been made that SED can predict double-slit electron diffraction. In order to explain the experimentally observed electron double-slit diffraction, different mechanisms motivated by SED were proposed, but no concrete calculation has been given except for a detailed account of the slit-diffracted vacuum field. In 1999, Kracklauer suggested that particles steered by the modulating waves of the SED vacuum field should display a diffraction pattern when passing through a slit, since the vacuum field itself is diffracted. In recent years, another diffraction mechanism is proposed by Cavalleri and colleagues in relation to a postulated electron spin motion. Despite these efforts, Boyer points out in a recent review article that at present there is still no concrete SED calculation on the double-slit diffraction. Boyer suggests that as the correlation function of the 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 vacuum field be the physical mechanism behind the electron double-slit diffraction? As Heisenberg’s uncertainty relation is a central feature in all matter diffraction phenomena, any proposed mechanism for electron double-slit diffraction must be able to account for Heisenberg’s uncertainty relation. In the physical system of the harmonic oscillator, SED demonstrates a mechanism that gives rise to the Heisenberg minimum uncertainty. We hope that the current simulation method may help in providing a detailed investigation on the proposed SED mechanisms for the electron slit-diffraction.

(The appendices on the vacuum field in unbounded and bounded space, the figures and the reference list are omitted for length; the complete text is at the source.)