Probability Calculations Within Stochastic Electrodynamics

1 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–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 de-accelerated 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⁻¹¹ 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–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.

2 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\le 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–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 $\langle {\left(X-\langle X\rangle \right)}^N\rangle =\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}xP\left(x\right){\left(x-\langle x\rangle \right)}^N$, for $N=\mathrm{1,2,3},\dots$, where the angled brackets mean “expectation value,” and where $P\left(x\right)$ 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 $\mu \equiv \langle X\rangle$, and the variance, ${\sigma }^2\equiv \langle {\left(X-\langle X\rangle \right)}^2\rangle$. All other moments $\langle {\left(X-\langle X\rangle \right)}^N\rangle$, for $N=\mathrm{1,2,3},\dots$, 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 $\langle {\left(X-\langle X\rangle \right)}^N\rangle$, whereas here we will directly calculate quantities like $P\left(x\right)$.

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\ge 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 $\mathbf{E}\left(\mathbf{x},t\right)$ and magnetic $\mathbf{B}\left(\mathbf{x},t\right)$ radiation fields in this large parallelepiped volume can be written as the following sum of plane waves [16]:

$\mathbf{E}\left(\mathbf{x},t\right)=\frac{1}{{\left(L_x{L}_y{L}_z\right)}^{1/2}}{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }^{\infty }{\displaystyle \sum }_{\lambda =\mathrm{1,2}}{\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\left[A_{{\mathbf{k}}_{\mathbf{n}},\lambda }\text{cos}\left({\mathbf{k}}_{\mathbf{n}}\cdot \mathbf{x}-{\omega }_{\mathbf{n}}t\right)+B_{{\mathbf{k}}_{\mathbf{n}},\lambda }\text{sin}\left({\mathbf{k}}_{\mathbf{n}}\cdot \mathbf{x}-{\omega }_{\mathbf{n}}t\right)\right],(1)$

$\mathbf{B}\left(\mathbf{x},t\right)=\frac{1}{{\left(L_x{L}_y{L}_z\right)}^{1/2}}{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }^{\infty }{\displaystyle \sum }_{\lambda =\mathrm{1,2}}\left({\widehat{k}}_{\mathbf{n}}\times {\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\right)\left[A_{{\mathbf{k}}_{\mathbf{n}},\lambda }\text{cos}\left({\mathbf{k}}_{\mathbf{n}}\cdot \mathbf{x}-{\omega }_{\mathbf{n}}t\right)+B_{{\mathbf{k}}_{\mathbf{n}},\lambda }\text{sin}\left({\mathbf{k}}_{\mathbf{n}}\cdot \mathbf{x}-{\omega }_{\mathbf{n}}t\right)\right],(2)$

where

${\mathbf{k}}_{\mathbf{n}}=\frac{2\pi n_x}{L_x}\widehat{\mathbf{x}}+\frac{2\pi n_y}{L_y}\widehat{\mathbf{y}}+\frac{2\pi n_z}{L_z}\widehat{\mathbf{z}},(3)$

and $n_x$, $n_y$, and $n_z$ are integers, and ${\omega }_{\mathbf{n}}=c\left|{\mathbf{k}}_{\mathbf{n}}\right|$, ${\mathbf{k}}_{\mathbf{n}}\cdot {\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }={\mathbf{k}}_{\mathbf{n}}\cdot {\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},{\lambda }^{\prime }}=0$, and ${\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\cdot {\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},{\lambda }^{\prime }}=0$ for $\lambda \ne {\lambda }^{\prime }$, where λ and ${\lambda }^{\prime }$ indicate the linear polarization direction. Specifically, λ might be represented by the values 1 or 2, and the same for ${\lambda }^{\prime }$. Also, ${\widehat{\mathbf{k}}}_{\mathbf{n}}={\mathbf{k}}_{\mathbf{n}}/\left|{\mathbf{k}}_{\mathbf{n}}\right|$. Equations 1 and 2 satisfy the wave equations of ${\nabla }^2\mathbf{E}\left(\mathbf{x},t\right)=\frac{1}{c^2}\frac{{\partial }^2}{\partial t^2}\mathbf{E}\left(\mathbf{x},t\right)$ and ${\nabla }^2\mathbf{B}\left(\mathbf{x},t\right)=\frac{1}{c^2}\frac{{\partial }^2}{\partial t^2}\mathbf{B}\left(\mathbf{x},t\right)$, which can be deduced from Maxwell’s equations for free space (charge density and current charge density both equal to zero). Moreover, the presence of ${\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ and $\left({\widehat{\mathbf{k}}}_{\mathbf{n}}\times {\widehat{\mathbf{\epsilon }}}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\right)$ in Eqs. 1, 2 respectively, and the cited relationships of ${\mathbf{k}}_{\mathbf{n}}\cdot {\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }={\mathbf{k}}_{\mathbf{n}}\cdot {\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},{\lambda }^{\prime }}=0$, and ${\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\cdot {\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},{\lambda }^{\prime }}=0$ for $\lambda \ne {\lambda }^{\prime }$, provide the other needed relationships for satisfying the four Maxwell’s equations for free fields, such as Faraday’s law of $\nabla \times \mathbf{E}=-\frac{\partial }{c\partial t}\mathbf{B}$.

Following more or less the lead of Planck’s first half of Ref. [25], the above radiation fields represent the stochastic fluctuations of thermal radiation, for $0\le T$ (i.e., including $T=0$), with the following assumptions. The coefficients of the expressions for $\mathbf{E}\left(\mathbf{x},t\right)$ and $\mathbf{B}\left(\mathbf{x},t\right)$ in Eqs. 1 and 2, namely, $A_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ and $B_{{\mathbf{k}}_{\mathbf{n}},\lambda }$, were assumed to be randomly distributed in the following way initially, but once fixed, they stay fixed in all physical analysis, such as in the interaction of charged oscillators and radiation, as in simulations of [28–31]. In the case of simulations, the physical 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_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ and $B_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ 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:

$\langle A_{{\mathbf{k}}_{\mathbf{n}},\lambda }\rangle =\langle B_{{\mathbf{k}}_{\mathbf{n}},\lambda }\rangle =0.(4)$

Moreover, the “A and B coefficients” are considered to be both independent and uncorrelated random variables in this ensemble, so

$\langle A_{{\mathbf{k}}_{\mathbf{n}},\lambda }{B}_{{\mathbf{k}}_{{\mathbf{n}}^{\prime }},{\lambda }^{\prime }}\rangle =0,(5)$

as are the “A coefficients” with different indices, and the same for the “B coefficients”:

$\langle A_{{\mathbf{k}}_{\mathbf{n}},\lambda }{A}_{{\mathbf{k}}_{{\mathbf{n}}^{\prime }},{\lambda }^{\prime }}\rangle =\langle B_{{\mathbf{k}}_{\mathbf{n}},\lambda }{B}_{{\mathbf{k}}_{{\mathbf{n}}^{\prime }},{\lambda }^{\prime }}\rangle =0,\text{if }\mathbf{n}\ne {\mathbf{n}}^{\prime }\text{ or }\lambda \ne {\lambda }^{\prime }.(6)$

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:

$\langle A_{{\mathbf{k}}_{\mathbf{n}},\lambda }{A}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\rangle =\langle B_{{\mathbf{k}}_{\mathbf{n}},\lambda }{B}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\rangle ={\left[\sigma \left({\omega }_{\mathbf{n}},T\right)\right]}^2.(7)$

Although Eqs. 4–7 are indeed assumed in SED, the more general relationship that includes all these relationships, plus more, is that the $A^{\prime }s$ and $B^{\prime }s$ coefficients are assumed to be independent random variables, with zero mean as in Eq. 4, with variance ${\left[\sigma \left({\omega }_{\mathbf{n}},T\right)\right]}^2,$ and with probability density distributions characterizing the ensemble of possible radiation situations characterizing thermal radiation at temperature T, for $0\le T$, as being Gaussian distributions. Specifically:

$P\left(A_{{\mathbf{k}}_{\mathbf{n}},\lambda }\right)=\frac{1}{\sqrt{2\pi {\left[\sigma \left({\omega }_{\mathbf{n}},T\right)\right]}^2}}\text{exp}\left\{-\frac{1}{2}{\left[\frac{A_{{\mathbf{k}}_{\mathbf{n}},\lambda }}{\sigma \left({\omega }_{\mathbf{n}},T\right)}\right]}^2\right\},(8)$

with the same also holding for $P\left(B_{{\mathbf{k}}_{\mathbf{n}},\lambda }\right)$. 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 ${\left[\sigma \left({\omega }_{\mathbf{n}},T\right)\right]}^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 ${\left[\sigma \left({\omega }_{\mathbf{n}},T\right)\right]}^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:

${\left[\sigma \left({\omega }_{\mathbf{n}},T\right)\right]}^2=2\pi \hslash {\omega }_{\mathbf{n}}+\frac{4\pi \hslash {\omega }_{\mathbf{n}}}{\text{exp}\left(\frac{\hslash {\omega }_{\mathbf{n}}}{k_{\text{B}}T}\right)-1}=2\pi \hslash {\omega }_{\mathbf{n}}\text{coth}\left(\frac{\hslash {\omega }_{\mathbf{n}}}{2k_{\text{B}}T}\right).(9)$

Note: $\underset{T\to 0}{\text{lim}}\text{coth}\left(\frac{\hslash {\omega }_{\mathbf{n}}}{2k_{\text{B}}T}\right)=1$. Consequently, the term $2\pi \hslash {\omega }_{\mathbf{n}}$ constitutes the ZP $\left(T\to 0\right)$ spectrum contribution, while the $\frac{4\pi \hslash {\omega }_{\mathbf{n}}}{\text{exp}\left(\frac{\hslash {\omega }_{\mathbf{n}}}{k_{\text{B}}T}\right)-1}$ 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\le T$, within the $L_x\times L_y\times 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

$\mathrm{ℰ}={{\displaystyle \int }}^{\text{}}\text{d}V\frac{1}{8\pi }\langle {\mathbf{E}}^2\left(\mathbf{x},t\right)+{\mathbf{B}}^2\left(\mathbf{x},t\right)\rangle ={\displaystyle \sum }_{\mathbf{n}}\left[\frac{\hslash {\omega }_{\mathbf{n}}}{2}+\frac{\hslash {\omega }_{\mathbf{n}}}{\text{exp}\left(\frac{\hslash {\omega }_{\mathbf{n}}}{k_{\text{B}}T}\right)-1}\right],(10)$

where ${\omega }_{\mathbf{n}}=c\left|{\mathbf{k}}_{\mathbf{n}}\right|$ follows from Eq. 3; $\mathbf{n}$ is composed of $\left\{n_x,n_y,n_z\right\},$ where $n_x,n_y,n_z$ are each integers, ranging from $-\infty$ to $\infty$. The term in Eq. 10 of $\frac{\hslash {\omega }_{\mathbf{n}}}{2}$ is considered the ZP radiation contribution, since as $T\to 0$, the second term of $\frac{\hslash {\omega }_{\mathbf{n}}}{\text{exp}\left(\frac{\hslash {\omega }_{\mathbf{n}}}{k_{\text{B}}T}\right)-1}$ 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 $\mathbf{E}\left(\mathbf{x},t\right)$ and $\mathbf{B}\left(\mathbf{x},t\right)$ in Eqs. 1 and 2, as well as consider much more complicated joint probabilities involving $\mathbf{E}\left(\mathbf{x},t\right)$ and $\mathbf{B}\left(\mathbf{x},t\right)$. 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^{\prime }s$ and $B^{\prime }s$ in Eqs. 1 and 2 each time a new radiation situation is considered, then new $A^{\prime }s$ and $B^{\prime }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 $\mathbf{x}$ and time t in Eq. 1 is

$\begin{array}{l}P\left(\mathbf{E}\text{ at }\mathbf{x},t\right)\\ =\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{A}_1\cdots \stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{A}_N\cdots \stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{B}_1\cdots \stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{B}_N\cdots P\left(A_1,\dots ,A_N,\dots ,B_N,\dots \right){\delta }^3\left[\mathbf{E}-{\mathbf{E}}_{\text{ZPP}}\left(\mathbf{x},t\right)\right],\end{array}(11)$

where $A^{\prime }s$ and $B^{\prime }s$ here are symbolically written to represent the coefficients in Eq. 1, but as labeled there by $A_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ and $B_{{\mathbf{k}}_{\mathbf{n}},\lambda }$. In Eq. 11, $P\left(A_1,\dots ,A_N,B_1,\dots ,B_N\right)$ represents the probability density function of all these coefficients. In the end, we would let $N\to \infty$. By ${\mathbf{E}}_{\text{ZPP}}\left(\mathbf{x},t\right)$ 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_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ and $B_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ variables. The integrations from $-\infty$ to $+\infty$ cover the range of their full possible values, while $P\left(A_1,\dots ,A_N,B_1,\dots ,B_N\right)$ provides the probability density associated with those values, and ${\delta }^3\left[\mathbf{E}-{\mathbf{E}}_{\text{ZPP}}\left(\mathbf{x},t\right)\right]$ selects the values such that the probability density $P\left(\mathbf{E}\text{ at }\mathbf{x},t\right)$ arises from all the possible matches of ${\mathbf{E}}_{\text{ZPP}}\left(\mathbf{x},t\right)$ to the electric field value in question of $\mathbf{E}$ at $\mathbf{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_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ and $B_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ coefficients in the ensemble of radiation possibilities, that result in the condition $\mathbf{E}$ at $\mathbf{x},t$. While Eq. 11 directly involves probabilities, the Feynman path integral involves the QM wave function $\text{Ψ}$, with $|\text{Ψ}{|}^2$ more indirectly providing the probability aspect.

Returning back to our present calculation involving Eq. 11, if either $\mathbf{n}\ne {\mathbf{n}}^{\prime },$ or $\lambda \ne {\lambda }^{\prime }$, then $A_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ and $A_{{\mathbf{k}}_{\mathbf{n}\prime },{\lambda }^{\prime }}$ represent independent random variables, as do $B_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ and $B_{{\mathbf{k}}_{\mathbf{n}\prime },{\lambda }^{\prime }}$; moreover, $A_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ and $B_{{\mathbf{k}}_{\mathbf{n}\prime },{\lambda }^{\prime }}$ are also independent random variables, even when $\mathbf{n}={\mathbf{n}}^{\prime }$ and $\lambda ={\lambda }^{\prime }$. Using the Fourier representation for the Dirac delta function in Eq. 11 of

${\delta }^3\left[\mathbf{E}-{\mathbf{E}}_{\text{ZPP}}\left(\mathbf{x},t\right)\right]=\frac{1}{2\pi }\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}{s}_x{e}^{i{s}_x\left(E_x-E_{x,\text{ZPP}}\right)}\frac{1}{2\pi }\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}{s}_y{e}^{i{s}_y\left(E_y-E_{y,\text{ZPP}}\right)}\frac{1}{2\pi }\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}{s}_z{e}^{i{s}_z\left(E_z-E_{z,\text{ZPP}}\right)},(12)$

