Probability Calculations Within Stochastic Electrodynamics

Several stochastic situations in stochastic electrodynamics (SED) are analytically calculated from first principles. These situations include probability density functions, as well as correlation functions at multiple points of time and space, for the zero-point (ZP) electromagnetic fields, as well as for ZP plus Planckian (ZPP) electromagnetic fields. More lengthy analytical calculations are indicated, using similar methods, for the simple harmonic electric dipole oscillator bathed in ZP as well as ZPP electromagnetic fields. The method presented here makes an interesting contrast to Feynman’s path integral approach in quantum electrodynamics (QED). The present SED approach directly entails probabilities, while the QED approach involves summing weighted paths for the wave function.


INTRODUCTION
This article is largely intended for providing a new calculational method for deducing the stochastic properties of electromagnetic radiation and classical charged particles in the theory called stochastic electrodynamics (SED). This new calculational method to be described here might well be useful in some contexts outside of SED; nevertheless, the main focus here will indeed be SED.
Summarizing quickly, SED is a classical physical theory involving classical charged particles and classical electromagnetic (E&M) radiation, where Maxwell's classical, microscopic electromagnetic equations hold. The motion of point charges is assumed to be described by the relativistic Lorentz-Dirac classical equation of motion. What is particularly interesting about SED is that the basic assumptions of SED are few, they are not complicated, and their basis makes clear physical sense. A number of physicists over the years, including the author, have felt that SED might not only be a substitute for the part of quantum theory (QT) consisting of quantum mechanics (QM) and quantum electrodynamics (QED), but much more so, provide the basis to derive or deduce QM and QED, or, at the very least, to provide a deeper physical understanding of QM and QED.
To emphasize this point, many well-known systems traditionally analyzed in QM, such as the simple harmonic oscillator (SHO) in either one, two, or three dimensions, fluctuating electric dipole SHOs, including interacting systems of such electric dipoles, plus van der Waals forces, Casimir forces, the thermal-like behavior of electrodynamic systems uniformly accelerated through the "vacuum," diamagnetism, aspects of hydrogen, and blackbody radiation dynamics [1,2], have all been analyzed within the classical theory of SED. This range of "QM" phenomena, that has always been considered outside the domain of classical physics, became understandable in a coherent, consistent, and logical manner in SED, without needing to draw on any extraneous "physical or phenomenological" concepts. Moreover, not only did these classical physics, calculations with ZP or ZPP classical E&M radiation provide close connections with QED results, in some cases, the SED results also preceded QED calculations, such as pioneered by Boyer in the case of uniformly electrodynamic system through the "vacuum" [3][4][5][6][7], or even in the case of Casimir and van der Waals forces, where the physics in SED seems clearer, and changing from T 0 calculations to T > 0 calculations is much simpler in the case of SED than QED [8].
Recently, Boyer made a strong argument [9] that SED is the best classical physical theory for describing physical phenomena, as it contains a range of QM and QED predictions, in addition to the expected classical physics. Moreover, these QT predictions by SED are remarkable, not just in their prediction via classical physics, but also that the supposed faults of classical physics, such as the collapse of Rutherford's orbital model of the atom, or no explanation for van der Waals or Casimir forces, are "repaired" by taking into account the full interaction between classical charges and classical electromagnetic radiation.
Some of the key original articles on SED were in the 1960s, by Marshall [10,11] and Boyer [12,13]. As expected for a classical physical theory, SED recognizes that accelerated and deaccelerated charges create E&M radiation; in turn, E&M radiation effects the motion of charged particles via the Lorentz force term. However, what is significantly different for SED from classical physics as typically taught, is the recognition that if classical charged particles can exist in a thermodynamic equilibrium state, such as for atoms and molecules, then this can only be done in a stochastic equilibrium between charges and E&M fields. The E&M fields fluctuate, as expected in blackbody radiation, but the charges also fluctuate in position, since the two entities are connected together. Fluctuations of one result in fluctuations of the other. In other words, the dilemma that Rutherford immediately recognized after proposing his "miniature solar system" model of the atom, with electrons orbiting the nucleus, is addressed in SED. If only classical charges are present, then Rutherford's model will result in E&M energy radiated off as the electrons orbit, and the orbits will collapse. We now know that the time for decay, for a classical hydrogen atom, is about 1.3 × 10 −11 s, starting from the Bohr radius [14]. Moreover, if the charges were attempted to be held in some static configuration (no orbiting), we know from Earnshaw's theorem that a stable stationary equilibrium configuration is also not possible [15]. Thus, classical E&M radiation and classical charged particles must both be present if there is any hope for equilibrium, with of course radiation effecting the particles, and particles creating radiation, and a stochastic balance resulting.
What the early SED researchers recognized is that to obtain thermodynamic equilibrium between radiation and charges, there must be special stochastic properties of the radiation, and consequently also of the charges. They deduced that these interesting relationships must logically be deductible at all temperatures, indeed, even at T 0, which gave rise to the notion of classical electromagnetic zero-point (ZP) radiation.
Several good reviews exist on all of this work: Ref. [16] provides an excellent history on the development of SED. Other reviews of interest are [9,[17][18][19]. These reviews discuss the deductions made by researchers about the properties of classical electromagnetic ZP radiation, such as Lorentz invariant [12,20], and that the fundamental definition of T 0 must be obeyed by ZP radiation [1,21,22]. Some of the more recent work, such as on hydrogen in SED, is briefly outlined in Ref. [23].
The outline of this article is as follows. In Sections 2.1 and 2.2, certain stochastic properties will be calculated for the E&M radiation fields in SED, first using a new technique that was covered to a lesser extent in Section 3 of Ref. [23]. Some comparisons of this approach will also be made to earlier n-point correlation function approaches by others, particularly by Boyer [24], as well as Marshall [20], the results of which have been used extensively, and extended, by others, including this author. Section 2 has several subsections, including Section 2.3, which provides checks on the results derived in Section 2.2. Section 2.4 relates results in Section 2.2 to the multivariate normal distribution, and why the latter applies to the stochastic fields described here.
Section 3 ends with some concluding remarks. In particular, some brief comments are first made about extending this calculational method to the electric dipole oscillator immersed in these stochastic fields, meaning either ZP or ZPP fields as treated in SED. The electric dipole oscillator, one of the first systems analyzed in SED, is discussed in terms of why this system "easily" lends itself to the stochastic method discussed here. The calculations are long with this method, but can be carried out. In contrast, a system like the classical hydrogen atom is far more difficult and has not yet been shown to be tractable. The point is made that a similar situation existed for the Feynman path integral approach.

PROBABILITY DENSITY FUNCTION CALCULATIONS FOR ELECTRIC AND MAGNETIC RADIATION FIELDS 2.1 Introduction to the Calculational Method
We will begin the calculations in this article by determining the probability density functions for various stochastic properties of the classical electromagnetic ZP radiation fields in SED, as well as the zero-point plus Planckian (ZPP) fields. Actually, both situations can be treated at once, with the temperature T in the expressions allowing the distinction, with 0 ≤ T.
The present Section 2 will address how to determine the probabilistic functions of the E&M stochastic fields in SED. The subsequent section, 3, will then turn to calculating the stochastic properties of a classical electric dipole SHO within ZPP radiation. The main difference between what will be done in this article and what is typically done in SED is that normally the mean, the variance, and correlation quantities of the stochastic E&M fields, and the position and momentum of the oscillator's fluctuating particle are directly calculated. These quantities are sometimes, although not always, calculated at different times and positions in SED. For some problems, the correlations at different times and positions are critical, although this is not the case for many other types of SED problems. One place where clearly the different times and positions were essential in the system analysis, had to do with the uniform acceleration of electrodynamic systems through the vacuum, such as in Refs. [3][4][5][6][7].
Of course, given any probability distribution or density function associated with a single random variable X, one can calculate an infinite number of moments, such as where the angled brackets mean "expectation value," and where P(x) is the probability density function associated with X. However, for a Gaussian probability density distribution of a single random variable, there are only two defining parameters, namely, the mean, or μ ≡ 〈X〉, and the variance, σ 2 ≡ 〈(X − 〈X〉) 2 〉. All other moments 〈(X − 〈X〉) N 〉, for N 1, 2, 3, . . ., can be expressed in terms of just µ and σ.
For a multivariate Gaussian distribution, which actually defines the classical ZP and ZPP E&M radiation, each frequency and polarization component of the radiation is assumed to be governed by an independent Gaussian process. Thus, the situation is more complicated. Nevertheless, the previous paragraph helps to identify the difference between past work in SED involving stochastic process calculations, which largely dealt with moment calculations like 〈(X − 〈X〉) N 〉, whereas here we will directly calculate quantities like P(x).
Early on in SED, the "two-point correlation functions" of the stochastic fields, at different times and positions, were calculated in detail, and were used to deduce "n point" correlation functions, again where different times and spatial positions were included when calculating these correlation functions. An excellent source for investigating deeply these "n point" correlation functions was by Boyer in Ref. [24], where not only were the SED correlation functions determined, but also compared to the expectation values of the corresponding QED functions, involving annihilation and creation operators. The quantities were shown to be in agreement, provided the quantum operators were symmetrized. In contrast, in this article, the probabilities of these quantities will be directly calculated, rather than calculating individual moments of fields and oscillator coordinates. This same approach can also be used in a similar manner for the probability distribution for linear nonrelativistic electric dipole SHOs, although the calculations become even longer than for the fields. Brief comments will be made on this topic in Section 3.
As often done in SED, where the ZP and ZPP fields are critically important to the final physics results, the radiation fields at temperature T ≥ 0 are characterized by T of course, but must also be thought of as having a rapid variation in space and time. This aspect was tackled by Planck using classical physics, covered in the first half of his famous book [25], "The Theory of Heat Radiation," which is still basically represented by the SED theory. The second half of his book, however, introduces quantum concepts involving energy and frequency related to what QT now treats as photons. SED certainly avoids this direction, but the first part of Planck's work, which was also used later by Einstein and Hopf [26,27], still applies to the beginnings of SED. Indeed, the work by Marshall and Boyer in SED has parallels to this early work by Planck and Einstein and Hopf, with the important caveat that equilibrium radiation must exist at T 0 [12,13], and the recognition that ZP radiation is key to getting the stochastic thermodynamic behavior of classical charged particles and classical E&M correct.
To adequately describe the "radiation dynamics" in SED, usually a large region of space is considered, where "large" means compared to the size that any charged particles representing atomic systems are encompassing or traversing. Thus, in the same vein as Planck, Einstein, and Hopf, SED typically considers a rectangular parallelepiped region in space, with dimensions L x , L y , and L z , along the x, y, and z axes. Other shapes can in principle be used, but a rectangular parallelepiped offers mathematical simplicity, without effecting the physical description if the volume is large. The radiation fields representing ZP or ZPP conditions are typically expressed as an infinite sum of plane waves, with periodic boundary conditions (bcs) imposed. The imposition of periodic bcs makes use of the Fourier analysis process for representing the fields, such that if the region is large enough, then the imposition of periodicity does not affect the physical analysis, but does simplify the subsequent mathematical analysis.
Thus, the following expressions for the "free" electric E(x, t) and magnetic B(x, t) radiation fields in this large parallelepiped volume can be written as the following sum of plane waves [16]: where and n x , n y , and n z are integers, and ω n c|k n |, k n · ε kn,λ k n · ε k n ,λ′ 0, and ε kn,λ · ε k n ,λ′ 0 for λ ≠ λ ′ , where λ and λ ′ indicate the linear polarization direction. Specifically, λ might be represented by the values 1 or 2, and the same for λ ′ . Also, k n k n /|k n |. Equations 1 and 2 satisfy the wave equations of ∇ 2 E(x, t) 1 picture is conceptually fairly simple, although of course computationally intensive. Each time a charged particle system, such as an "oscillator or atom," undergoes its motion due to an atomic binding force, and due to the radiation fields in Eqs. 1 and 2, the subsequent scenario of particle motion and radiation fluctuations would represent a real situation of say, the classical atom inside a large cavity kept at temperature T. However, "redoing" the simulation or "experiment" would be carried out with a different set of A kn,λ and B kn,λ coefficients, to represent a similar but different initial set of conditions. Each time a similar "experiment" of radiation and charges is considered, the experiment is treated as another member of the ensemble of similar experiments. This part of SED coincides with the thoughts in the first half of Planck's major treatise [25] and later by Einstein's and Hopf [26,27]. Related SED discussions can be found in Refs. [16,17]. However, most of the following relationships make physical sense without even referring to these references.
Specifically, the expectation value of these coefficients characterizing the ensemble of ZPP radiation fields is of course zero: Moreover, the "A and B coefficients" are considered to be both independent and uncorrelated random variables in this ensemble, so as are the "A coefficients" with different indices, and the same for the "B coefficients": However, for two "A coefficients" with the same indices, and similarly for the "B coefficients," then of course, these quantities cannot be zero, but are assumed to be functions of the frequency of the radiation and of the temperature T: Although Eqs. 4-7 are indeed assumed in SED, the more general relationship that includes all these relationships, plus more, is that the A ′ s and B ′ s coefficients are assumed to be independent random variables, with zero mean as in Eq. 4, with variance [σ(ω n , T)] 2 , and with probability density distributions characterizing the ensemble of possible radiation situations characterizing thermal radiation at temperature T, for 0 ≤ T, as being Gaussian distributions. Specifically: with the same also holding for P(B kn,λ ). From these independent Gaussian distributions, Eqs. 4-7 also follow.
Initial understanding of the importance of the statistical properties of ZP and ZPP in SED, and how these properties relate to the resulting fluctuating and equilibrium properties of charged particles interacting with this radiation, focused a fair bit on [σ(ω n , T)] 2 [16]. The Lorentz invariant property of ZP found independently by Marshall [20] and Boyer [12], is due to the functional form connected to this function. Similarly, the thermodynamic connection of the meaning of T 0 to this radiation and interacting particles is also tied to [σ(ω n , T)] 2 [1,21,22]. Other work by Boyer deduced additional symmetry properties of the required classical E&M nature of ZP and ZPP radiation that involved scaling and conformal invariances [32].
All of these analyses have led to the following in SED for the ZPP spectrum: Note: lim term is the Planckian part. Using Eqs. 9, 1, and 2, the ensemble average of the net energy due to these thermal radiation fields for 0 ≤ T, within the L x × L y × L z rectilinear parallelepiped, can be calculated. Specifically, using the relationships above, and the usual relationship between electromagnetic energy in free space and the E&M fields [15], yields where ω n c|k n | follows from Eq. 3; n is composed of {n x , n y , n z }, where n x , n y , n z are each integers, ranging from −∞ to ∞. The term in Eq. 10 of Zωn 2 is considered the ZP radiation contribution, since as T → 0, the second term of Zωn

vanishes. This second
term is what Planck concentrated his efforts upon, and of course is connected to the Planck spectrum. We will refer to Eq. 10 as being due to the ZP plus Planckian spectrum, or as the ZPP spectrum. Now, we are in a position to calculate the probability distributions of E(x, t) and B(x, t) in Eqs. 1 and 2, as well as consider much more complicated joint probabilities involving E(x, t) and B(x, t). We will carry this analysis out now; again, in Section 3, we will apply these ideas to electric dipole oscillators in SED.
We start by calculating the probability density function of realizing a specific value of the electric field, for the ZPP situation. Our ensemble varies of course due to its ensemble members, meaning by the ensemble distribution of the A ′ s and B ′ s in Eqs. 1 and 2 each time a new radiation situation is considered, then new A ′ s and B ′ s are realized according to the probability density distribution in Eq. 8, that then remain of constant values over the course of the subsequent physical analysis involving charged particles and fields.
As a start, the probability density distribution at position x and time t in Eq. 1 is where A ′ s and B ′ s here are symbolically written to represent the coefficients in Eq. 1, but as labeled there by A k n ,λ and B k n ,λ . In Eq. 11, P(A 1 , . . . , A N , B 1 , . . . , B N ) represents the probability density function of all these coefficients. In the end, we would let N → ∞. By E ZPP (x, t) in the Dirac delta function, we mean Eq. 1, but where the ZPP conditions of Eqs. 8-9 hold.
The key variables being integrated over in Eq. 11 are A kn,λ and B kn,λ variables. The integrations from −∞ to +∞ cover the range of their full possible values, while P(A 1 , . . . , A N , B 1 , . . . , B N ) provides the probability density associated with those values, and selects the values such that the probability density P(E at x, t) arises from all the possible matches of E ZPP (x, t) to the electric field value in question of E at x, t. As a side comment, in a sense, Eq. 11 has something in common with the Feynman path integral method in QM and QED, as the latter integrates over all weighted "path" contributions of a wave function evolving from one state to another. In contrast, Eq. 11 considers all the "paths," or allowed values of the A kn,λ and B kn,λ coefficients in the ensemble of radiation possibilities, that result in the condition E at x, t. While Eq. 11 directly involves probabilities, the Feynman path integral involves the QM wave function Ψ, with Ψ| 2 more indirectly providing the probability aspect.
Returning back to our present calculation involving Eq. 11, if either n ≠ n ′ , or λ ≠ λ ′ , then A k n ,λ and A kn ′ ,λ′ represent independent random variables, as do B k n ,λ and B kn ′ ,λ′ ; moreover, A k n ,λ and B kn ′ ,λ′ are also independent random variables, even when n n ′ and λ λ ′ . Using the Fourier representation for the Dirac delta function in Eq. 11 of in addition to the Gaussian distribution in Eq. 8, then Eq. 11 becomes: where to simplify notation, ω n and T will be suppressed here: To evaluate Eq. 13, Eq. 1 needs to be substituted in three places on the last line. To simplify notation yet again, let us replace Eq. 1 via where q represents all the indices of n x , n y , n z , λ, with their appropriate ranges, A q still represents A kn,λ , and likewise for B q and B kn,λ , and we will also refer to σ 2 n as σ 2 q from here on, again to simplify notation. Also, in Eq. 15, the expression has been abbreviated using and Hence: These integrals can be done by completing the squares of the A q and B q variables, then integrating over the resulting Gaussian expressions, followed by the integrals over s 1 , s 2 , s 3 . For example, completing the square: results in: since the Gaussian integral in the second line equals unity.
Continuing for each A q and B q results in: Frontiers in Physics | www.frontiersin.org April 2021 | Volume 8 | Article 580869 5

Cole
Probability Calculations in SED Three integrals remain, namely, over s x , s y , s z . The arguments of the two exponential terms in the second and third lines of Eq. 21, become, after making use of Eqs. 16 and 17: Still, the three integrals in Eq. 21 are nontrivial to evaluate, because of the cross terms in Eq. 22. However, the integrals can be greatly simplified by first summing over the polarization indices of λ 1, 2 as part of the "q" set of indices, and making use of the following identities for the three perpendicular unit vectors of ε k n ,1 , ε k n ,2 , and k n : After summing over the λ part of the q indices, one obtains: Although the "cross terms" involving s x s y , s x s z , s y s z still remain, upon summing over n in the last three terms, we obtain since n x , n y , n z each vary as integers symmetrically from −∞ to +∞, where k n is given in Eq. 3.
Consequently, from Eqs. 21, 25 and 26: where (28) Simplifying notation, let By then completing the square in Eq. 28 and carrying out the integral, yields More insight into Eq. 31 can be gained by relating α i in Eq. 29, to 〈E 2 ZPP,i (x, t)〉, using Eqs. 6, 7, and 23: Combining Eqs. 29, 31 and 32 and noting we obtain: Note that the mathematically detailed development of Eq. 34, and shortly Eq. 37 for the magnetic field case, agree nicely with the less detailed, but still the same result from Ref. [16].
In the symmetrical situation, with L x L y L z chosen for the rectilinear parallelepiped, then Thus, the probability density for the radiation electric field E at position x and time t, from either Eqs. 34 or 36, equals the product of three Gaussian functions. Moreover, the probability density of E is independent of position x and time t. If material walls existed, as in a cavity of arbitrary shape, as opposite to this free space situation treated by periodic bcs, then the probability density function for the fields could well be dependent on position. In Planck's original treatment of blackbody radiation [25], he considered a cavity with smooth walls and a size that was large compared to the key wavelengths of interest.
Since his work, researchers have probed on variations of these concerns, including small cavities, often referred to as the areas of quantum cavity electrodynamics [33,34]. To treat these problems in SED, one would need to take into account the precise nodal structure due to the cavity shape, and likely not take continuum approximation limits.
Looking back at the calculations, it is fairly easy to show that when analyzing magnetic fields, but now using Eq. 2, that: Finally, a caution needs to be made upon understanding Eqs. 32-37. When ZP radiation is included in the analysis, which is indeed a cornerstone of SED, 〈E 2 ZPP,i (x, t)〉 is infinite, as the energy spectrum monotonically grows with larger values of frequency. If one only considers the Planckian part of the spectrum, then this infinity does not happen. However, for calculating quantities like Casimir forces, van der Waals forces, and prevention of hydrogen collapse, it is absolutely essential to include the ZP spectrum. Cutoffs of the spectrum have been considered, but to date, the usual treatment has been to examine changes in regions between material boundaries, such as plates or cavity walls, when these are displaced. Such changes in energy due to wall displacements are finite, even with ZP fields [2]. Moreover, the results agree with experiments carried out to date. However, when considering the calculation in the next section, involving the probability of the electric field at two different positions and/or times, this infinity problem does not occur, unless the two points are chosen to be the same point both in space and time, or if |Δx| c|Δt|.

Joint Probability Density for Two Electric Field Values
The method just used can in principle be carried out for a wide range of probabilistic situations, with the key starting point being a similar condition to Eq. 11. A second calculation for a more complicated situation will be carried out here to illustrate this point. Before beginning, it is interesting to note that a number of "two-point" correlation functions of fields in SED have been calculated before by researchers, with the key reference being [24], but also [20], as well as by the present author in Ref. [35] and even for two-point correlation functions involving points in space and time following uniformly accelerated trajectories [7]. Clearly, probability density distributions such as in Eqs. 34 and 37 are more general, since they can be used to deduce all possible moments of the probability distribution; however, their calculation is in general much more involved.
Here, we will calculate the joint probability density function of where the semicolon in the first line is intended as a shortened meaning for the logical "AND" symbol of ∩ . Again, we make use of the random variable independence of the A ′ s and B ′ s for a normal thermodynamic radiation situation, and impose the distribution Eq. 8, plus use Eq. 12, to obtain: where the positions x 1 and x 2 and times t 1 and t 2 are contained in the radiation field expressions of Eq. 1 or 15. Thus, in Eq. 40, E 1x,ZPP , E 1y,ZPP , and E 1z,ZPP refer to E ZPP (x 1 , t 1 ) as in Eq. 15, and similarly for E ZPP (x 2 , t 2 ). Again abbreviating expressions, as in 16 and 17, with a 1, 2, as below: Collecting the A q -related terms in Eq. 40 as in the following manner, then later doing similarly for the B q terms: Now integrating over each A q term, followed later by integrating over the related B q term expression, results in: (43) Integrating over each A q and B q results in: sq,x + s 1y E 1,sq,y + s 1z E 1,sq,z + s 2x E 2,sq,x + s 2y E 2,sq,y + s 2z E 2,sq,z 2 σ 2 q 2 ⎤ ⎥ ⎥ ⎦ (44) Six integrals remain, namely, over s 1x , s 1y , s 1z , s 2x , s 2y , and s 2z . The arguments of the exponential terms in the last lines of Eq. 44, using Eq. 41 as well as Eqs. 16 and 17, and recognizing that many of the terms below have the simplification factor of cos 2 (k n · x a − ω n t a ) + sin 2 (k n · x a − ω n t a ) 1 , (45) for either a 1 or 2, then: s 1x E 1,cq,x + s 1y E 1,cq,y + s 1z E 1,cq,z + s 2x E 2,cq,x + s 2y E 2,cq,y + s 2z E 2,cq,z 2 + s 1x E 1,sq,x + s 1y E 1,sq,y + s 1z E 1,sq,z + s 2x E 2,sq,x + s 2y E 2,sq,y + s 2z E 2,sq,z 2 1 L x L y L z s 2 1x ε 2 q,x + s 2 1y ε 2 q,y + s 2 1z ε 2 q,z + s 2 2x ε 2 q,x + s 2 2y ε 2 q,y + s 2 2z ε 2 q,z + 2 L x L y L z s 1x ε q,x s 1y ε q,y + s 1x ε q,x s 1z ε q,z + s 1y ε q,y s 1z ε q,z + 2 L x L y L z s 2x ε q,x s 2y ε q,y + s 2x ε q,x s 2z ε q,z + s 2y ε q,y s 2z ε q,z + 2s 1x ε q,x L x L y L z s 2x ε q,x + s 2y ε q,y + s 2z ε q,z (C 1 C 2 + S 1 S 2 ) + 2s 1y ε q,y L x L y L z s 2x ε q,x + s 2y ε q,y + s 2z ε q,z (C 1 C 2 + S 1 S 2 ) + 2s 1z ε q,z L x L y L z s 2x ε q,x + s 2y ε q,y + s 2z ε q,z (C 1 C 2 + S 1 S 2 ) The meaning of the abbreviated terms (C 1 C 2 + S 1 S 2 ) is: Again summing over the polarization indices of λ 1, 2 as part of the "q" set of indices in the last lines of Eq. 44, and making use Eqs. 23 and 24, plus noting that the cross terms of ε q,i ε q,j for i ≠ j, will drop out due to the sum in Eq. 26, results in: q s 1x E 1,cq,x + s 1y E 1,cq,y + s 1z E 1,cq,z + s 2x E 2,cq,x + s 2y E 2,cq,y + s 2z E 2,cq,z 2 + q s 1x E 1,sq,x + s 1y E 1,sq,y + s 1z E 1,sq,z + s 2x E 2,sq,x + s 2y E 2,sq,y + s 2z E 2,sq,z Hence: n,x k 2 n + s 2 1y + 2s 1y s 2y C 12,ñ + s 2 2y ⎛ ⎝ 1 − k 2 n,y k 2 n ⎞ ⎠ + s 2 1z + 2s 1z s 2z C 12,ñ + s 2 Now to evaluate one of these integrals, it does not matter which we pick, as they all have the same form. Choosing I 1x2x , we can first complete the square in s 1x , then integrate over s 1x , followed by completing the square in s 2x , and then integrating over s 2x . Let and where C 12,n was defined in Eq. 47. Then: Completing the square with results in: (59) To now carry out the integration over s 2, B x in Eq. 56 must be expanded, as B x contains s 2 . Also, let so that (61) Then: Applying Eq. 58 again: Frontiers in Physics | www.frontiersin.org April 2021 | Volume 8 | Article 580869 As will be discussed in more detail in Section 2.4, our deduction of Eq. 62 is actually a multivariate normal (Gaussian) distribution involving E 1x and E 2x . Moreover, this distribution depends on the spatial and time differences, (x 1 − x 2 ) and (t 1 − t 2 ), between the two space time points, through the quantity C x . Moreover, since I 1x2x , I 1y2y , and I 1z2z will all have the same form as in Eq. 62, and the final probability density P(E 1 at x 1 , t 1 ; E 2 at x 2 , t 2 ) in Eq. 49 is just the product I 1x2x · I 1y2y · I 1z2z , then we will have obtained a multivariate normal distribution involving six field values at two points in space and time: E 1x , E 1y , E 1z , E 2x , E 2y , and E 2z . Again, these points will be made clearer in Section 2.4.

Checks on Behavior of I 1x2x
To be a probability density for E 1x and E 2x , as given by I 1x2x , certain probabilistic properties must hold. We will examine some of them here, such as should equal unity. Checking: so this is fine.

Another check is whether
and of course the opposite situation of ∞ −∞ dE 1i I 1i2i P(E 2i ). In Eq. 64, we have already deduced from earlier work, including Eq. 34, that As shown earlier, Eq. 65 turns out to be independent of x 1 , t 1 , where i 1, 2, 3 refers to x, y, z, respectively. Mathematically, this independence on x 1 , t 1 arises because 〈E 2 ZPP,i 〉 is independent of x 1 , t 1 , as seen in Eq. 32. A more "physical" view of this result is that the stochastic properties of the ZP and ZPP fields are homogeneous and isotropic in space and independent of time origin. In any case, from Eqs. 32, 54, and 55, and if we generalize A x to A i for i 1, 2, and 3, to include all three x, y, z cases, then: Returning to Eq. 64 and using 50: Ci Ai In the integral on the far right, bottom line, let E 1x − E 2x Cx Ax u and du dE 1x : since the integral is odd in E 2x . Now trying a check that is more involved, we will compute a 2point correlation function in the fields, using the probability density I 1x2x . We will calculate an example, along the lines of Ref. [24], but also carrying the calculations to a final analytical expression, more along [7,35]. Hence, this will be an example, with of course many other two point sets of coordinates in time and space that could be carried out, but even this single example is nontrivial to carry out. Most importantly, however, this example shows how to carry out the analysis in Ref. [24] via the probability density method discussed here.
Following roughly along Refs. [24,35], and [7], then: Making the change of discrete to continuous variables, with k n Ly y + 2πn z Lz z, for large values of L x , L y , L z , with enables integrals to be carried out. Moreover, although the ZPP spectrum in the integral could be evaluated, the ZP spectrum is certainly much easier to do so analytically. Since this is just an example, we will proceed with restriction to the ZP case, or [σ(ω, T)] 2 → 2πZω 2πZkc: × cos k y R cos(ωt) + sin k y R sin(ωt) .
The second term of sin(k y R)sin(ωt) makes the integrand odd in k y . Hence: cos k y R cos(kct).
Next, we will make k y be the axis where the polar angle is measured from, so that k y kcosθ. Consequently, k x ksinθsinϕ and k z ksinθcosϕ. Hence: Consequently, Eq. 72 becomes: Substituting ct R b, and using an integral table [36] (p. 504, No. 8), also discussed in the limiting sense in Ref. [7], Appendix C: enables Eq. 76 to be evaluated: As mentioned earlier, unless R ct, this two-point correlation function is not singular.
The above calculation has typically been, roughly, the means for calculating such "two-point correlation" functions in SED, or even "npoint correlation functions" [24]. We will proceed to calculate the same quantity as in Eq. 81, but by using the joint probability density function for two electric field values, I 1x2x , Eq. 62, deduced in Section 2.2. Of course the two results should agree, but it is interesting to see the difference in methods.
Here, the meaning of the subscript at the end of x 1 , t 1 ; x 2 , t 2 is that the two electric field points E 1x and E 2x are to be evaluated at the two space and time points, x 1 , t 1 and x 2 , t 2 , respectively, in the function C 12,n , Eq. 47, contained within C x , in Eq. 60.
The integral over E 1x on the right can be broken up as: The first integral equals zero, as it is odd in x. Hence: 1 L x L y L z nx ,ny ,nz −∞ ∞ C 12,n x1,t1;x2,t2 K 2 n,x σ 2 n .
The last expression for C x came from Eq. 60. Hence, using Eqs. 47 and 54 replacing the coordinates x 1 and t 1 with 0, 0, and x 2 and t 2 with yR and t, respectively, and σ 2 n with ZP of 2πZω n,λ : Implementing the same continuum approximation as with the other method, Eq. 70, then results in: The second integral with sin(k y R) is odd in k y and equals zero by symmetry. Hence: However, this result agrees exactly with the earlier result, Eq. 71, obtained partway through the ensemble derivation of 〈E ZP,x (0, 0)E ZP,x ( yR, t)〉. Thus, continuing with further steps in evaluating Eq. 83 will provide a final result, using the joint probability density approach of

Multivariate Normal Distribution
Much of the work carried out here can be generalized using the multivariate normal distribution. The two key expressions for us here are the Fourier decomposition of the radiation fields in Eqs. 1 and 2 not just because they are Fourier decompositions, but also that they are a linear sum of the random variables A kn,λ and B kn,λ . To put this in better perspective, if we imagine an ensemble of boxes L x × L y × L z , each the same size, but existing at different points in space and/or in time, then electromagnetic field fluctuations of E ZPP and B ZPP will occur at each point within each box. However, for each box, there is only one set of coefficients A k n ,λ and B k n ,λ , as these coefficients do not change from the initial point of field evolution. However, as viewed over the entire ensemble of boxes, the coefficients are assumed to be independent random variables obeying Gaussian distributions. The mean for each, over the ensemble, is zero, as in Eq. 4, and the normal distribution for either A kn,λ or B kn,λ is given by Eq. 8, while the variance of each is given by [σ(ω n , T)] 2 , as in Eq. 8.
Thus, E ZPP and B ZPP can be viewed as the linear transformation of the random variables A k n ,λ and B k n ,λ . However, the coefficients multiplying A k n ,λ and B k n ,λ in these linear sums in Eqs. 1 and 2 are not constants, as they depend on time and space. In particular, 1 (LxLyLz) 1/2 ε k n ,λ cos(k n · x − ω n t) and 1 (LxLyLz) 1/2 ε k n ,λ sin(k n · x − ω n t) are the coefficients multiplying the random variables of A k n ,λ and B k n ,λ for E(x, t) in Eq. 1. An exactly similar situation occurs for A kn,λ and B kn,λ regarding B(x, t) in Eq. 2, except that ε kn,λ is replaced by k n × ε kn,λ . Thus, although the probabilistic properties for the random variables A kn,λ and B kn,λ are independent of time and space, as expressed by Eq. 8, the same is not true for the probabilistic/stochastic properties of E ZPP and B ZPP , as seen for example in Eq. 81 and other related results discussed in this subsection.
A multivariate normal expression for the probability density of a set of field values would be represented by: . . . ; E n at x n , t n ; B n+1 at x n+1 , t n+1 ; B n+2 at x n+2 , t n+2 ; . . . ; B n+m at x n+m , t n+m , where there are n electric field vector values (i.e., 3 × n component values, as indicated below) and m magnetic field vector values (3 × m component values), at respective positions in space and time, as indicated. However, since this is a multivariate normal distribution, where all 〈E i 〉 and all 〈B i 〉 ensemble averages equal zero, then by probability theory, the above would be represented by [37]: where X is the vector of values. Moreover, Σ is the covariant matrix as expressed by since 〈X i 〉 0 for all E i and B i values, due to Eqs. 4, 1 and 2. Also, |Σ| and Σ −1 represent the determinant and inverse matrix of the covariant matrix, Σ, respectively. If we calculated all the components of the covariant matrix, meaning all combinations Σ ij 〈X i X j 〉 of pairs of electric and magnetic field expectation values, then the probability density, Eq. 85, could be evaluated for any vector X of electric and magnetic field components at different space and time points. Two comparisons can immediately be made with work already covered here. In Section 2.1, the probability density was deduced for P(E at x, t) in Eq. 34, but this also follows from Eq. 85 with X being the vector (E 1x , E 1y , E 1z ), then using Eq. 86 plus earlier relations in Section 2.1, and which leads to Eq. 34. Similarly, Eq. 85 can be used to deduce the probability density function for two electric field points that was covered in Section 2.2. For a multivariate normal distribution for two field points, although Eq. 85 is certainly the correct equation to use, it is usually rewritten in the following simplified form in probability textbooks, and referred to as the "bivariate" normal distribution [38]: Here, μ 1x ≡ 〈E 1x 〉 and μ 2x ≡ 〈E 2x 〉 are both zero for our ZP and ZPP cases. Also, σ 2 1x ≡ 〈(E 1x − μ 1x ) 2 〉 〈E 2 1x 〉, for our situation, and similarly for (σ 2x ) 2 〈(E 2x ) 2 〉, while is called the Pearson's correlation coefficient of E 1x and E 2x . With μ 1x μ 2x 0 for our case, then Putting these expressions into Eq. 87 results in: This is then readily related to our result of I 1x2x in Section 2.2, Eq. 62, since from Eq. 66, A i 〈E 2 ZPP,i 〉 for (i → x, y, z), and is independent of space and time, so 〈E 2 1x 〉 〈E 2 2x 〉, for example. Moreover, C x in Eq. 60 can be shown to be the two-point correlation function of the x component of the electric field at two different space/time points, or, 〈E 1x E 2x 〉. Thus, which is just ρ 1x2x in Eq. 89. Thus, Eq. 90 agrees with I 1x2x in Eq. 62.
Extending this result to the discussion at the end of Section 2.2 involving x, y, and z components of two electric field values: P(E 1 at x 1 , t 1 ; E 2 at x 2 , t 2 )

CONCLUDING REMARKS
The technique used here in Eq. 11 for one electric field value E, or three component values E x , E y , E z , or with Eq. 40 for P(E 1 at x 1 , t 1 ; E 2 at x 2 , t 2 ), for two electric field values, or six component values, was clearly understood to be extendable to n electric and/or m magnetic field values. Various tests and examinations were carried out in Section 2.3 to provide further understanding of the results derived here. Section 2.4 showed how to more easily use the multivariate normal distribution to obtain the same probability densities. However, although the methods of Eqs. 11 and 40 and the obvious generalizations to far more electric and magnetic radiation field values, result in long calculations, there are a few aspects that should be noted. As briefly discussed in Ref. [23], these techniques can readily be applied to the electric dipole simple harmonic oscillator, in either one, two, or three oscillatory degrees of freedom. For example, for a 1-D oscillator, P(x 1 at t 1 ; x 2 at t 2 ) represents the probability density of finding an oscillator extension of x 1 at time t 1 , and extension of x 2 at time t 2 . The subtle point here is the expressions for x(t 1 ) and x(t 2 ) must be inserted into Eq. 92, and these depend on the A n and B n values. Equation 92 is similar to Eq. 40, but more complicated. The electric dipole oscillator system, often phrased in terms of the simple harmonic oscillator (SHO), was a key early system that was studied in SED. However, the drawback is that the oscillator computations are even longer than in Sections 2.1-2.3, which only involved the probability states of electric radiation values in ZP and ZPP conditions.
Another interesting aspect of the present method is that, in principle, the method can begin to tackle other systems, particularly one that has not yet been solved analytically in SED: the classical hydrogen problem. Now, for that system, there is no simplification, such as the multivariate normal distribution, to provide a simpler method of solution. The classical hydrogen atom is not a linear SHO system; it is nonlinear, and the multivariate normal distribution only becomes possible for linear sums of random variables that each obey the normal distribution. Possibly, the classical hydrogen system in SED is intractable with the present method, but as far as the author knows, this approach has not yet been tried. This brings us back to brief comments made in Section 2.1 about the Feynman path integral for QM and QED, and how there is a slight connection to the present method for SED. Although Feynman developed the technique by 1948, it was initially applied only to some relatively simple systems, such as discussed in [39]. The hydrogen atom escaped solution by Feynman and others until about 1979, by Duru and Kleinert. In analogy, at first blush, although the present method for classical hydrogen in SED may be too difficult, still it seems an interesting perspective to consider.
Finally, a last comment: many of the computations shown here for field values could likely be extended to far more field points, simply by writing the correct code to make a versatile program of n-point correlation functions, such as with the aid of a symbolic mathematic program, as has been done in some cases for Feynman diagram calculations.