in addition to the Gaussian distribution in Eq. 8, then Eq. 11 becomes:

$\begin{array}{l}P\left(\mathbf{E}\text{ at }\mathbf{x},t\right)=\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{A}_1\cdots \stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{A}_N\dots \frac{1}{\sqrt{2\pi {\sigma }_{{\mathbf{n}}_1}^2}}\text{exp}\left[-\frac{{\left(A_1\right)}^2}{2{\sigma }_{{\mathbf{n}}_1}^2}\right]\dots \frac{1}{\sqrt{2\pi {\sigma }_{{\mathbf{n}}_N}^2}}\text{exp}\left[-\frac{{\left(A_N\right)}^2}{2{\sigma }_{{\mathbf{n}}_N}^2}\right]\dots \\ \times \stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{B}_1\cdots \stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{B}_N\dots \frac{1}{\sqrt{2\pi {\sigma }_{{\mathbf{n}}_1}^2}}\text{exp}\left[-\frac{{\left(B_1\right)}^2}{2{\sigma }_{{\mathbf{n}}_1}^2}\right]\dots \frac{1}{\sqrt{2\pi {\sigma }_{{\mathbf{n}}_N}^2}}\text{exp}\left[-\frac{{\left(B_N\right)}^2}{2{\sigma }_{{\mathbf{n}}_N}^2}\right]\dots \\ \times \frac{1}{2\pi }\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_x{e}^{i{s}_x\left(E_x-E_{x,\text{ZPP}}\right)}\frac{1}{2\pi }\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_y{e}^{i{s}_y\left(E_y-E_{y,\text{ZPP}}\right)}\frac{1}{2\pi }\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_z{e}^{i{s}_z\left(E_z-E_{z,\text{ZPP}}\right)},\end{array}(13)$

where to simplify notation, ${\omega }_{\mathbf{n}}$ and T will be suppressed here:

${\left[\sigma \left({\omega }_{\mathbf{n}},T\right)\right]}^2\equiv {\sigma }_{\mathbf{n}}^2.(14)$

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. 1via

${\mathbf{E}}_{\text{ZPP}}\left(\mathbf{x},t\right)={\displaystyle \sum }_q{A}_q{\mathbf{E}}_{cq}+{\displaystyle \sum }_q{B}_q{\mathbf{E}}_{sq},(15)$

where q represents all the indices of $n_x,n_y,n_z,\lambda$, with their appropriate ranges, $A_q$ still represents $A_{{\mathbf{k}}_{\mathbf{n}},\lambda }$, and likewise for $B_q$ and $B_{{\mathbf{k}}_{\mathbf{n}},\lambda }$, and we will also refer to ${\sigma }_{\mathbf{n}}^2$ as ${\sigma }_q^2$ from here on, again to simplify notation. Also, in Eq. 15, the expression has been abbreviated using

${\mathbf{E}}_{cq}\equiv \frac{1}{{\left(L_x{L}_y{L}_z\right)}^{1/2}}{\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\text{cos}\left({\mathbf{k}}_{\mathbf{n}}\cdot \mathbf{x}-{\omega }_{\mathbf{n}}t\right),(16)$

and

${\mathbf{E}}_{sq}\equiv \frac{1}{{\left(L_x{L}_y{L}_z\right)}^{1/2}}{\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\text{sin}\left({\mathbf{k}}_{\mathbf{n}}\cdot \mathbf{x}-{\omega }_{\mathbf{n}}t\right).(17)$

Hence:

$\begin{array}{l}P\left(\mathbf{E}\text{ at }\mathbf{x},t\right)={{\displaystyle \int }}^{\text{}}\text{d}{A}_1\cdots {{\displaystyle \int }}^{\text{}}\text{d}{A}_N\dots \frac{1}{\sqrt{2\pi {\sigma }_1^2}}\text{exp}\left[-\frac{{\left(A_1\right)}^2}{2{\sigma }_1^2}\right]\dots \frac{1}{\sqrt{2\pi {\sigma }_N^2}}\text{exp}\left[-\frac{{\left(A_N\right)}^2}{2{\sigma }_N^2}\right]\dots \\ \times {{\displaystyle \int }}^{\text{}}\text{d}{B}_1\cdots {{\displaystyle \int }}^{\text{}}\text{d}{B}_N\dots \frac{1}{\sqrt{2\pi {\sigma }_1^2}}\text{exp}\left[-\frac{{\left(B_1\right)}^2}{2{\sigma }_1^2}\right]\dots \frac{1}{\sqrt{2\pi {\sigma }_N^2}}\text{exp}\left[-\frac{{\left(B_N\right)}^2}{2{\sigma }_N^2}\right]\dots \\ \times \frac{1}{2\pi }\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_x{e}^{i{s}_x\left(E_x-{\displaystyle \sum }_q\left(A_q{E}_{cq,x}+B_q{E}_{sq,x}\right)\right)}\times \frac{1}{2\pi }\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_y{e}^{i{s}_y\left(E_y-{\displaystyle \sum }_q{A}_q{E}_{cq,y}+B_q{E}_{sq,y}\right)}\\ \times \frac{1}{2\pi }\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_z{e}^{i{s}_z\left(E_z-{\displaystyle \sum }_q{A}_q{E}_{cq,z}+B_q{E}_{sq,z}\right)}.\end{array}(18)$

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:

$\begin{array}{l}-\frac{{\left(A_q\right)}^2}{2{\sigma }_q^2}-i{s}_x{A}_q{E}_{cq,x}-i{s}_y{A}_q{E}_{cq,y}-i{s}_z{A}_q{E}_{cq,z}\\ =-\frac{1}{2{\sigma }_q^2}{\left[A_q+i\left(s_x{E}_{cq,x}+s_y{E}_{cq,y}+s_z{E}_{cq,z}\right){\sigma }_q^2\right]}^2-{\left(s_x{E}_{cq,x}+s_y{E}_{cq,y}+s_z{E}_{cq,z}\right)}^2\frac{{\sigma }_q^2}{2}\end{array}(19)$

results in:

$\begin{array}{l}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{A}_q\frac{1}{\sqrt{2\pi {\sigma }_q^2}}\text{exp}\left[-\frac{{\left(A_q\right)}^2}{2{\sigma }_q^2}\right]e^{-i{s}_x{A}_q{E}_{cq,x}}{e}^{-i{s}_y{A}_q{E}_{cq,y}}{e}^{-i{s}_z{A}_q{E}_{cq,z}}\\ =\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{A}_q\frac{1}{\sqrt{2\pi {\sigma }_q^2}}\text{exp}\left\{-\frac{1}{2{\sigma }_q^2}{\left[A_q+i\left(s_x{E}_{cq,x}+s_y{E}_{cq,y}+s_z{E}_{cq,z}\right){\sigma }_q^2\right]}^2\right\}\\ \times \exp \left[-{\left(s_x{E}_{cq,x}+s_y{E}_{cq,y}+s_z{E}_{cq,z}\right)}^2\frac{{\sigma }_q^2}{2}\right]\\ =\text{exp}\left[-{\left(s_x{E}_{cq,x}+s_y{E}_{cq,y}+s_z{E}_{cq,z}\right)}^2\frac{{\sigma }_q^2}{2}\right],\end{array}(20)$

since the Gaussian integral in the second line equals unity.

Continuing for each $A_q$ and $B_q$ results in:

$\begin{array}{l}P\left(\mathbf{E}\text{ at }\mathbf{x},t\right)=\frac{1}{2\pi }\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_x{e}^{i{s}_x{E}_x}\frac{1}{2\pi }\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_y{e}^{i{s}_y{E}_y}\frac{1}{2\pi }\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_z{e}^{i{s}_z{E}_z}\\ \times \exp \left[-{\displaystyle \sum }_q{\left(s_x{E}_{cq,x}+s_y{E}_{cq,y}+s_z{E}_{cq,z}\right)}^2\frac{{\sigma }_q^2}{2}\right]\\ \times \exp \left[-{\displaystyle \sum }_q{\left(s_x{E}_{sq,x}+s_y{E}_{sq,y}+s_z{E}_{sq,z}\right)}^2\frac{{\sigma }_q^2}{2}\right].\end{array}(21)$

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:

$\begin{array}{l}{\left(s_x{E}_{cq,x}+s_y{E}_{cq,y}+s_z{E}_{cq,z}\right)}^2+{\left(s_x{E}_{sq,x}+s_y{E}_{sq,y}+s_z{E}_{sq,z}\right)}^2\\ =\frac{1}{\left(L_x{L}_y{L}_z\right)}\left(s_x^2{\epsilon }_{q,x}^2+s_y^2{\epsilon }_{q,y}^2+s_z^2{\epsilon }_{q,z}^2+2s_x{s}_y{\epsilon }_{q,x}{\epsilon }_{q,y}+2s_x{s}_z{\epsilon }_{q,x}{\epsilon }_{q,z}+2s_y{s}_z{\epsilon }_{q,y}{\epsilon }_{q,z}\right).\end{array}(22)$

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 $\lambda =\mathrm{1,2}$ as part of the “q” set of indices, and making use of the following identities for the three perpendicular unit vectors of ${\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},1}$, ${\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},2}$, and ${\mathbf{k}}_{\mathbf{n}}$:

${\displaystyle \sum }_{\lambda =\mathrm{1,2}}{\left[{\left({\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\right)}_i\right]}^2=1-{\left[\frac{{\left({\mathbf{k}}_{\mathbf{n}}\right)}_i}{k_{\mathbf{n}}}\right]}^2,(23)$

${\displaystyle \sum }_{\lambda =\mathrm{1,2}}{\left({\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\right)}_i{\left({\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\right)}_j={\delta }_{ij}-\frac{{\left({\mathbf{k}}_{\mathbf{n}}\right)}_i{\left({\mathbf{k}}_{\mathbf{n}}\right)}_j}{k_{\mathbf{n}}^2}.(24)$

After summing over the λ part of the q indices, one obtains:

$\begin{array}{l}{\displaystyle \sum }_q\left(s_x^2{\epsilon }_{q,x}^2+s_y^2{\epsilon }_{q,y}^2+s_z^2{\epsilon }_{q,z}^2+2s_x{s}_y{\epsilon }_{q,x}{\epsilon }_{q,y}+2s_x{s}_z{\epsilon }_{q,x}{\epsilon }_{q,z}+2s_y{s}_z{\epsilon }_{q,y}{\epsilon }_{q,z}\right){\sigma }_q^2\\ ={\displaystyle \sum }_{\mathbf{n}}\left\{s_x^2\left[1-{\left(\frac{k_{\mathbf{n},x}}{k_{\mathbf{n}}}\right)}^2\right]+s_y^2\left[1-{\left(\frac{k_{\mathbf{n},y}}{k_{\mathbf{n}}}\right)}^2\right]+s_z^2\left[1-{\left(\frac{k_{\mathbf{n},z}}{k_{\mathbf{n}}}\right)}^2\right]\right\}{\sigma }_{\mathbf{n}}^2\\ -2{\displaystyle \sum }_{\mathbf{n}}\left[s_x{s}_y\frac{k_{\mathbf{n},x}{k}_{\mathbf{n},y}}{k_{\mathbf{n}}^2}+s_x{s}_z\frac{k_{\mathbf{n},x}{k}_{\mathbf{n},z}}{k_{\mathbf{n}}^2}+s_y{s}_z\frac{k_{\mathbf{n},y}{k}_{\mathbf{n},z}}{k_{\mathbf{n}}^2}\right]{\sigma }_{\mathbf{n}}^2.\end{array}(25)$

Although the “cross terms” involving $s_x{s}_y$, $s_x{s}_z$, $s_y{s}_z$ still remain, upon summing over $\mathbf{n}$ in the last three terms, we obtain

${\displaystyle \sum }_{\mathbf{n}}\frac{k_{\mathbf{n},i}{k}_{\mathbf{n},j}}{k_{\mathbf{n}}^2}{\sigma }_{\mathbf{n}}^2=0\text{ for }i\ne j,(26)$

since $n_x$, $n_y$, $n_z$ each vary as integers symmetrically from $-\infty$ to $+\infty$, where ${\mathbf{k}}_{\mathbf{n}}$ is given in Eq. 3.

Consequently, from Eqs. 21, 25 and 26:

$P\left(\mathbf{E}\text{ at }\mathbf{x},t\right)=I_1{I}_2{I}_3,(27)$

where

$I_i\equiv \frac{1}{2\pi }\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}{s}_i\text{exp}\left(i{s}_i{E}_i-\frac{1}{2}\frac{s_i^2}{\left(L_x{L}_y{L}_z\right)}{\displaystyle \sum }_{\mathbf{n}}\left[1-{\left(\frac{k_{\mathbf{n},i}}{k_{\mathbf{n}}}\right)}^2\right]{\sigma }_{\mathbf{n}}^2\right).(28)$

Simplifying notation, let

${\alpha }_i\equiv \frac{1}{\left(L_x{L}_y{L}_z\right)}{\displaystyle \sum }_{\mathbf{n}}\left[1-{\left(\frac{k_{\mathbf{n},i}}{k_{\mathbf{n}}}\right)}^2\right]{\sigma }_{\mathbf{n}}^2.(29)$

By then completing the square in Eq. 28 and carrying out the integral, yields

$\begin{array}{l}{I}_i=\frac{1}{2\pi }\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_i\text{exp}\left(-s_i^2{\alpha }_i\frac{1}{2}+i{s}_i{E}_i\right)\\ =\frac{1}{2\pi }\text{exp}\left(-\frac{E_i^2}{2{\alpha }_i}\right)\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_i\text{exp}\left[-\frac{{\alpha }_i}{2}{\left(s_i-\frac{i{E}_i}{{\alpha }_i}\right)}^2\right]\\ =\frac{1}{2\pi }\text{exp}\left(-\frac{E_i^2}{2{\alpha }_i}\right){\left(\frac{2\pi }{{\alpha }_i}\right)}^{1/2},\end{array}(30)$

resulting in

$P\left(\mathbf{E}\text{ at }\mathbf{x},t\right)=I_1{I}_2{I}_3=\frac{1}{{\left(2\pi \right)}^{\frac{3}{2}}}\frac{1}{{\left({\alpha }_1{\alpha }_2{\alpha }_3\right)}^{1/2}}\text{exp}\left(-\frac{E_x^2}{2{\alpha }_1}-\frac{E_y^2}{2{\alpha }_2}-\frac{E_z^2}{2{\alpha }_3}\right).(31)$

More insight into Eq. 31 can be gained by relating ${\alpha }_i$ in Eq. 29, to $\langle {\mathbf{E}}_{\text{ZPP},i}^2\left(\mathbf{x},t\right)\rangle$, using Eqs. 6,7, and 23:

$\begin{array}{c}\langle {\mathbf{E}}_{\text{ZPP},i}^2\left(\mathbf{x},t\right)\rangle \\ =\langle {\left\{\frac{1}{{\left(L_x{L}_y{L}_z\right)}^{1/2}}\stackrel{\infty }{{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }}{\displaystyle \sum }_{\lambda =\mathrm{1,2}}{\left({\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\right)}_i\times \left[A_{{\mathbf{k}}_{\mathbf{n}},\lambda }\text{cos}\left({\mathbf{k}}_{\mathbf{n}}\cdot \mathbf{x}-{\omega }_{\mathbf{n}}t\right)+B_{{\mathbf{k}}_{\mathbf{n}},\lambda }\text{sin}\left({\mathbf{k}}_{\mathbf{n}}\cdot \mathbf{x}-{\omega }_{\mathbf{n}}t\right)\right]\right\}}^2\rangle \\ =\frac{1}{\left(L_x{L}_y{L}_z\right)}\stackrel{\infty }{{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }}{\displaystyle \sum }_{\lambda =\mathrm{1,2}}{\left[{\left({\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\right)}_i\right]}^2{\sigma }_{\mathbf{n}}^2\left[{\text{cos}}^2\left({\mathbf{k}}_{\mathbf{n}}\cdot \mathbf{x}-{\omega }_{\mathbf{n}}t\right)+{\text{sin}}^2\left({\mathbf{k}}_{\mathbf{n}}\cdot \mathbf{x}-{\omega }_{\mathbf{n}}t\right)\right]\\ =\frac{1}{\left(L_x{L}_y{L}_z\right)}\stackrel{\infty }{{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }}\left[1-{\left(\frac{k_{\mathbf{n},i}}{k_{\mathbf{n}}}\right)}^2\right]{\sigma }_{\mathbf{n}}^2 .\end{array}(32)$

Combining Eqs. 29, 31 and 32 and noting

${\alpha }_i=\langle {\mathbf{E}}_{\text{ZPP},i}^2\left(\mathbf{x},t\right)\rangle ,(33)$

we obtain:

$P\left(\mathbf{E}\text{ at }\mathbf{x},t\right)=\frac{\text{exp}\left[-\frac{E_x^2}{2\langle E_{\text{ZPP},x}^2\rangle }\right]}{\sqrt{2\pi \langle E_{\text{ZPP},x}^2\rangle }}\frac{\text{exp}\left[-\frac{E_y^2}{2\langle E_{\text{ZPP},y}^2\rangle }\right]}{\sqrt{2\pi \langle E_{\text{ZPP},y}^2\rangle }}\frac{\text{exp}\left[-\frac{E_z^2}{2\langle E_{\text{ZPP},z}^2\rangle }\right]}{\sqrt{2\pi \langle E_{\text{ZPP},z}^2\rangle }}.(34)$

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

$\begin{array}{l}\langle {\mathbf{E}}_{\text{ZPP},x}^2\rangle =\langle {\mathbf{E}}_{\text{ZPP},y}^2\rangle =\langle {\mathbf{E}}_{\text{ZPP},z}^2\rangle =\langle {\mathbf{E}}_{\text{ZPP},i}^2\rangle \\ =\frac{1}{\left(L_x{L}_y{L}_z\right)}\stackrel{\infty }{{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }}\left\{1-{\left[\frac{{\left({\mathbf{k}}_{\mathbf{n}}\right)}_i}{k_{\mathbf{n}}}\right]}^2\right\}{\sigma }_{\mathbf{n}}^2=\frac{2}{3\left(L_x{L}_y{L}_z\right)}\stackrel{\infty }{{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }}{\sigma }_{\mathbf{n}}^2\end{array}(35)$

and

$P_{L_x=L_y=L_z}\left(\mathbf{E}\text{ at }\mathbf{x},t\right)=\frac{1}{{\left(2\pi \langle E_i^2\rangle \right)}^{3/2}}\text{exp}\left[-\frac{\left(E_x^2+E_y^2+E_z^2\right)}{2\langle E_i^2\rangle }\right].(36)$

Thus, the probability density for the radiation electric field $\mathbf{E}$ at position $\mathbf{x}$ and time t, from either Eqs. 34 or 36, equals the product of three Gaussian functions. Moreover, the probability density of $\mathbf{E}$ is independent of position $\mathbf{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:

$P\left(\mathbf{B}\text{ at }\mathbf{x},t\right)=\frac{\text{exp}\left[-\frac{B_x^2}{2\langle B_{\text{ZPP},x}^2\rangle }\right]}{\sqrt{2\pi \langle B_{\text{ZPP},x}^2\rangle }}\frac{\text{exp}\left[-\frac{B_y^2}{2\langle B_{\text{ZPP},y}^2\rangle }\right]}{\sqrt{2\pi \langle B_{\text{ZPP},y}^2\rangle }}\frac{\text{exp}\left[-\frac{B_z^2}{2\langle B_{\text{ZPP},z}^2\rangle }\right]}{\sqrt{2\pi \langle B_{\text{ZPP},z}^2\rangle }}.(37)$

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, $\langle {\mathbf{E}}_{\text{ZPP},i}^2\left(\mathbf{x},t\right)\rangle$ 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 $\left|\text{Δ}\mathbf{x}\right|=c\left|\text{Δ}t\right|$.

2.2 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

$\begin{array}{l}P\left({\mathbf{E}}_1\text{ at }{\mathbf{x}}_1,t_1;{\mathbf{E}}_2\text{ at }{\mathbf{x}}_2,t_2\right)\\ ={{\displaystyle \int }}^{\text{}}\text{d}{A}_1\cdots {{\displaystyle \int }}^{\text{}}\text{d}{A}_N\cdots {{\displaystyle \int }}^{\text{}}\text{d}{B}_1\cdots {{\displaystyle \int }}^{\text{}}\text{d}{B}_N\cdots P\left(A_1,\dots ,A_N,\dots ,B_N,\dots \right)\end{array}(38)$

$\times {\delta }^3\left[{\mathbf{E}}_1-{\mathbf{E}}_{\text{ZPP}}\left({\mathbf{x}}_1,t_1\right)\right]{\delta }^3\left[{\mathbf{E}}_2-{\mathbf{E}}_{\text{ZPP}}\left({\mathbf{x}}_2,t_2\right)\right],(39)$

where the semicolon in the first line is intended as a shortened meaning for the logical “AND” symbol of ${{\displaystyle \cap }}^{\text{}}.$

Again, we make use of the random variable independence of the $A^{\prime }s$ and $B^{\prime }s$ for a normal thermodynamic radiation situation, and impose the distribution Eq. 8, plus use Eq. 12, to obtain:

$\begin{array}{l}P\left({\mathbf{E}}_1\text{ at }{\mathbf{x}}_1,t_1;{\mathbf{E}}_2\text{ at }{\mathbf{x}}_2,t_2\right)\\ ={{\displaystyle \int }}^{\text{}}\text{d}{A}_1\cdots {{\displaystyle \int }}^{\text{}}\text{d}{A}_N\cdots \frac{1}{\sqrt{2\pi {\sigma }_{{\mathbf{n}}_1}^2}}\text{exp}\left[-\frac{{\left(A_1\right)}^2}{2{\sigma }_{{\mathbf{n}}_1}^2}\right]\cdots \frac{1}{\sqrt{2\pi {\sigma }_{{\mathbf{n}}_N}^2}}\text{exp}\left[-\frac{{\left(A_N\right)}^2}{2{\sigma }_{{\mathbf{n}}_N}^2}\right]\cdots \\ \times {{\displaystyle \int }}^{\text{}}\text{d}{B}_1\cdots {{\displaystyle \int }}^{\text{}}\text{d}{B}_N\cdots \frac{1}{\sqrt{2\pi {\sigma }_{{\mathbf{n}}_1}^2}}\text{exp}\left[-\frac{{\left(B_1\right)}^2}{2{\sigma }_{{\mathbf{n}}_1}^2}\right]\cdots \frac{1}{\sqrt{2\pi {\sigma }_{{\mathbf{n}}_N}^2}}\text{exp}\left[-\frac{{\left(B_N\right)}^2}{2{\sigma }_{{\mathbf{n}}_N}^2}\right]\cdots \\ \times \frac{1}{2\pi }\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_{1x}{e}^{i{s}_{1x}\left(E_{1x}-E_{1x,\text{ZPP}}\right)}\frac{1}{2\pi }\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_{1y}{e}^{i{s}_{1y}\left(E_{1y}-E_{1y,\text{ZPP}}\right)}\frac{1}{2\pi }\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_{1z}{e}^{i{s}_{1z}\left(E_{1z}-E_{1z,\text{ZPP}}\right)}\\ \times \frac{1}{2\pi }\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_{2x}{e}^{i{s}_{2x}\left(E_{2x}-E_{2x,\text{ZPP}}\right)}\frac{1}{2\pi }\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_{2y}{e}^{i{s}_{2y}\left(E_{2y}-E_{2y,\text{ZPP}}\right)}\frac{1}{2\pi }\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_{2z}{e}^{i{s}_{2z}\left(E_{2z}-E_{2z,\text{ZPP}}\right)},\end{array}(40)$

where the positions ${\mathbf{x}}_1$ and ${\mathbf{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,\text{ZPP}}$, $E_{1y,\text{ZPP}}$, and $E_{1z,\text{ZPP}}$ refer to ${\mathbf{E}}_{\text{ZPP}}\left({\mathbf{x}}_1,t_1\right)$ as in Eq. 15, and similarly for ${\mathbf{E}}_{\text{ZPP}}\left({\mathbf{x}}_2,t_2\right)$. Again abbreviating expressions, as in 16 and 17, with $a=\mathrm{1,2}$, as below:

$\begin{array}{c}{\mathbf{E}}_{\text{ZPP}}\left({\mathbf{x}}_a,t_a\right)\\ =\frac{1}{{\left(L_x{L}_y{L}_z\right)}^{1/2}}\stackrel{\infty }{{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }}{\displaystyle \sum }_{\lambda =\mathrm{1,2}}{\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\left[A_{{\mathbf{k}}_{\mathbf{n}},\lambda }\text{cos}\left({\mathbf{k}}_{\mathbf{n}}\cdot {\mathbf{x}}_{\mathbf{a}}-{\mathbf{\omega }}_{\mathbf{n}}{t}_a\right)+B_{{\mathbf{k}}_{\mathbf{n}},\lambda }\text{sin}\left({\mathbf{k}}_{\mathbf{n}}\cdot {\mathbf{x}}_{\mathbf{a}}-{\omega }_{\mathbf{n}}{t}_a\right)\right]\\ ={\displaystyle \sum }_q{A}_q{\mathbf{E}}_{a,cq}+{\displaystyle \sum }_q{B}_q{\mathbf{E}}_{a,sq}.\end{array}(41)$

Collecting the $A_q$-related terms in Eq. 40 as in the following manner, then later doing similarly for the $B_q$ terms:

$\begin{array}{l}-\frac{{\left(A_q\right)}^2}{2{\sigma }_q^2}-i{s}_{1x}{A}_q{E}_{1,cq,x}-i{s}_{1y}{A}_q{E}_{1,cq,y}-i{s}_{1z}{A}_q{E}_{1,cq,z}-i{s}_{2x}{A}_q{E}_{2,cq,x}-i{s}_{2y}{A}_q{E}_{2,cq,y}-i{s}_{2z}{A}_q{E}_{2,cq,z}=-\frac{1}{2{\sigma }_q^2}{\left[A_q+i\left(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}\right){\sigma }_q^2\right]}^2-{\left(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}\right)}^2\frac{{\sigma }_q^2}{2}.\end{array}(42)$

Now integrating over each $A_q$ term, followed later by integrating over the related $B_q$ term expression, results in:

$\begin{array}{l}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{A}_q\frac{1}{\sqrt{2\pi {\sigma }_q^2}}\text{exp}\left[-\frac{{\left(A_q\right)}^2}{2{\sigma }_q^2}\right]e^{-i\left(s_{1x}{A}_q{E}_{1,cq,x}+s_{1y}{A}_q{E}_{1,cq,y}+s_{1z}{A}_q{E}_{1,cq,z}+s_{2x}{A}_q{E}_{2,cq,x}+s_{2y}{A}_q{E}_{2,cq,y}+s_{2z}{A}_q{E}_{2,cq,z}\right)}\\ =\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{A}_q\frac{1}{\sqrt{2\pi {\sigma }_q^2}}\text{exp}\left[-\frac{1}{2{\sigma }_q^2}{\left[A_q+i\left(\begin{array}{c}{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}\end{array}\right){\sigma }_q^2\right]}^2\right]\\ \times \exp \left[-{\left(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}\right)}^2\frac{{\sigma }_q^2}{2}\right]\\ =\text{exp}\left[-{\left(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}\right)}^2\frac{{\sigma }_q^2}{2}\right].\end{array}(43)$

Integrating over each $A_q$ and $B_q$ results in:

$\begin{array}{l}P\left({\mathbf{E}}_1\text{ at }{\mathbf{x}}_1,t_1;{\mathbf{E}}_2\text{ at }{\mathbf{x}}_2,t_2\right)\\ =\frac{1}{{\left(2\pi \right)}^6}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_{1x}{e}^{i{s}_{1x}{E}_{1x}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_{1y}{e}^{i{s}_{1y}{E}_{1y}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_{1z}{e}^{i{s}_{1z}{E}_{1z}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_{2x}{e}^{i{s}_{2x}{E}_{2x}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_{2y}{e}^{i{s}_{2y}{E}_{2y}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_{2z}{e}^{i{s}_{2z}{E}_{2z}}\\ \times \exp \left[-{\displaystyle \sum }_q{\left(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}\right)}^2\frac{{\sigma }_q^2}{2}\right]\\ \times \exp \left[-{\displaystyle \sum }_q{\left(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}\right)}^2\frac{{\sigma }_q^2}{2}\right]\end{array}(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

${\text{cos}}^2\left({\mathbf{k}}_{\mathbf{n}}\cdot {\mathbf{x}}_a-{\omega }_{\mathbf{n}}{t}_a\right)+{\text{sin}}^2\left({\mathbf{k}}_{\mathbf{n}}\cdot {\mathbf{x}}_a-{\omega }_{\mathbf{n}}{t}_a\right)=1,(45)$

for either $a=1$ or 2, then:

$\begin{array}{l}{\left(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}\right)}^2\\ +{\left(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}\right)}^2\\ =\frac{1}{\left(L_x{L}_y{L}_z\right)}\left(s_{1x}^2{\epsilon }_{q,x}^2+s_{1y}^2{\epsilon }_{q,y}^2+s_{1z}^2{\epsilon }_{q,z}^2+s_{2x}^2{\epsilon }_{q,x}^2+s_{2y}^2{\epsilon }_{q,y}^2+s_{2z}^2{\epsilon }_{q,z}^2\right)\\ +\frac{2}{\left(L_x{L}_y{L}_z\right)}\left(s_{1x}{\epsilon }_{q,x}{s}_{1y}{\epsilon }_{q,y}+s_{1x}{\epsilon }_{q,x}{s}_{1z}{\epsilon }_{q,z}+s_{1y}{\epsilon }_{q,y}{s}_{1z}{\epsilon }_{q,z}\right)\\ +\frac{2}{\left(L_x{L}_y{L}_z\right)}\left(s_{2x}{\epsilon }_{q,x}{s}_{2y}{\epsilon }_{q,y}+s_{2x}{\epsilon }_{q,x}{s}_{2z}{\epsilon }_{q,z}+s_{2y}{\epsilon }_{q,y}{s}_{2z}{\epsilon }_{q,z}\right)\\ +\frac{2s_{1x}{\epsilon }_{q,x}}{\left(L_x{L}_y{L}_z\right)}\left(s_{2x}{\epsilon }_{q,x}+s_{2y}{\epsilon }_{q,y}+s_{2z}{\epsilon }_{q,z}\right)\left(C_1{C}_2+S_1{S}_2\right)\\ +\frac{2s_{1y}{\epsilon }_{q,y}}{\left(L_x{L}_y{L}_z\right)}\left(s_{2x}{\epsilon }_{q,x}+s_{2y}{\epsilon }_{q,y}+s_{2z}{\epsilon }_{q,z}\right)\left(C_1{C}_2+S_1{S}_2\right)\\ +\frac{2s_{1z}{\epsilon }_{q,z}}{\left(L_x{L}_y{L}_z\right)}\left(s_{2x}{\epsilon }_{q,x}+s_{2y}{\epsilon }_{q,y}+s_{2z}{\epsilon }_{q,z}\right)\left(C_1{C}_2+S_1{S}_2\right)\end{array}(46)$

The meaning of the abbreviated terms $\left(C_1{C}_2+S_1{S}_2\right)$ is:

$\begin{array}{l}\left(C_1{C}_2+S_1{S}_2\right)=\text{cos}\left({\mathbf{k}}_{\mathbf{n}}\cdot {\mathbf{x}}_1-{\omega }_{\mathbf{n}}{t}_1\right)\text{cos}\left({\mathbf{k}}_{\mathbf{n}}\cdot {\mathbf{x}}_2-{\omega }_{\mathbf{n}}{t}_2\right)+\text{sin}\left({\mathbf{k}}_{\mathbf{n}}\cdot {\mathbf{x}}_1-{\omega }_{\mathbf{n}}{t}_1\right)\text{sin}\left({\mathbf{k}}_{\mathbf{n}}\cdot {\mathbf{x}}_2-{\omega }_{\mathbf{n}}{t}_2\right)\\ =\text{cos}\left[{\mathbf{k}}_{\mathbf{n}}\cdot \left({\mathbf{x}}_1-{\mathbf{x}}_2\right)-{\omega }_{\mathbf{n}}\left(t_1-t_2\right)\right]\equiv C_{12,\mathbf{n}}.\end{array}(47)$

Again summing over the polarization indices of $\lambda =\mathrm{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 ${\epsilon }_{q,i}{\epsilon }_{q,j}$ for $i\ne j$, will drop out due to the sum in Eq. 26, results in:

$\begin{array}{c}{\displaystyle \sum }_q{\left(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}\right)}^2\\ +{\displaystyle \sum }_q{\left(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}\right)}^2\\ =\frac{1}{\left(L_x{L}_y{L}_z\right)}\stackrel{\infty }{{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }}\left[\left(s_{1x}^2+s_{2x}^2\right)\left(1-\frac{{\mathbf{k}}_{\mathbf{n},x}^{\mathbf{2}}}{k_{\mathbf{n}}^{\mathbf{2}}}\right)+\left(s_{1y}^2+s_{2y}^2\right)\left(1-\frac{{\mathbf{k}}_{\mathbf{n},y}^{\mathbf{2}}}{k_{\mathbf{n}}^{\mathbf{2}}}\right)+\left(s_{1z}^2+s_{2z}^2\right)\left(1-\frac{{\mathbf{k}}_{\mathbf{n},z}^{\mathbf{2}}}{k_{\mathbf{n}}^{\mathbf{2}}}\right)\right]\\ +\frac{2}{\left(L_x{L}_y{L}_z\right)}\stackrel{\infty }{{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }}\left[s_{1x}{s}_{2x}\left(1-\frac{{\mathbf{k}}_{\mathbf{n},x}^{\mathbf{2}}}{k_{\mathbf{n}}^{\mathbf{2}}}\right)+s_{1y}{s}_{2y}\left(1-\frac{{\mathbf{k}}_{\mathbf{n},y}^{\mathbf{2}}}{k_{\mathbf{n}}^{\mathbf{2}}}\right)+s_{1z}{s}_{2z}\left(1-\frac{{\mathbf{k}}_{\mathbf{n},z}^{\mathbf{2}}}{k_{\mathbf{n}}^{\mathbf{2}}}\right)\right]C_{12,\mathbf{n}} .\end{array}(48)$

Hence:

$\begin{array}{c}P\left({\mathbf{E}}_1\text{ at }{\mathbf{x}}_1,t_1;{\mathbf{E}}_2\text{ at }{\mathbf{x}}_2,t_2\right)\\ =\frac{1}{{\left(2\pi \right)}^6}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{s}_{1x}{e}^{i{s}_{1x}{E}_{1x}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{s}_{1y}{e}^{i{s}_{1y}{E}_{1y}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{s}_{1z}{e}^{i{s}_{1z}{E}_{1z}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{s}_{2x}{e}^{i{s}_{2x}{E}_{2x}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{s}_{2y}{e}^{i{s}_{2y}{E}_{2y}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{s}_{2z}{e}^{i{s}_{2z}{E}_{2z}}\\ \times \exp \left[-\frac{1}{\left(L_x{L}_y{L}_z\right)}{\displaystyle \sum _{n_x,n_y,n_z=-\infty }^{\infty }\left\{\begin{array}{c}\left(s_{1x}^2+2s_{1x}{s}_{2x}{C}_{12,\tilde{\mathbf{n}}}+s_{\mathbf{2}x}^2\right)\left(1-\frac{{\mathbf{k}}_{\mathbf{n},x}^{\mathbf{2}}}{k_{\mathbf{n}}^{\mathbf{2}}}\right)\\ +\left(s_{1y}^2+2s_{1y}{s}_{2y}{C}_{12,\tilde{\mathbf{n}}}+s_{2y}^2\right)\left(1-\frac{{\mathbf{k}}_{\mathbf{n},y}^{\mathbf{2}}}{k_{\mathbf{n}}^{\mathbf{2}}}\right)\\ +\left(s_{1z}^2+2s_{1z}{s}_{2z}{C}_{12,\tilde{\mathbf{n}}}+s_{2z}^2\right)\left(1-\frac{{\mathbf{k}}_{\mathbf{n},z}^{\mathbf{2}}}{k_{\mathbf{n}}^{\mathbf{2}}}\right)\end{array}\right\}}\frac{{\sigma }_q^2}{2}\right]\\ =I_{1x2x}{I}_{1y2y}{I}_{1z2z}\end{array}(49)$

where

$I_{1x2x}\equiv \frac{1}{{\left(2\pi \right)}^2}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{s}_{1x}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{s}_{2x}\text{exp}\left[\begin{array}{c}i{s}_{1x}{E}_{1x}+i{s}_{2x}{E}_{2x}\\ -\frac{1}{\left(L_x{L}_y{L}_z\right)}{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }^{\infty }\left(s_{1x}^2+2s_{1x}{s}_{2x}{C}_{12,\mathbf{n}}+s_{2x}^2\right)\left(1-\frac{{\mathbf{k}}_{\mathbf{n},x}^{\mathbf{2}}}{k_{\mathbf{n}}^{\mathbf{2}}}\right)\frac{{\sigma }_{\mathbf{n}}^{\mathbf{2}}}{2}\end{array}\right]},(50)$

$I_{1y2y}\equiv \frac{1}{{\left(2\pi \right)}^2}\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}d{s}_{1y}\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}d{s}_{2y}\text{exp}\left[\begin{array}{c}i{s}_{1y}{E}_{1y}+i{s}_{2y}{E}_{2y}\\ -\frac{1}{\left(L_x{L}_y{L}_z\right)}{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }^{\infty }\left(s_{1y}^2+2s_{1y}{s}_{2y}{C}_{12,\mathbf{n}}+s_{2y}^2\right)\left(1-\frac{\frac{{\mathbf{k}}_{\mathbf{n},y}^{\mathbf{2}}}{k_{\mathbf{n}}^2}}{k_{\mathbf{n}}^{\mathbf{2}}}\right)\frac{{\sigma }_{\mathbf{n}}^2}{2}\end{array}\right] ,(51)$

$I_{1z2z}\equiv \frac{1}{{\left(2\pi \right)}^2}\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}d{s}_{1z}\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}d{s}_{2z}\text{exp}\left[\begin{array}{c}i{s}_{1z}{E}_{1z}+i{s}_{2z}{E}_{2z}\\ -\frac{1}{\left(L_x{L}_y{L}_z\right)}{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }^{\infty }\left(s_{1z}^2+2s_{1z}{s}_{2z}{C}_{12,\mathbf{n}}+s_{2z}^2\right)\left(1-\frac{{\mathbf{k}}_{\mathbf{n},z}^{\mathbf{2}}}{k_{\mathbf{n}}^{\mathbf{2}}}\right)\frac{{\sigma }_{\mathbf{n}}^{\mathbf{2}}}{2}\end{array}\right] .(52)$

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}$.

$I_{1x2x}=\frac{1}{{\left(2\pi \right)}^2}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d}{s}_{1x}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d}{s}_{2x}\exp \left[\begin{array}{c}-s_{1x}^2\frac{1}{\left(L_x{L}_y{L}_z\right)}{\displaystyle \sum _{n_x,n_y,n_z=-\infty }^{\infty }{K}_{\mathbf{n},x}^{\mathbf{2}}}\frac{{\sigma }_{\mathbf{n}}^{\mathbf{2}}}{2}\\ +s_{1x}\left(i{E}_{1x}-2s_{2x}\frac{1}{\left(L_x{L}_y{L}_z\right)}{\displaystyle \sum _{n_x,n_y,n_z=-\infty }^{\infty }{C}_{12,\mathbf{n}}}{K}_{\mathbf{n}}^{\mathbf{2}}\frac{{\sigma }_{\mathbf{n}}^{\mathbf{2}}}{2}\right)\\ +i{s}_{2x}{E}_{2x}-s_{2x}^2\frac{1}{\left(L_x{L}_y{L}_z\right)}{\displaystyle \sum _{n_x,n_y,n_z=-\infty }^{\infty }{K}_{\mathbf{n},\mathbf{x}}^{\mathbf{2}}}\frac{{\sigma }_{\mathbf{n}}^{\mathbf{2}}}{2}\end{array}\right](53)$

where

$K_{\mathbf{n},x}^2\equiv \left(1-\frac{{\mathbf{k}}_{\mathbf{n},x}^{\mathbf{2}}}{k_{\mathbf{n}}^{\mathbf{2}}}\right).(54)$

Let

$A_x\equiv \frac{1}{\left(L_x{L}_y{L}_z\right)}{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }^{\infty }{K}_{\mathbf{n},x}^2{\sigma }_{\mathbf{n}}^2(55)$

and

$B_x\equiv \left(i{E}_{1x}-\frac{s_{2x}}{\left(L_x{L}_y{L}_z\right)}{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }^{\infty }{C}_{12,\mathbf{n}}{K}_{\mathbf{n},x}^2{\sigma }_{\mathbf{n}}^2\right),(56)$

where $C_{12,\mathbf{n}}$ was defined in Eq. 47. Then:

$\begin{array}{c}{I}_{1x2x}=\frac{1}{{\left(2\pi \right)}^2}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{s}_{1x}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{s}_{2x}\text{exp}\left[-\frac{1}{2}{s}_{1x}^2{A}_x+s_{1x}{B}_x+i{s}_{2x}{E}_{2x}-\frac{1}{2}{s}_{2x}^2{A}_x\right]}\\ =\frac{1}{{\left(2\pi \right)}^2}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{s}_{2x}\text{exp}\left(i{s}_{2x}{E}_{2x}-\frac{1}{2}{s}_{2x}^2{A}_x\right)}\times {\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{s}_{1x}\text{exp}\left(-\frac{1}{2}{s}_{1x}^2{A}_x+s_{1x}{B}_x\right)}\end{array}(57)$

Completing the square with

$-\mathfrak{u}{u}^2+\mathfrak{B}u=-\left(\mathfrak{u}{u}^2-\mathfrak{B}u\right)=-{\left(\sqrt{\mathfrak{u}}u-\frac{\mathfrak{B}}{2\sqrt{\mathfrak{u}}}\right)}^2+\frac{{\mathfrak{B}}^2}{4\mathfrak{u}} ,(58)$

results in:

$\begin{array}{l}{I}_{1x2x}=\frac{1}{{\left(2\pi \right)}^2}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_{2x}\text{exp}\left(i{s}_{2x}{E}_{2x}-\frac{1}{2}{s}_{2x}^2{A}_x\right)\\ \times \left\{\text{exp}\left(\frac{B_x^2}{2A_x}\right)\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_{1x}\text{exp}\left[-{\left(\sqrt{\frac{A_x}{2}}{s}_{1x}-\frac{B_x}{2\sqrt{\frac{A_x}{2}}}\right)}^2\right]\right\}\\ =\frac{1}{{\left(2\pi \right)}^2}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{s}_{2x}\text{exp}\left(i{s}_{2x}{E}_{2x}-\frac{1}{2}{s}_{2x}^2{A}_x\right)\left[\text{exp}\left(\frac{B_x^2}{2A_x}\right)\frac{\sqrt{\pi }}{{\left(\frac{1}{2}{A}_x\right)}^{1/2}}\right].\end{array}(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

$C_x\equiv \frac{1}{\left(L_x{L}_y{L}_z\right)}{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }^{\infty }{C}_{12,\mathbf{n}}{K}_{\mathbf{n},x}^2{\sigma }_{\mathbf{n}}^2,(60)$

so that

$B_x\equiv i{E}_{1x}-s_{2x}{C}_x.(61)$

Then:

$\begin{array}{c}{I}_{1x2x}=\frac{\sqrt{2\pi }}{{\left(2\pi \right)}^2}\frac{1}{A_x^{1/2}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{s}_{2x}\text{exp}\left(i{s}_{2x}{E}_{2x}-\frac{1}{2}{s}_{2x}^2{A}_x\right)\times \text{exp}\left[\frac{1}{2A_x}{\left(i{E}_{1x}-s_{2x}{C}_x\right)}^2\right]}\\ =\frac{1}{{\left(2\pi \right)}^{3/2}{A}_x^{1/2}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{s}_{2x}\text{exp}\left(i{s}_{2x}{E}_{2x}-\frac{1}{2}{s}_{2x}^2{A}_x\right)\times \text{exp}\left[-\frac{E_{1x}^2}{2A_x}-\frac{s_{2x}i{E}_{1x}{C}_x}{A_x}+\frac{s_{2x}^2{C}_x^2}{2A_x}\right]}\\ =\frac{1}{{\left(2\pi \right)}^{3/2}{A}_x^{1/2}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{s}_{2x}\text{exp}\left\{-s_{2x}^2\left[\frac{1}{2}{A}_x-\frac{C_x^2}{2A_x}\right]+s_{2x}\left(i{E}_{2x}-\frac{i{E}_{1x}{C}_x}{A_x}\right)-\frac{E_{1x}^2}{2A_x}\right\}}.\end{array}$

Applying Eq. 58 again:

$\begin{array}{c}{I}_{1x2x}=\frac{\exp \left(-\frac{E_{1x}^2}{2A_x}\right)}{{\left(2\pi \right)}^{3/2}{A}_x^{1/2}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d}{s}_{2x}\times \exp \left(-{\left[s_{2x}{\left(\frac{A_x}{2}-\frac{C_x^2}{2A_x}\right)}^{1/2}-\frac{\left(i{E}_{2x}-\frac{i{E}_{1x}{C}_x}{A_x}\right)}{2{\left(\frac{A_x}{2}-\frac{C_x^2}{2A_x}\right)}^{1/2}}\right]}^2+\frac{{\left(i{E}_{2x}-\frac{i{E}_{1x}{C}_x}{A_x}\right)}^2}{4\left(\frac{A_x}{2}-\frac{C_x^2}{2A_x}\right)}\right)\\ =\frac{\text{exp}\left(-\frac{E_{1x}^2}{2A_x}\right)}{{\left(2\pi \right)}^{3/2}{A}_x^{1/2}}\frac{\sqrt{\pi }}{{\left[\frac{1}{2}{A}_x-\frac{C_x^2}{2A_x}\right]}^{1/2}}\text{exp}\left\{-\frac{{\left(E_{2x}-\frac{E_{1x}{C}_x}{A_x}\right)}^2}{4\left[\frac{1}{2}{A}_x-\frac{C_x^2}{2A_x}\right]}\right\}\\ =\frac{1}{2\pi A_x{\left(1-\frac{C_x^2}{A_x^2}\right)}^{1/2}}\text{exp}\left[\frac{-E_{1x}^2-E_{2x}^2+2E_{2x}{E}_{1x}\frac{C_x}{A_x}}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\end{array}(62)$

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, $\left({\mathbf{x}}_1-{\mathbf{x}}_2\right)$ and $\left(t_1-t_2\right)$, 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\left({\mathbf{E}}_1\text{ at }{\mathbf{x}}_1,t_1;{\mathbf{E}}_2\text{ at }{\mathbf{x}}_2,t_2\right)$ in Eq. 49 is just the product $I_{1x2x}\cdot I_{1y2y}\cdot 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},\text{ and }{E}_{2z}$. Again, these points will be made clearer in Section 2.4.

2.3 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

$\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}{E}_{2x}\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}{E}_{1x}{I}_{1x2x}$

should equal unity. Checking:

$\begin{array}{c}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}{I}_{1x2x}\\ =\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}\frac{1}{2\pi A_x{\left(1-\frac{C_x^2}{A_x^2}\right)}^{1/2}}\times \text{exp}\left\{\frac{-E_{1x}^2-E_{2x}^2+2E_{2x}{E}_{1x}\frac{C_x}{A_x}}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right\}\\ =\frac{1}{2\pi A_x{\left(1-\frac{C_x^2}{A_x^2}\right)}^{1/2}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}\text{exp}\left[-\frac{E_{2x}^2}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\times \stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}\text{exp}\left[\frac{-{\left(E_{1x}-E_{2x}\frac{C_x}{A_x}\right)}^2+E_{2x}^2\frac{C_x^2}{A_x^2}}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\\ =\frac{1}{2\pi A_x{\left(1-\frac{C_x^2}{A_x^2}\right)}^{1/2}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}\text{exp}\left[\frac{-E_{2x}^2}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\times \text{exp}\left[\frac{+E_{2x}^2\frac{C_x^2}{A_x^2}}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\\ \times \stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}\text{exp}\left[\frac{-{\left(E_{1x}-E_{2x}\frac{C_x}{A_x}\right)}^2}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\\ =\frac{1}{2\pi A_x{\left(1-\frac{C_x^2}{A_x^2}\right)}^{1/2}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}\text{exp}\left(-\frac{E_{2x}^2}{2A_x}\right)\left[\sqrt{\pi }\sqrt{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\\ =\frac{1}{\sqrt{2\pi A_x}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}\text{exp}\left(-\frac{E_{2x}^2}{2A_x}\right)=\frac{1}{\sqrt{2\pi A_x}}\sqrt{\pi 2A_x}=1 ,\end{array}(63)$

so this is fine.

Another check is whether

$\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}{E}_{2i}{I}_{1i2i}=P\left(E_{1i}\right),(64)$

and of course the opposite situation of $\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}{E}_{1i}{I}_{1i2i}=P\left(E_{2i}\right)$. In Eq. 64, we have already deduced from earlier work, including Eq. 34, that

$P\left(E_i\text{ at }{\mathbf{x}}_1,t_1\right)=\frac{\text{exp}\left[-\frac{E_i^2}{2\langle E_{\text{ZPP},i}^2\rangle }\right]}{\sqrt{2\pi \langle E_{\text{ZPP},i}^2\rangle }}.(65)$

As shown earlier, Eq. 65 turns out to be independent of ${\mathbf{x}}_1,t_1$, where $i=\mathrm{1,2,3}$ refers to $x,y,z$, respectively. Mathematically, this independence on ${\mathbf{x}}_1,t_1$ arises because $\langle E_{\text{ZPP},i}^2\rangle$ is independent of ${\mathbf{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=\mathrm{1,2},\text{ and }3$, to include all three $x,y,z$ cases, then:

$\langle E_{\text{ZPP},i}^2\rangle =\frac{1}{\left(L_x{L}_y{L}_z\right)}{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }^{\infty }\left(1-\frac{{\mathbf{k}}_{\mathbf{n},i}^2}{k_{\mathbf{n}}^2}\right){\sigma }_{\mathbf{n}}^2=A_i.(66)$

Returning to Eq. 64 and using 50:

$\begin{array}{c}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{E}_{2i}{I}_{1i2i}}\\ =\frac{\text{exp}\left[\frac{-E_{1i}^2}{2A_i\left(1-\frac{C_i^2}{A_i^2}\right)}\right]}{2\pi A_i{\left(1-\frac{C_i^2}{A_i^2}\right)}^{1/2}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{E}_{2i}\text{exp}\left[\frac{-E_{2i}^2+2E_{2i}{E}_{1i}\frac{C_i}{A_i}}{2A_i\left(1-\frac{C_i^2}{A_i^2}\right)}\right]}\\ =\frac{\exp \left[\frac{-E_{1i}^2}{2A_i\left(1-\frac{C_i^2}{A_i^2}\right)}\right]}{2\pi A_i{\left(1-\frac{C_i^2}{A_i^2}\right)}^{1/2}}\exp \left[\frac{E_{1i}^2\frac{C_i^2}{A_i^2}}{2A_i\left(1-\frac{C_i^2}{A_i^2}\right)}\right]{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{E}_{2i}\times \exp \left[-{\left(\frac{1}{\sqrt{2A_i\left(1-\frac{C_i^2}{A_i^2}\right)}}{E}_{2i}-\frac{E_{1i}\frac{C_i}{A_i}}{\sqrt{2A_i\left(1-\frac{C_i^2}{A_i^2}\right)}}\right)}^2\right]}\\ =\frac{1}{2\pi A_i{\left(1-\frac{C_i^2}{A_i^2}\right)}^{1/2}}\text{exp}\left[\frac{-E_{1i}^2}{2A_i}\right]\sqrt{\pi 2A_i\left(1-\frac{C_i^2}{A_i^2}\right)}=\frac{1}{\sqrt{2\pi A_i}}\text{exp}\left[\frac{-E_{1i}^2}{2A_i}\right]=\frac{\text{exp}\left[-\frac{E_{1i}^2}{2\langle E_{\text{ZPP},1i}^2\rangle }\right]}{\sqrt{2\pi \langle E_{\text{ZPP},1i}^2\rangle }},\end{array}$

using Eq. 66 at the end. By symmetry, $\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}{E}_{1i}{I}_{1i2i}=P\left(E_{2i}\right)$ then also holds.

Another obvious item to check is whether $\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}{E}_{2x}\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}{E}_{1x}{E}_{1x}{I}_{1x2x}=0$ (or vice versa$\langle E_{2x}\rangle =0$, by symmetry):

$\begin{array}{c}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}{E}_{1x}{I}_{1x2x}\\ =\frac{1}{2\pi A_x{\left(1-\frac{C_x^2}{A_x^2}\right)}^{1/2}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}\times \text{exp}\left(-\frac{E_{2x}^2}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right)\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}{E}_{1x}\text{exp}\left[\frac{-E_{1x}^2+2E_{2x}{E}_{1x}\frac{C_x}{A_x}}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\\ =\frac{1}{2\pi A_x{\left(1-\frac{C_x^2}{A_x^2}\right)}^{1/2}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}\text{exp}\left(-\frac{E_{2x}^2}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right)\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}{E}_{1x}\times \text{exp}\left[\frac{-{\left(E_{1x}-E_{2x}\frac{C_x}{A_x}\right)}^2+E_{2x}^2\frac{C_x^2}{A_x^2}}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\\ =\frac{1}{2\pi A_x{\left(1-\frac{C_x^2}{A_x^2}\right)}^{1/2}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}\times \text{exp}\left(-\frac{E_{2x}^2}{2A_x}\right)\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}{E}_{1x}\text{exp}\left[\frac{-{\left(E_{1x}-E_{2x}\frac{C_x}{A_x}\right)}^2}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right] .\end{array}(67)$

In the integral on the far right, bottom line, let $E_{1x}-E_{2x}\frac{C_x}{A_x}=u$ and $\text{d}u=\text{d}{E}_{1x}$:

$\begin{array}{c}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}{E}_{1x}{I}_{1x2x}\\ =\frac{1}{2\pi A_x{\left(1-\frac{C_x^2}{A_x^2}\right)}^{1/2}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}\text{exp}\left(-\frac{E_{2x}^2}{2A_x}\right)\times \stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}du\left(u+E_{2x}\frac{C_x}{A_x}\right)\text{exp}\left[\frac{-u^2}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\\ =\frac{1}{2\pi A_x{\left(1-\frac{C_x^2}{A_x^2}\right)}^{1/2}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}\times \text{exp}\left(-\frac{E_{2x}^2}{2A_x}\right)\left\{0+E_{2x}\frac{C_x}{A_x}{\left[\pi 2A_x\left(1-\frac{C_x^2}{A_x^2}\right)\right]}^{1/2}\right\}\\ =\frac{\frac{C_x}{A_x}}{{\left(2\pi A_x\right)}^{1/2}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}{E}_{2x}\text{exp}\left(-\frac{E_{2x}^2}{2A_x}\right)=0 ,\end{array}(68)$

since the integral is odd in $E_{2x}$.

Now trying a check that is more involved, we will compute a 2-point 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:

$\begin{array}{c}\langle E_{\text{ZPP},x}\left(\mathrm{0,0}\right)E_{\text{ZPP},x}(\widehat{\mathbf{y}}R,t)\rangle \\ ={\left(\frac{1}{{\left(L_x{L}_y{L}_z\right)}^{1/2}}\right)}^2\stackrel{\infty }{{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }}{\displaystyle \sum }_{\lambda =\mathrm{1,2}}{\left({\widehat{\mathit{\epsilon }}}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\right)}_x\stackrel{\infty }{{\displaystyle \sum }_{n_x\prime ,n_y\prime ,n_z\prime =-\infty }}{\displaystyle \sum }_{\lambda \prime =\mathrm{1,2}}{\left({\widehat{\mathit{\epsilon }}}_{{\mathbf{k}}_{\mathbf{n}\prime },\lambda \prime }\right)}_x\times \\ \times \langle \begin{array}{c}\left[A_{{\mathbf{k}}_{\mathbf{n}},\lambda }\text{cos}\left({\mathbf{k}}_{\mathbf{n}}\cdot 0-{\omega }_{\mathbf{n}}0\right)+B_{{\mathbf{k}}_{\mathbf{n}},\lambda }\text{sin\hspace{0.17em}}\left({\mathbf{k}}_{\mathbf{n}}\cdot 0-{\omega }_{\mathbf{n}}0\right)\right]\times \\ \times \left[A_{{\mathbf{k}}_{\mathbf{n}\prime },\lambda \prime }\text{cos\hspace{0.17em}}({\mathbf{k}}_{\mathbf{n}\prime }\cdot \widehat{\mathbf{y}}R-{\omega }_{\mathbf{n}\prime }t)+B_{{\mathbf{k}}_{\mathbf{n}\prime },\lambda \prime }\text{sin\hspace{0.17em}}\left({\mathbf{k}}_{\mathbf{n}\prime }\cdot \widehat{\mathbf{y}}R-{\omega }_{\mathbf{n}\prime }t\right)\right]\end{array}\rangle \\ =\frac{1}{\left(L_x{L}_y{L}_z\right)}{\displaystyle \sum _{\begin{array}{c}{n}_x,n_y,n_z=-\infty \\ \lambda =\mathrm{1,2}\end{array}}^{\infty }{\left({\widehat{\mathit{\epsilon }}}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\right)}_x}{\displaystyle \sum _{n_x^{\prime },n_y^{\prime },n_z^{\prime }=-\infty }^{\infty }{\displaystyle \sum _{{\lambda }^{\prime }=\mathrm{1,2}}{\left({\widehat{\mathit{\epsilon }}}_{{\mathbf{k}}_{{\mathbf{n}}^{\prime }},{\lambda }^{\prime }}\right)}_x}}\delta {\mathbf{n}}_{\mathbf{n}\prime }{\delta }_{\lambda ,{\lambda }^{\prime }}\mathrm{cos\hspace{0.17em}}\left({\mathbf{k}}_{{\mathbf{n}}^{\prime }}·\widehat{\mathbf{y}}R-{\omega }_{{\tilde{\mathbf{n}}}^{\prime }}t\right){\sigma }^2\left({\omega }_{\tilde{\mathbf{n}}},T\right)\\ =\frac{1}{\left(L_x{L}_y{L}_z\right)}\stackrel{\infty }{{\displaystyle \sum }_{\begin{array}{c}{n}_x,n_y,n_z=-\infty \\ \lambda =\mathrm{1,2}\end{array}}}{\left({\widehat{\mathit{\epsilon }}}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\right)}_x^2\text{cos}\left({\mathbf{k}}_{\mathbf{n}}\cdot \widehat{\mathbf{y}}R-{\omega }_{\mathbf{n}}t\right){\sigma }^2\left({\omega }_{\mathbf{n}},T\right)\\ =\frac{1}{\left(L_x{L}_y{L}_z\right)}\stackrel{\infty }{{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }}\left[1-{\left(\frac{{\mathbf{k}}_{\mathbf{n},x}}{k_{\mathbf{n}}}\right)}^2\right]\text{cos\hspace{0.17em}}\left({\mathbf{k}}_{\mathbf{n},y}R-{\omega }_{\mathbf{n}}t\right){\sigma }^2\left({\omega }_{\mathbf{n}},T\right)\end{array}(69)$

Making the change of discrete to continuous variables, with ${\mathbf{k}}_{\mathbf{n}}=\frac{2\pi n_x}{L_x}\widehat{\mathbf{x}}+\frac{2\pi n_y}{L_y}\widehat{\mathbf{y}}+\frac{2\pi n_z}{L_z}\widehat{\mathbf{z}}$, for large values of $L_x$, $L_y$, $L_z$, with

$\text{Δ}{n}_x\approx \frac{L_x}{2\pi }\text{d}{k}_x,\text{Δ}{n}_y\approx \frac{L_y}{2\pi }\text{d}{k}_y,\text{Δ}{n}_z\approx \frac{L_z}{2\pi }\text{d}{k}_z(70)$

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 ${\left[\sigma \left(\omega ,T\right)\right]}^2\to 2\pi \hslash \omega =2\pi \hslash kc$:

$\begin{array}{c}\langle E_{\text{ZP},x}\left(\mathrm{0,0}\right)E_{\text{ZP},x}\left(\widehat{\mathbf{y}}R,t\right)\rangle \\ \approx \frac{1}{\left(2\pi \right)}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{k}_x\frac{1}{\left(2\pi \right)}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{k}_y\frac{1}{\left(2\pi \right)}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{k}_z\left[1-{\left(\frac{{\mathbf{k}}_{\mathbf{n},x}}{k_{\mathbf{n}}}\right)}^2\right]\times \text{cos}\left(k_{y}R-\omega t\right)2\pi \hslash \omega }\\ =\frac{2\pi \hslash }{{\left(2\pi \right)}^3}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{k}_x}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{k}_y}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{k}_{z}kc\left[1-{\left(\frac{k_x}{k}\right)}^2\right]\times [\text{cos}\left(k_{y}R\right)\text{cos}\left(\omega t\right)+\text{sin}\left(k_{y}R\right)\text{sin}\left(\omega t\right)]}.\end{array}$

The second term of $\text{sin}\left(k_{y}R\right)\text{sin}\left(\omega t\right)$ makes the integrand odd in $k_y$. Hence:

$\langle E_{\text{ZP},x}\left(\mathrm{0,0}\right)E_{\text{ZP},x}\left(\widehat{\mathbf{y}}R,t\right)\rangle =\frac{\hslash c}{{\left(2\pi \right)}^2}\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}{k}_x\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}{k}_y\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}{k}_z\times K\left[1-{\left(\frac{k_x}{k}\right)}^2\right]\text{cos}\left(k_{y}R\right)\text{cos}\left(kct\right).$

Next, we will make $k_y$ be the axis where the polar angle is measured from, so that $k_y=k\text{cos}\theta$. Consequently, $k_x=k\text{sin}\theta \text{sin}\varphi$ and $k_z=k\text{sin}\theta \text{cos}\varphi$. Hence:

$\langle E_{\text{ZP},x}\left(\mathrm{0,0}\right)E_{\text{ZP},x}\left(\widehat{\mathbf{y}}R,t\right)\rangle =\frac{\hslash c}{{\left(2\pi \right)}^2}\underset{0}{\overset{\infty }{{\displaystyle \int }}}\text{d}k{k}^2\underset{0}{\overset{\pi }{{\displaystyle \int }}}\text{d}\theta \text{sin}\theta \times \underset{0}{\overset{2\pi }{{\displaystyle \int }}}\text{d}\varphi k\left[1-{\text{sin}}^2\theta {\text{sin}}^2\varphi \right]\text{cos}\left(Rk\text{cos}\theta \right)\text{cos}\left(kct\right).(71)$

Since

$\frac{1}{2\pi }\underset{0}{\overset{2\pi }{{\displaystyle \int }}}\text{d}\varphi {\text{sin}}^2\varphi =\frac{1}{2},$

then:

$\langle E_{\text{ZP},x}\left(\mathrm{0,0}\right)E_{\text{ZP},x}\left(\widehat{\mathbf{y}}R,t\right)\rangle =\frac{\hslash c\pi }{{\left(2\pi \right)}^2}\underset{0}{\overset{\infty }{{\displaystyle \int }}}\text{d}k{k}^3\text{cos}\left(kct\right)\underset{0}{\overset{\pi }{{\displaystyle \int }}}\text{d}\theta \text{sin}\theta \left[2-{\text{sin}}^2\theta \right]\text{cos}\left(Rk\text{cos}\theta \right).(72)$

Using

$\underset{0}{\overset{\pi }{{\displaystyle \int }}}\text{d}\theta \text{sin}\theta \left[2-{\text{sin}}^2\theta \right]\text{cos}\left(Rk\text{cos}\theta \right)=2\underset{0}{\overset{\pi }{{\displaystyle \int }}}\text{d}\theta \text{sin}\theta \text{cos}\left(Rk\text{cos}\theta \right)-\underset{0}{\overset{\pi }{{\displaystyle \int }}}\text{d}\theta {\text{sin}}^3\theta \text{cos}\left(Rk\text{cos}\theta \right),(73)$

with $u=Rk\text{cos}\theta$, $du=-Rk\text{sin}\theta d\theta$, the first term in Eq. 73 becomes:

$2\underset{0}{\overset{\pi }{{\displaystyle \int }}}\text{d}\theta \text{sin}\theta \text{cos}\left(Rk\text{cos}\theta \right)=2\underset{Rk}{\overset{-Rk}{{\displaystyle \int }}}\left(-\frac{\text{d}u}{Rk}\right)\text{cos}\left(u\right)=-\frac{2}{Rk}\text{sin}u\underset{Rk}{\overset{-Rk}{{\displaystyle |}}}=\frac{4}{Rk}\text{sin}\left(Rk\right).(74)$

The second term in Eq. 73 becomes:

$\begin{array}{c}-\stackrel{\pi }{\underset{0}{{\displaystyle \int }}}d\theta {\text{sin}}^3\theta \text{cos}\left(Rk\text{cos}\theta \right)\\ =-\stackrel{+Rk}{\underset{-Rk}{{\displaystyle \int }}}\frac{du}{Rk}\left[1-{\text{cos}}^2\theta \right]\text{cos}\left(u\right)=-\frac{1}{Rk}\stackrel{+Rk}{\underset{-Rk}{{\displaystyle \int }}}du\text{cos}\left(u\right)+\frac{1}{{\left(Rk\right)}^3}\stackrel{+Rk}{\underset{-Rk}{{\displaystyle \int }}}du\cdot u^2\text{cos}\left(u\right)\\ =-\frac{1}{Rk}\text{sin}\left(u\right){|}_{-Rk}^{Rk}+\frac{1}{{\left(Rk\right)}^3}\left[4Rk\text{cos}\left(Rk\right)-4\text{sin}\left(Rk\right)+2{\left(Rk\right)}^2\text{sin}\left(Rk\right)\right]\\ =+\frac{4}{{\left(Rk\right)}^2}\text{cos}\left(Rk\right)-\frac{4}{{\left(Rk\right)}^3}\text{sin}\left(Rk\right) .\end{array}(75)$

From Eqs. 73–75:

$\begin{array}{l}\stackrel{\pi }{\underset{0}{{\displaystyle \int }}}\text{d}\theta \text{sin}\theta \left[2-{\text{sin}}^2\theta \right]\text{cos}\left(Rk\text{cos}\theta \right)=2\stackrel{\pi }{\underset{0}{{\displaystyle \int }}}\text{d}\theta \text{sin}\theta \text{cos}\left(Rk\text{cos}\theta \right)-\stackrel{\pi }{\underset{0}{{\displaystyle \int }}}\text{d}\theta {\text{sin}}^3\theta \text{cos}\left(Rk\text{cos}\theta \right)\\ =\frac{4}{Rk}\text{sin}\left(Rk\right)+\frac{4}{{\left(Rk\right)}^2}\text{cos}\left(Rk\right)-\frac{4}{{\left(Rk\right)}^3}\text{sin}\left(Rk\right).\end{array}$

Consequently, Eq. 72 becomes:

$\begin{array}{l}\langle E_{\text{ZP},x}\left(\mathrm{0,0}\right)E_{\text{ZP},x}\left(\widehat{\mathbf{y}}R,t\right)\rangle \\ =\frac{\hslash c\pi }{{\left(2\pi \right)}^2}\stackrel{\infty }{\underset{0}{{\displaystyle \int }}}\text{d}k{k}^3\text{cos}\left(kct\right)\left[\frac{4}{Rk}\text{sin}\left(Rk\right)+\frac{4}{{\left(Rk\right)}^2}\text{cos}\left(Rk\right)-\frac{4}{{\left(Rk\right)}^3}\text{sin}\left(Rk\right)\right]\\ =\frac{\hslash c\pi }{{\left(2\pi \right)}^2}\frac{1}{R^4}\stackrel{\infty }{\underset{0}{{\displaystyle \int }}}\text{d}w{w}^3\text{cos}\left(w\frac{ct}{R}\right)\times \left[\frac{4}{w}\text{sin}\left(w\right)+\frac{4}{w^2}\text{cos}\left(w\right)-\frac{4}{w^3}\text{sin}\left(w\right)\right].\end{array}(76)$

Substituting $\frac{ct}{R}=b$, and using an integral table [36] (p. 504, No. 8), also discussed in the limiting sense in Ref. [7], Appendix C:

$\underset{0}{\overset{\infty }{{\displaystyle \int }}}\text{d}w\text{cos}\left(wb\right)\text{sin}\left(w\right)=\frac{1}{\left(1-b^2\right)},(77)$

$\underset{0}{\overset{\infty }{{\displaystyle \int }}}\text{d}ww\text{cos}\left(wb\right)\text{cos}\left(w\right)=-\left[\frac{1}{\left(1-b^2\right)}+\frac{2b^2}{{\left(1-b^2\right)}^2}\right],(78)$

$\underset{0}{\overset{\infty }{{\displaystyle \int }}}\text{d}w{w}^2\text{cos}\left(wb\right)\text{sin}\left(w\right)=-2\left[\frac{1}{\left(1-b^2\right)}+\frac{5b^2}{{\left(1-b^2\right)}^2}+\frac{4b^4}{{\left(1-b^2\right)}^3}\right],(79)$

enables Eq. 76 to be evaluated:

$\begin{array}{l}\langle E_{\text{ZP},x}\left(\mathrm{0,0}\right)E_{\text{ZP},x}\left(\widehat{\mathbf{y}}R,t\right)\rangle \\ =\frac{\hslash c\pi }{{\left(2\pi \right)}^2}\frac{1}{R^4}4\left\{-2\left[\frac{1}{\left(1-b^2\right)}+\frac{5b^2}{{\left(1-b^2\right)}^2}+\frac{4b^4}{{\left(1-b^2\right)}^3}]-[\frac{1}{\left(1-b^2\right)}+\frac{2b^2}{{\left(1-b^2\right)}^2}\right]-\frac{1}{\left(1-b^2\right)}\right\}\end{array}(80)$

$\begin{array}{l}=\frac{\hslash c\pi }{{\left(2\pi \right)}^2}\frac{1}{R^4}4\frac{1}{{\left(1-b^2\right)}^3}\left\{-2{\left(1-b^2\right)}^2-10b^2\left(1-b^2\right)-8b^4-{\left(1-b^2\right)}^2-2b^2\left(1-b^2\right)-{\left(1-b^2\right)}^2\right\}\\ =\frac{\hslash c\pi }{{\left(2\pi \right)}^2}\frac{1}{R^4}4\frac{1}{{\left(1-b^2\right)}^3}\left(-4b^2-4\right)=-4\frac{\hslash c\pi }{{\pi }^2}\frac{1}{R^4}\frac{\left[1+{\left(\frac{ct}{R}\right)}^2\right]}{{\left(1-{\left(\frac{ct}{R}\right)}^2\right)}^3}\\ =-\frac{4\hslash c}{\pi }\frac{\left[R^2+{\left(ct\right)}^2\right]}{{[R^2-{\left(ct\right)}^2]}^3}\end{array}$

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 “n-point 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.

$\begin{array}{c}\langle E_{\text{ZP},x}\left(\mathrm{0,0}\right)E_{\text{ZP},x}\left(\widehat{\mathbf{y}}R,t\right)\rangle \\ =\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}{d{E}_{2x}{E}_{1x}{E}_{2x}{I}_{1x2x}|}_{\text{at }{\mathbf{x}}_1,t_1\& {\mathbf{x}}_2,t_2}\\ =\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}{E}_{1x}{E}_{2x}{\left\{\frac{1}{2\pi A_x{\left(1-\frac{C_x^2}{A_x^2}\right)}^{1/2}}\left[\frac{-E_{1x}^2-E_{2x}^2+2E_{2x}{E}_{1x}\left(\frac{C_x}{A_x}\right)}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\right\}}_{{\mathbf{x}}_1,t_1;{\mathbf{x}}_2,t_2} .\\ \end{array}$

Here, the meaning of the subscript at the end of ${\mathbf{x}}_1,t_1$; ${\mathbf{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, ${\mathbf{x}}_1,t_1$ and ${\mathbf{x}}_2,t_2$, respectively, in the function $C_{12,\mathbf{n}}$, Eq. 47, contained within $C_x$, in Eq. 60.

Thus,

$\begin{array}{c}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}\left(E_{1x}{E}_{2x}\right)I_{1x2x}\\ ={\left\{\begin{array}{c}\frac{1}{2\pi A_x{\left(1-\frac{C_x^2}{A_x^2}\right)}^{1/2}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d}{E}_{2x}{E}_{2x}\exp \left[-\frac{E_{2x}^2}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\times \\ \times {\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d}{E}_{1x}{E}_{1x}\exp \left[\frac{-E_{1x}^2+2E_{2x}{E}_{1x}\frac{C_x}{A_x}}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\end{array}\right\}}_{{\mathbf{x}}_1,t_1\text{;}{\mathbf{x}}_2,t_2}\\ ={\left\{\begin{array}{c}\frac{1}{2\pi A_x{\left(1-\frac{C_x^2}{A_x^2}\right)}^{1/2}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d}{E}_{2x}{E}_{2x}\exp \left[\frac{-E_{2x}^2+E_{2x}^2\frac{C_x^2}{A_x^2}}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\times \\ \times {\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d}{E}_{1x}{E}_{1x}\exp \left[\frac{-{\left(E_{1x}-E_{2x}\frac{C_x}{A_x}\right)}^2}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\end{array}\right\}}_{{\mathbf{x}}_1,t_1\text{;}{\mathbf{x}}_2,t_2}\\ ={\left\{\frac{1}{2\pi A_x{\left(1-\frac{C_x^2}{A_x^2}\right)}^{1/2}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}{E}_{2x}\text{exp}\left(\frac{-E_{2x}^2}{2A_x}\right)\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}{E}_{1x}\text{exp}\left[\frac{-{\left(E_{1x}-E_{2x}\frac{C_x}{A_x}\right)}^2}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\right\}}_{{\mathbf{x}}_1,t_1;{\mathbf{x}}_2,t_2}\end{array}$

The integral over $E_{1x}$ on the right can be broken up as:

$\begin{array}{c}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}{E}_{1x}\text{exp}{\left[\frac{-{\left(E_{1x}-E_{2x}\frac{C_x}{A_x}\right)}^2}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]}_{{\mathbf{x}}_1,t_1;{\mathbf{x}}_2,t_2}\\ =\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}{\left\{\left(E_{1x}-E_{2x}\frac{C_x}{A_x}\right)\text{exp}\left[\frac{-{\left(E_{1x}-E_{2x}\frac{C_x}{A_x}\right)}^2}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\right\}}_{{\mathbf{x}}_1,t_1;{\mathbf{x}}_2,t_2}\\ +\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}{E}_{2x}{\left\{\frac{C_x}{A_x}\text{exp}\left[\frac{-{\left(E_{1x}-E_{2x}\frac{C_x}{A_x}\right)}^2}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\right\}}_{{\mathbf{x}}_1,t_1;{\mathbf{x}}_2,t_2}\\ =\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}dx\cdot x\text{exp}{\left[\frac{-x^2}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]}_{{\mathbf{x}}_1,t_1;{\mathbf{x}}_2,t_2}+{\left\{E_{2x}\frac{C_x}{A_x}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}\text{exp}\left[\frac{-{\left(E_{1x}-E_{2x}\frac{C_x}{A_x}\right)}^2}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]\right\}}_{{\mathbf{x}}_1,t_1;{\mathbf{x}}_2,t_2}\end{array}$

The first integral equals zero, as it is odd in x. Hence:

$\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}{E}_{1x}{E}_{1x}\text{exp}{\left[\frac{-{\left(E_{1x}-E_{2x}\frac{C_x}{A_x}\right)}^2}{2A_x\left(1-\frac{C_x^2}{A_x^2}\right)}\right]}_{{\mathbf{x}}_1,t_1;{\mathbf{x}}_2,t_2}=E_{2x}{\left\{\frac{C_x}{A_x}{\left[\pi 2A_x\left(1-\frac{C_x^2}{A_x^2}\right)\right]}^{1/2}\right\}}_{{\mathbf{x}}_1,t_1;{\mathbf{x}}_2,t_2}$

and

$\begin{array}{c}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}\left(E_{1x}{E}_{2x}\right)I_{1x2x}\\ ={\left[\frac{1}{2\pi A_x{\left(1-\frac{C_x^2}{A_x^2}\right)}^{1/2}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}{E}_{2x}\text{exp}\left(-\frac{E_{2x}^2}{2A_x}\right)\times \left\{E_{2x}\frac{C_x}{A_x}{\left[\pi 2A_x\left(1-\frac{C_x^2}{A_x^2}\right)\right]}^{1/2}\right\}\right]}_{{\mathbf{x}}_1,t_1;{\mathbf{x}}_2,t_2}\\ ={\left\{\frac{1}{2\pi A_x{\left(1-\frac{C_x^2}{A_x^2}\right)}^{1/2}}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}{E}_{2x}^2\text{exp}\left(-\frac{E_{2x}^2}{2A_x}\right)\frac{C_x}{A_x}{\left[\pi 2A_x\left(1-\frac{C_x^2}{A_x^2}\right)\right]}^{1/2}\right\}}_{{\mathbf{x}}_1,t_1;{\mathbf{x}}_2,t_2}\\ =\frac{1}{{\left(2\pi A_x\right)}^{1/2}}\frac{{\left(C_x\right)}_{{\mathbf{x}}_1,t_1;{\mathbf{x}}_2,t_2}}{A_x}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}{E}_{2x}^2\text{exp}\left(-\frac{E_{2x}^2}{2A_x}\right)\\ =\frac{1}{{\left(2\pi A_x\right)}^{1/2}}\frac{{\left(C_x\right)}_{{\mathbf{x}}_1,t_1;{\mathbf{x}}_2,t_2}}{A_x}\frac{{\pi }^{1/2}}{2}{\left(2A_x\right)}^{3/2}={\left(C_x\right)}_{{\mathbf{x}}_1,t_1;{\mathbf{x}}_2,t_2}\\ =\frac{1}{\left(L_x{L}_y{L}_z\right)}\stackrel{\infty }{{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }}{\left(C_{12,\mathbf{n}}\right)}_{{\mathbf{x}}_1,t_1;{\mathbf{x}}_2,t_2}{K}_{\mathbf{n},x}^2{\sigma }_{\mathbf{n}}^2.\end{array}$

The last expression for $C_x$ came from Eq. 60.

Hence, using Eqs. 47 and 54 replacing the coordinates ${\mathbf{x}}_1$ and $t_1$ with 0, 0, and ${\mathbf{x}}_2$ and $t_2$ with $\widehat{\mathbf{y}}R$ and t, respectively, and ${\sigma }_{\mathbf{n}}^2$ with ZP of $2\pi \hslash {\omega }_{\mathbf{n},\lambda }$:

$\begin{array}{l}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{E}_{2x}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}\text{d}{E}_{1x}\left(E_{1x}{E}_{2x}\right)I_{1x2x}\\ =\frac{1}{\left(L_x{L}_y{L}_z\right)}\stackrel{\infty }{{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }}\text{cos}\left[{\mathbf{k}}_{\mathbf{n}}\cdot \left({\mathbf{x}}_1-{\mathbf{x}}_2\right)-{\omega }_{\mathbf{n}}\left(t_1-t_2\right)\right]\left[1-{\left(\frac{{\mathbf{k}}_{\mathbf{n},x}}{k_{\mathbf{n}}}\right)}^2\right]{\sigma }_{\mathbf{n}}^2\\ =\frac{1}{\left(L_x{L}_y{L}_z\right)}\stackrel{\infty }{{\displaystyle \sum }_{n_x,n_y,n_z=-\infty }}\text{cos}\left[k_{y}R-ckt\right]\left[1-{\left(\frac{{\mathbf{k}}_{\mathbf{n},x}}{k_{\mathbf{n}}}\right)}^2\right]2\pi \hslash {\omega }_{\mathbf{n},\lambda }.\end{array}$

Implementing the same continuum approximation as with the other method, Eq. 70, then results in:

$\begin{array}{c}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}\left(E_{1x}{E}_{2x}\right)I_{1x2x}\\ \approx \frac{1}{\left(2\pi \right)}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{k}_x\frac{1}{\left(2\pi \right)}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{k}_y\frac{1}{\left(2\pi \right)}\times \stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{k}_z\text{cos}\left(k_{y}R-ckt\right)\left[1-{\left(\frac{k_x}{k}\right)}^2\right]2\pi \hslash \omega \\ =\frac{1}{{\left(2\pi \right)}^3}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{k}_x\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{k}_y\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{k}_z\left[\text{cos}\left(k_{y}R\right)\text{cos}\left(ckt\right)+\text{sin}\left(k_{y}R\right)\text{sin}\left(ckt\right)\right]\left[1-{\left(\frac{k_x}{k}\right)}^2\right]2\pi \hslash \omega \end{array}$

The second integral with $\text{sin}\left(k_{y}R\right)$ is odd in $k_y$ and equals zero by symmetry. Hence:

$\begin{array}{c}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{2x}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{E}_{1x}\left(E_{1x}{E}_{2x}\right)I_{1x2x}\\ =\frac{1}{{\left(2\pi \right)}^3}\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{k}_x\stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{k}_y\times \stackrel{\infty }{\underset{-\infty }{{\displaystyle \int }}}d{k}_z\text{cos}\left(k_{y}R\right)\text{cos}\left(\omega t\right)\left[1-{\left(\frac{k_x}{k}\right)}^2\right]2\pi \hslash \omega \end{array}(82)$

$\begin{array}{c}=\frac{1}{{\left(2\pi \right)}^3}\stackrel{\infty }{\underset{0}{{\displaystyle \int }}}{k}^{2}dk\stackrel{\pi }{\underset{0}{{\displaystyle \int }}}d\theta \text{sin}\theta \stackrel{2\pi }{\underset{0}{{\displaystyle \int }}}d\varphi \text{cos}\left(kR\text{cos}\theta \right)\text{cos}\left(\omega t\right)\times \left(1-{\text{sin}}^2\theta {\text{sin}}^2\varphi \right)\left(2\pi \hslash kc\right)\\ =\frac{\hslash c}{{\left(2\pi \right)}^2}\stackrel{\infty }{\underset{0}{{\displaystyle \int }}}{k}^{3}dk\stackrel{\pi }{\underset{0}{{\displaystyle \int }}}d\theta \text{sin}\theta \times \stackrel{2\pi }{\underset{0}{{\displaystyle \int }}}d\varphi \text{cos}\left(kR\text{cos}\theta \right)\text{cos}\left(\omega t\right)\left(1-{\text{sin}}^2\theta {\text{sin}}^2\varphi \right).\end{array}(83)$

However, this result agrees exactly with the earlier result, Eq. 71, obtained partway through the ensemble derivation of $\langle E_{\text{ZP},x}\left(\mathrm{0,0}\right)E_{\text{ZP},x}\left(\widehat{\mathbf{y}}R,t\right)\rangle$. Thus, continuing with further steps in evaluating Eq. 83 will provide a final result, using the joint probability density approach of $\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}{E}_{2x}\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\text{d}{E}_{1x}\left(E_{1x}{E}_{2x}\right)I_{1x2x}$, that exactly agrees with the ensemble derivation $\langle E_{\text{ZP},x}\left(\mathrm{0,0}\right)E_{\text{ZP},x}\left(\widehat{\mathbf{y}}R,t\right)\rangle$ of $-\frac{4\hslash c}{\pi }\frac{\left[R^2+{\left(ct\right)}^2\right]}{{[R^2-{\left(ct\right)}^2]}^3}$ in Eq. 81.

2.4 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_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ and $B_{{\mathbf{k}}_{\mathbf{n}},\lambda }$. To put this in better perspective, if we imagine an ensemble of boxes $L_x\times L_y\times L_z$, each the same size, but existing at different points in space and/or in time, then electromagnetic field fluctuations of ${\mathbf{E}}_{\text{ZPP}}$ and ${\mathbf{B}}_{\text{ZPP}}$ will occur at each point within each box. However, for each box, there is only one set of coefficients $A_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ and $B_{{\mathbf{k}}_{\mathbf{n}},\lambda }$, 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_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ or $B_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ is given by Eq. 8, while the variance of each is given by ${\left[\sigma \left({\omega }_{\mathbf{n}},T\right)\right]}^2$, as in Eq. 8.

Thus, ${\mathbf{E}}_{\text{ZPP}}$ and ${\mathbf{B}}_{\text{ZPP}}$ can be viewed as the linear transformation of the random variables $A_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ and $B_{{\mathbf{k}}_{\mathbf{n}},\lambda }$. However, the coefficients multiplying $A_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ and $B_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ in these linear sums in Eqs. 1 and 2 are not constants, as they depend on time and space. In particular, $\frac{1}{{\left(L_x{L}_y{L}_z\right)}^{1/2}}{\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\text{cos}\left({\mathbf{k}}_{\mathbf{n}}\cdot \mathbf{x}-{\omega }_{\mathbf{n}}t\right)$ and $\frac{1}{{\left(L_x{L}_y{L}_z\right)}^{1/2}}{\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }\text{sin}\left({\mathbf{k}}_{\mathbf{n}}\cdot \mathbf{x}-{\omega }_{\mathbf{n}}t\right)$ are the coefficients multiplying the random variables of $A_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ and $B_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ for $\mathbf{E}\left(\mathbf{x},t\right)$ in Eq. 1. An exactly similar situation occurs for $A_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ and $B_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ regarding $\mathbf{B}\left(\mathbf{x},t\right)$ in Eq. 2, except that ${\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ is replaced by ${\widehat{\mathbf{k}}}_{\mathbf{n}}\times {\widehat{\epsilon }}_{{\mathbf{k}}_{\mathbf{n}},\lambda }$. Thus, although the probabilistic properties for the random variables $A_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ and $B_{{\mathbf{k}}_{\mathbf{n}},\lambda }$ are independent of time and space, as expressed by Eq. 8, the same is not true for the probabilistic/stochastic properties of ${\mathbf{E}}_{\text{ZPP}}$ and ${\mathbf{B}}_{\text{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:

$P\left(\begin{array}{c}{\mathbf{E}}_1\text{ at }{\mathbf{x}}_1,t_1;{\mathbf{E}}_2\text{ at }{\mathbf{x}}_2,t_2;\dots ;{\mathbf{E}}_n\text{ at }{\mathbf{x}}_n,t_n;{\mathbf{B}}_{n+1}\text{ at }{\mathbf{x}}_{n+1},t_{n+1};\\ {\mathbf{B}}_{n+2}\text{ at }{\mathbf{x}}_{n+2},t_{n+2};\dots ;{\mathbf{B}}_{n+m}\text{ at }{\mathbf{x}}_{n+m},t_{n+m}\end{array}\right),(84)$

where there are n electric field vector values (i.e., $3\times n$ component values, as indicated below) and m magnetic field vector values ($3\times m$ component values), at respective positions in space and time, as indicated. However, since this is a multivariate normal distribution, where all $\langle {\mathbf{E}}_i\rangle$ and all $\langle {\mathbf{B}}_i\rangle$ ensemble averages equal zero, then by probability theory, the above would be represented by [37]:

$P\left(\mathbf{X}\right)=\frac{\text{exp}\left[-\frac{1}{2}{\mathbf{X}}^T{\text{Σ}}^{-1}\mathbf{X}\right]}{\sqrt{{\left(2\pi \right)}^{n+m}\left|\text{Σ}\right|}},(85)$

where $\mathbf{X}$ is the vector of

$\left(E_{1x},E_{1y},E_{1z},\dots ,E_{nx},E_{ny},E_{nz},B_{\left(n+1\right)x},B_{\left(n+1\right)y},B_{\left(n+1\right)z},\dots ,B_{\left(n+m\right)x},B_{\left(n+m\right)y},B_{\left(n+m\right)z}\right)$

values. Moreover, $\text{Σ}$ is the covariant matrix as expressed by

${\text{Σ}}_{ij}=\langle X_i{X}_j\rangle ,(86)$

since $\langle X_i\rangle =0$ for all $E_i$ and $B_i$ values, due to Eqs. 4, 1 and 2. Also, $\left|\text{Σ}\right|$ and ${\text{Σ}}^{-1}$ represent the determinant and inverse matrix of the covariant matrix, $\text{Σ}$, respectively.

If we calculated all the components of the covariant matrix, meaning all combinations ${\text{Σ}}_{ij}=\langle X_i{X}_j\rangle$ of pairs of electric and magnetic field expectation values, then the probability density, Eq. 85, could be evaluated for any vector $\mathbf{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\left(\mathbf{E}\text{ at }\mathbf{x},t\right)$ in Eq. 34, but this also follows from Eq. 85 with $\mathbf{X}$ being the vector $\left(E_{1x},E_{1y},E_{1z}\right)$, then using Eq. 86 plus earlier relations in Section 2.1, and

${\text{Σ}}_{ij}=\left[\begin{array}{ccc}\langle E_{1x}^2\rangle & 0& 0\\ 0& \langle E_{1y}^2\rangle & 0\\ 0& 0& \langle E_{1z}^2\rangle \end{array}\right],$

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]:

$P\left(E_{1x},E_{2x}\right)=\frac{\text{exp}\left\{-\frac{1}{2\left(1-{\rho }_{1x2x}^2\right)}\left[\frac{{\left(E_{1x}-{\mu }_{1x}\right)}^2}{{\sigma }_{1x}^2}+\frac{{\left(E_{2x}-{\mu }_{2x}\right)}^2}{{\sigma }_{2x}^2}-\frac{2{\rho }_{1x2x}\left(E_{1x}-{\mu }_{1x}\right)\left(E_{2x}-{\mu }_{2x}\right)}{{\sigma }_{1x}{\sigma }_{2x}}\right]\right\}}{2\pi {\sigma }_{1x}{\sigma }_{2x}\sqrt{1-{\rho }_{1x2x}^2}}.(87)$

Here, ${\mu }_{1x}\equiv \langle E_{1x}\rangle$ and ${\mu }_{2x}\equiv \langle E_{2x}\rangle$ are both zero for our ZP and ZPP cases. Also, ${\sigma }_{1x}^2\equiv \langle {\left(E_{1x}-{\mu }_{1x}\right)}^2\rangle =\langle E_{1x}^2\rangle$, for our situation, and similarly for ${\left({\sigma }_{2x}\right)}^2=\langle {\left(E_{2x}\right)}^2\rangle ,$ while

${\rho }_{1x2x}\equiv \frac{\langle \left(E_{1x}-{\mu }_{1x}\right)\left(E_{2x}-{\mu }_{2x}\right)\rangle }{{\sigma }_{1x}{\sigma }_{2x}}(88)$

is called the Pearson’s correlation coefficient of $E_{1x}$ and $E_{2x}$. With ${\mu }_{1x}={\mu }_{2x}=0$ for our case, then

${\rho }_{1x2x}\to \frac{\langle E_{1x}{E}_{2x}\rangle }{\sqrt{\langle E_{1x}^2\rangle \langle E_{2x}^2\rangle }}.(89)$

Putting these expressions into Eq. 87 results in:

$P\left(E_{1x},E_{2x}\right)=\frac{\exp \left\{-\frac{1}{2\left(1-\frac{{\langle E_{1x}{E}_{2x}\rangle }^2}{\langle E_{1x}^2\rangle \langle E_{2x}^2\rangle }\right)}\left[\frac{E_{1x}^2}{\langle E_{1x}^2\rangle }+\frac{E_{2x}^2}{\langle E_{2x}^2\rangle }-2E_{1x}{E}_{2x}\frac{\langle E_{1x}{E}_{2x}\rangle }{\langle E_{1x}^2\rangle \langle E_{2x}^2\rangle }\right]\right\}}{2\pi \langle E_{1x}^2\rangle \langle E_{2x}^2\rangle \sqrt{1-\frac{{\langle E_{1x}{E}_{2x}\rangle }^2}{\langle E_{1x}^2\rangle \langle E_{2x}^2\rangle }}}(90)$

This is then readily related to our result of $I_{1x2x}$ in Section 2.2, Eq. 62, since from Eq. 66, $A_i=\langle E_{\text{ZPP},i}^2\rangle$ for ($i\to x,y,z$), and is independent of space and time, so $\langle E_{1x}^2\rangle =\langle E_{2x}^2\rangle$, 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, $\langle E_{1x}{E}_{2x}\rangle$. Thus,

$\frac{C_x}{A_x}=\frac{\langle E_{1x}{E}_{2x}\rangle }{\sqrt{\langle E_{1x}^2\rangle \langle E_{2x}^2\rangle }},$

which is just ${\rho }_{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:

$\begin{array}{l}P\left({\mathbf{E}}_1\text{ at }{\mathbf{x}}_1,t_1;{\mathbf{E}}_2\text{ at }{\mathbf{x}}_2,t_2\right)\\ =I_{1x2x}\cdot I_{1y2y}\cdot I_{1z2z}\\ =\frac{\exp \left[\frac{\mathrm{-}{E}_{1x}^2\mathrm{-}{E}_{2x}^2+2E_{1x}{E}_{2x}{\rho }_{1x,2x}}{2\left(1\mathrm{-}{\rho }_{1x,2x}^2\right)\langle E_{1x}^2\rangle }\right]}{2\pi \langle E_{1x}^2\rangle {\left(1\mathrm{-}{\rho }_{1x,2x}^2\right)}^{1/2}}\times \frac{\exp \left[\frac{\mathrm{-}{E}_{1y}^2\mathrm{-}{E}_{2y}^2+2E_{1y}{E}_{2y}{\rho }_{1y,2y}}{2\left(1\mathrm{-}{\rho }_{1y,2y}^2\right)\langle E_{1y}^2\rangle }\right]}{2\pi \langle E_{1y}^2\rangle {\left(1\mathrm{-}{\rho }_{1y,2y}^2\right)}^{1/2}}\times \frac{\exp \left[\frac{\mathrm{-}{E}_{1z}^2\mathrm{-}{E}_{2z}^2+2E_{1z}{E}_{2z}{\rho }_{1z,2z}}{2\left(1\mathrm{-}{\rho }_{1z,2z}^2\right)\langle E_{1z}^2\rangle }\right]}{2\pi \langle E_{1z}^2\rangle {\left(1\mathrm{-}{\rho }_{1z,2z}^2\right)}^{1/2}}\end{array}(91)$

3 Concluding Remarks

The technique used here in Eq. 11 for one electric field value $\mathbf{E}$, or three component values $E_x$, $E_y$, $E_z$, or with Eq. 40 for $P\left({\mathbf{E}}_1\text{ at }{\mathbf{x}}_1,t_1;{\mathbf{E}}_2\text{ at }{\mathbf{x}}_2,t_2\right)$, 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,

$\begin{array}{c}P\left(x_1\text{ at }{t}_1;x_2\text{ at }{t}_2\right)=\\ {{\displaystyle \int }}^{\text{}}d{A}_1\cdots {{\displaystyle \int }}^{\text{}}d{A}_N\dots {{\displaystyle \int }}^{\text{}}d{B}_1\cdots {{\displaystyle \int }}^{\text{}}d{B}_N\dots P\left(A_{1,}\cdots ,A_N,\cdots ,B_N,\dots \right)\delta \left[x_1-x\left(t_1\right)\right]\delta \left[x_2-x\left(t_2\right)\right] ,\end{array}(92)$

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\left(t_1\right)$ and $x\left(t_2\right)$ 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.

Data Availability Statement

All datasets presented in this study are included in the article/ Supplementary Material.

Author Contributions

The author confirms being the sole contributor of this work and has approved it for publication.

Conflict of Interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

1. Cole DC. Reinvestigation of the thermodynamics of blackbody radiation via classical physics. Phys Rev A (1992) 45:8471–89. doi:10.1103/physreva.45.8471

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Cole DC. Thermodynamics of blackbody radiation via classical physics for arbitrarily shaped cavities with perfectly conducting walls. Found Phys (2000) 30(11):1849–67. doi:10.1023/a:1003706320972

CrossRef Full Text | Google Scholar

3. Boyer TH. Thermal effects of acceleration for a classical dipole oscillator in classical electromagnetic zero-point radiation. Phys Rev D (1984) 29(6):1089–95. doi:10.1103/physrevd.29.1089

CrossRef Full Text | Google Scholar

4. Boyer TH. Thermal effects of acceleration through random classical radiation. Phys Rev D (1980) 21(8):2137–48. doi:10.1103/physrevd.21.2137

CrossRef Full Text | Google Scholar

5. Boyer TH. Thermal effects of acceleration for a classical spinning magnetic dipole in classical electromagnetic zero-point radiation. Phys Rev D (1984) 30(6):1228–32. doi:10.1103/physrevd.30.1228

CrossRef Full Text | Google Scholar

6. Cole DC. Properties of a classical charged harmonic oscillator accelerated through classical electromagnetic zero-point radiation. Phys Rev D (1985) 31(8):1972–81. doi:10.1103/physrevd.31.1972

CrossRef Full Text | Google Scholar

7. Cole DC. Thermal effects of acceleration for a spatially extended electromagnetic system in classical electromagnetic zero-point radiation: transversely positioned classical oscillators. Phys Rev D (1987) 35:562–83. doi:10.1103/physrevd.35.562

CrossRef Full Text | Google Scholar

8. Boyer TH. Temperature dependence of Van der Waals forces in classical electrodynamics with classical electromagnetic zero-point radiation. Phys Rev A (1975) 11:1650–63. doi:10.1103/physreva.11.1650

CrossRef Full Text | Google Scholar

9. Boyer T. Stochastic electrodynamics: the closest classical approximation to quantum theory. Atoms (2019) 7(1):29–39. doi:10.3390/atoms7010029

CrossRef Full Text | Google Scholar

10. Marshall TW. Random electrodynamics. Proc R Soc London, Ser A (1963) 276:475–91.

Google Scholar

11. Marshall TW. A classical treatment of blackbody radiation. Nuovo Cim (1965) 38:206–15. doi:10.1007/bf02750449

CrossRef Full Text | Google Scholar

12. Boyer TH. Derivation of the blackbody radiation spectrum without quantum assumptions. Phys Rev (1821) 1374–83. doi:10.1103/phyRev.182.1374

Google Scholar

13. Boyer TH. Classical statistical thermodynamics and electromagnetic zero-point radiation. Phys Rev (1861) 1304–18. doi:10.1103/phyRev.186.1304

Google Scholar

14. Cole DC, Zou Y. Simulation study of aspects of the classical hydrogen atom interacting with electromagnetic radiation: circular orbits. J Scientific Comput (2004) 20(1):43–68. doi:10.1023/a:1025846412872

CrossRef Full Text | Google Scholar

15. Jackson JD. Classical electrodynamics. 2nd ed.. New York , NY John Wiley & Sons (1975).

16. de la Peña L, Cetto AM. The quantum dice - an introduction to stochastic electrodynamics. Kluwer Dordrecht: Kluwer Acad. Publishers (1996).

17. Boyer TH. Random electrodynamics: the theory of classical electrodynamics with classical electromagnetic zero-point radiation. Phys Rev D (1975) 11(4):790–808. doi:10.1103/physrevd.11.790

CrossRef Full Text | Google Scholar

18. Cole DC, Lakhtakia A. Compendium book, “essays on formal aspects of electromagnetic theory. edited. Singapore: World Scientific (1993). p. 501–32.

CrossRef Full Text

19. Boyer TH. The classical vacuum. Sci Am (1985) 253:70–8. doi:10.1038/scientificamerican0885-70

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Marshall TW. Statistical electrodynamics. Math Proc Camb Phil Soc (1965) 61:537–46. doi:10.1017/s0305004100004114

CrossRef Full Text | Google Scholar

21. Cole DC. Derivation of the classical electromagnetic zero-point radiation spectrum via a classical thermodynamic operation involving van der Waals forces. Phys Rev A (1990) 42:1847–62. doi:10.1103/physreva.42.1847

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Cole DC. Entropy and other thermodynamic properties of classical electromagnetic thermal radiation. Phys Rev A (1990) 42:7006–24. doi:10.1103/physreva.42.7006

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Cole DC. Energy considerations of classical electromagnetic zero-point radiation and a specific probability calculation in stochastic electrodynamics. Atoms (2019) 7(2):50. doi:10.3390/atoms7020050

CrossRef Full Text | Google Scholar

24. Boyer TH. `General connection between random electrodynamics and quantum electrodynamics for free electromagnetic fields and for dipole oscillator systems. Phys Rev D (1975) 11(4):809–30. doi:10.1103/physrevd.11.809

CrossRef Full Text | Google Scholar

25. Planck M. The theory of Heat radiation. New York: Dover (1959). This publication is an English translation of the second edition of Planck’s work entitled Waermestrahlung, published in 1913. A more recent republication of this work is Vol. 11 of the series The History of Modern Physics 1800–1950 (AIP, New York, 1988).

26. Einstein A, Hopf L. Über einen Satz der Wahrscheinlichkeitsrechnung und seine Anwendung in der Strahlungstheorie. Ann Phys (1910) 338:1096–104.doi:10.1002/andp.19103381603

CrossRef Full Text | Google Scholar

27. Einstein A, Hopf L. Statistische Untersuchung der Bewegung eines Resonators in einem Strahlungsfeld. Ann Phys (1910) 338:1105–15. doi:10.1002/andp.19103381604

CrossRef Full Text | Google Scholar

28. Cole DC, Zou Y. Quantum mechanical ground state of hydrogen obtained from classical electrodynamics. Phys Lett A (2003) 317(1–2):14–20. doi:10.1016/j.physleta.2003.08.022

CrossRef Full Text | Google Scholar

29. Cole DC. Simulation results related to stochastic electrodynamics. In AIP Conference Proceedings, USA. 810 99 – 113(.) 2006.

Google Scholar

30. Nieuwenhuizen TM, Liska MTP. Simulation of the hydrogen ground state in stochastic electrodynamics. Phys Scr (2015) T165:014006. doi:10.1088/0031-8949/2015/t165/014006

CrossRef Full Text | Google Scholar

31. Nieuwenhuizen TM, Liska MTP. Simulation of the hydrogen ground state in stochastic electrodynamics-2: inclusion of relativistic corrections. Found Phys (2015) 45(10):1190–202. doi:10.1007/s10701-015-9919-0

CrossRef Full Text | Google Scholar

32. Boyer TH. Conformal symmetry of classical electromagnetic zero-point radiation. Found Phys (1989) 19(4):349–65. doi:10.1007/bf00731830

CrossRef Full Text | Google Scholar

33. Haroche S, Kleppner D. Cavity quantum electrodynamics. Phys Today (1989) 42(1):24–30. doi:10.1063/1.881201

CrossRef Full Text | Google Scholar

34. Haroche S, Raimond J-M. Cavity quantum electrodynamics. Sci Am (1993) 268(4):26–33. doi:10.1038/scientificamerican0493-54

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Cole DC. Correlation functions for homogeneous, isotropic random classical electromagnetic radiation and the electromagnetic fields of a fluctuating classical electric dipole. Phys Rev D (1986) 33:2903–15. doi:10.1103/physrevd.33.2903

CrossRef Full Text | Google Scholar

36. Gradshteyn IS, Ryzhik IM. Tables of integrals, series, and products. New York: Academic (1980).

37. Taboga M. Lectures on probability theory and mathematical statistics. 3rd ed.. Scotts Valley, California, US: CreateSpace Independent Publishing Platform (2017).

38. Shirazi AN, Fleck I. Bivariate normal distribution for finding inliers in hough space for a time projection chamber, 150. EPJ Web Conf (2017). 00010.

39. Feynman RP, Hibbs AR. Quantum mechanics and path integrals. New York: McGraw-Hill (1965).