Space and Time Averaged Quantum Stress Tensor Fluctuations

We extend previous work on the numerical diagonalization of quantum stress tensor operators in the Minkowski vacuum state, which considered operators averaged in a finite time interval, to operators averaged in a finite spacetime region. Since real experiments occur over finite volumes and durations, physically meaningful fluctuations may be obtained from stress tensor operators averaged by compactly supported sampling functions in space and time. The direct diagonalization, via a Bogoliubov transformation, gives the eigenvalues and the probabilities of measuring those eigenvalues in the vacuum state, from which the underlying probability distribution can be constructed. For the normal-ordered square of the time derivative of a massless scalar field in a spherical cavity with finite degrees of freedom, analysis of the tails of these distributions confirms previous results based on the analytical treatment of the high moments. We find that the probability of large vacuum fluctuations is reduced when spatial averaging is included, but the tail still decreases more slowly than exponentially as the magnitude of the measured eigenvalues increases, suggesting vacuum fluctuations may not always be subdominant to thermal fluctuations and opening up the possibility of experimental observation under the right conditions.

averaged by compactly supported sampling functions in space and time. The direct diagonalization, via a Bogoliubov transformation, gives the eigenvalues and the probabilities of measuring those eigenvalues in the vacuum state, from which the underlying probability distribution can be constructed. For the normal-ordered square of the time derivative of a massless scalar field in a spherical cavity with finite degrees of freedom, analysis of the tails of these distributions confirms previous results based on the analytical treatment of the high moments. We find that the probability of large vacuum fluctuations is reduced when spatial averaging is included, but the tail still decreases more slowly than exponentially as the magnitude of the measured eigenvalues increases, suggesting vacuum fluctuations may not always be subdominant to thermal fluctuations and opening up the possibility of experimental observation under the right conditions.

I. INTRODUCTION
The semiclassical theory of gravity, where matter fields are treated as quantum fields, whereas the gravitational field is treated as a classical field, deals with the expectation value of the energy-momentum tensor operator of the matter fields as an approximation to a full theory of quantum gravity [1]. This semiclassical theory offers a plausible description of the backreaction of the Hawking radiation on the gravitational background of a black hole [2] and opens the possibility of quantum particle creation [3] through higher order derivative terms of the metric [4].
Nevertheless, the semiclassical theory does not provide a description of the expected quantum fluctuations of the stress tensor around its expectation value. Such fluctuations may have observable physical effects and, in recent years, has captured a great amount of attention from the physics community [5][6][7][8][9][10][11][12]. Generally speaking, physical effects of the quantum fluctuations of a stress tensor operator may be addressed via the calculation of the probability distribution of the time-or spacetime-averaged operator. Since these averages require normal ordering, the probability distribution in the vacuum state has zero mean and thus a nonzero probability of measuring negative components of the stress-energy tensor, such as the energy density. In some sense, the averaging process may be viewed as a consequence of a physical measurement probing the outcomes of the operator.
The exact probability distribution associated with measurements of the stress-energy tensor in a two-dimensional conformal field theory in the vacuum is calculated in Ref. [13].
Using a Gaussian temporal sampling function, the resulting distribution is a shifted Gamma distribution, with the shift given by the optimal quantum inequality bound [14]. In four dimensions, the situation is more involved. Qualitatively, the probability distributions may be inferred from the moments of the averaged operator, as done for several normal-ordered quadratic operators in the vacuum state using Lorentzian time averaging in Ref. [15] and later generalized in Ref. [16] for compactly supported functions, which are functions that are exactly zero outside a defined domain. The main prediction in both references is the asymptotic form of the distribution, which represents the probability of large fluctuations, falling more slowly than exponentially as where x refers to a dimensionless measurement of the stress tensor fluctuations and a, b, c, c 0 are constants that vary according to the smearing function and the specific operator. Take u to be the Lorentzian time average of the electromagnetic energy density in a timescale τ , so that a dimensionless measurement of the averaged energy density may be expressed as supported temporal sampling functions result in probability distributions that fall even more slowly [16], because c = α/p in Eq. (1), where 0 < α < 1 and p depends on the quantum operator. Here α, to be defined in Sec. II, is a parameter that describes the rate of switching.
For example, for the energy density we have p = 3. Note that compactly supported functions are more representative of measurements that occur in a finite period of time when compared to Lorentzian functions, which have an infinitely long tail.
An independent confirmation of the behavior of the tail of the probability distribution in Eq. (1) for the cases of Lorentzian and compactly supported temporal sampling functions is done in Ref. [17]. By performing a direct diagonalization of the time-averaged square of the time derivative of a massless scalar field in Minkowski spacetime, the authors numerically evaluate P (x) for large vacuum fluctuations in a spherical cavity with finite degrees of freedom. The fitted values of the parameters {c 0 , a, b, c} that govern the asymptotic behavior of P (x) are reported in Ref. [17] and are in good agreement with those obtained in Refs. [15,16] based on the moments approach.
The moments approach applied to quantum stress tensor operators averaged over finite time intervals is further developed in Ref. [18] to include averaging over finite regions of space. If the spatial sampling scale is small compared to the temporal scale, the asymptotic behavior of P (x) for a spacetime-averaged quadratic operator is expected to first decay as the worldline limit discussed earlier, with c = α/p in Eq. (1), before smoothly transitioning to a form with c = α instead. Although the inclusion of spatial averaging increases the decay rate of the probability distribution compared to time averaging alone, the distribution still falls more slowly than exponentially for large x. Under the right conditions, large vacuum fluctuations could then produce several physical effects that overshadow those of thermal fluctuations, opening up the possibility of experimental or observational confirmation. In the following paragraphs, we briefly mention some of the most striking effects.
Fluctuating gravity waves produced by quantum stress tensor fluctuations of a conformal field in inflationary models are studied in Ref. [19]. These gravity waves are potentially observable in the cosmic microwave background radiation and from gravity wave detectors, providing a probe of trans-Planckian physics.
Large vacuum radiation pressure fluctuations on particles with electric charge or nonzero polarizability may push the particles over potential barriers as shown in Ref. [20]. Depending on the details of the averaging over the finite time interval, the penetration rate via this mechanism may even surpass the known quantum tunneling rate.
The vacuum decay of a metastable state of a self-interacting scalar field is analyzed in Ref. [21]. Large quantum fluctuations of the time derivative of a scalar field averaged over a finite spacetime region lead to a decay rate comparable with the standard rate from the instanton approximation [22]. However, for operators that are quadratic in the time derivative of a scalar field, the probability distribution falls slower than an exponential function, in which case the decay rate is governed by these quadratic field fluctuations rather than quantum tunneling and linear field fluctuations.
Large fluctuations around the zero point density of a fluid as an analog model for quantum stress tensor fluctuations is studied in Ref. [23]. These density fluctuations may potentially be detectable in low-temperature light scattering experiments [24] by observing fluctuations in the number of scattered photons.
In this paper, we extend the diagonalization approach developed in Refs. [17,25] to treat probability distributions of quantum stress tensor operators averaged over a finite spacetime region. The paper is organized as follows. In Sec. II, we review the asymptotic behavior of the probability distributions of spacetime-averaged stress tensor operators, based on the high moments approach of Ref. [18]. In Sec. III, we consider the alternative diagonalization method, developed in Ref. [17], that allows us to numerically construct the probability distributions. We apply this numerical method to the square of the time derivative of a massless scalar field, and in Sec. IV we discuss the approximations and analyze the results.
In Sec. V, we review the key takeaways and remaining loose ends.
Units in which the reduced Planck constant and the speed of light are equal to unity, = c = 1, are used throughout the paper.

ATORS
The moments of a quantum operator can be used to infer the properties of the underlying probability distribution. However, the moments of a quadratic field operator, which composes the stress tensor operator, are not well-defined at a single spacetime point, complicating efforts to do so. One workaround is to investigate the moments of a quadratic field operator that has been averaged in time alone or space and time. 1 One is further led to consider averaging over a finite duration and volume, which are more representative of physical measurements in an experiment. Consider a normal-ordered, quadratic field operator, which can be expanded in terms of creation and annihilation operators in the form The moments of the quadratic operator T (t, x) generically diverge. In order to obtain finite moments, we need to average T (t, x) in time alone, to find or in space and time, to find where f (t) and g(x) are the temporal and spatial sampling functions, respectively, and µ n is the nth moment. We assume the integrals of f (t) in time and g(x) in space are normalized 1 Note that quadratic operators averaged in space alone still have diverging moments in four spacetime dimensions, as discussed in footnote 2 of Ref. [16].
The finite moments of the time-averaged or spacetime-averaged operators can be related to the moments of a probability density function, where x denotes the eigenvalues of the operator in question, and the lower integral bound −x 0 comes from quantum inequalities. Quantum inequalities are constraints on the expectation values of averaged stress tensor operators. Because these expectation values can be arbitrarily negative when evaluated at a single spacetime point [26], macroscopic violations of physical laws become possible, a problem that can be resolved by arguing for quantum inequality constraints [27]. Note that −x 0 is the lowest eigenvalue of the averaged operator, and hence is both the minimum expectation value and the lower bound on the probability distribution.
The case for a time-averaged quadratic operator is discussed in detail in Ref. [16], and here we proceed to summarize the main results. Let us consider the normal-ordered quadratic operator T (t, x) = τ 4 (:φ 2 (t, x) :) .
Here ϕ(t, x) is the quantized massless scalar field, soφ 2 (t, x) has dimensions of (length) −4 in units where = c = 1. We introduce the extra factor of τ 4 to make the operator dimensionless, where τ is the characteristic temporal sampling scale.
In rectangular coordinates, ϕ(t, x) has the usual solution where V is the quantization volume and ω = k. Averaging T (t, x) in time, we find We assume f (t) is a compactly supported, real, symmetric sampling function with a Fourier that asymptotically approachesf Here C f and β > 0 are constants and 0 < α < 1. Note that τ , the characteristic sampling scale, may be defined by Eq. (15) as the decay scale of the Fourier transform. For the functions which will be used in this paper, this is of the same order as the characteristic duration of f (t), but this need not be true in general. Assuming a fixed spatial location x = 0, we find that F k,k' (x = 0) and G k,k' (x = 0) in Eq. (4) are given by By calculating the moments in Eq. (5), we find that, for large n, there is one dominant term given by the expression Qualitatively, we may see this by noting that F kk (0) falls off more slowly than G kk (0) due to the arguments off (ω) in Eqs. (16) and (17). Since a † k a k F kk (0) annihilates the vacuum state, we need G kk (0) and G * kk (0) placed at the ends. M n contains the maximum possible number of factors of F kk (0), which leads to its dominance over other terms in µ n . Taking the continuum limit, done in detail in Sec. IV of Ref. [16], we find that M n when n 1 is of the order The Hamburger and Stieltjes moment theorems [28], applied to distributions on whole lines and half-lines, respectively, guarantee unique probability distributions provided the moments do not grow too quickly with n. For distributions that are bounded below, as in the case of stress tensor operators subject to quantum inequality constraints, the Stieltjes moment theorem may be more relevant. The moments in Eq. (19) grow faster than the criteria of either moment theorem, so we cannot guarantee these moments specify a unique distribution. However, it can be shown that distributions with moments given in Eq. (19) can be different from the previously referenced asymptotic form, Eq. (1), by merely some oscillatory function, leaving the salient features unaffected. Calculating the moments of the probability distribution in Eq. (1) using Eq. (10), we find Comparing Eqs. (19) and (20), we identify Numerical simulations performed in Ref. [17] based on the direct diagonalization of the averaged operator T (0) find good agreement with these predictions.
The moments approach may be readily extended for spacetime-averaged operators, which is discussed in Ref. [18]. In this case we consider the spacetime-averaged analog of Eq. (13), where the choices of ϕ(t, x) and f (t) are identical to the time-averaged case, and g(x) is a compactly-supported, real, spherically symmetric sampling function with a Fourier trans- that asymptotically approacheŝ Here C g is a constant, is the characteristic sampling length scale, and 0 < λ < 1 with λ ≤ α. Note that the factor of k λ−2 arises in a specific function constructed in Ref. [18], but need not appear more generally. The F kk and G kk matrix elements in Eq. (7) are The dominant contribution to the moments in Eq. (8) is assumed to be the form given in Eq. (18) for the same reasons as the time-averaged case. As in the time-averaged case, the high moments can be approximated and compared to the moments of the proposed P (x) in Eq. (1). Interestingly, the Stieltjes moment theorem holds when α > 1/2, suggesting a unique P (x) can be determined in those cases. The analytical calculation in Ref. [18], though similar to that of the time-averaged case in Ref. [16], is unable to precisely predict most of the parameters in Eq. (1) due to poor understanding of the regime where the approximations hold. The unambiguously predicted parameters in Eq. (1) are c and, with a caveat, a; the other parameters are not well-known. The spacetime-averaged distribution is expected to eventually decay as Here B is a constant predicted for a class of sampling functions in Sec. V D in Ref. [18] but otherwise not known in general. Note that we are ignoring possible overall factors in powers of x before the exponential in Eq. (27), which are well predicted for the time-averaged case in Eq. (1). If < τ , the effect of spatial averaging is not pronounced for the lower moments, suggesting that the worldline behavior found in the time-averaged case approximately holds in some regime. For a given nth moment of the probability distribution P (x) in Eq. (10), one can estimate the location x that contributes most to the integral, which depends on n.
Using an estimate of the highest moment for which the time averaging is dominant, we find that the worldline regime holds for x x * , where The asymptotic behavior of P (x) for a spacetime-averaged quadratic operator is then expected to first decay as the worldline limit in Eqs. (1) and (21) until around x ∼ x * , where the distribution transitions to the form in Eq. (27).

III. DIAGONALIZATION OF QUADRATIC BOSONIC OPERATORS
Here we proceed to generalize the diagonalization procedure in Ref. [17] for time-averaged quantum stress tensor operators to include averaging over finite spatial volumes.

A. General procedure
Our goal is to numerically evaluate x and P (x) for an arbitrary spacetime-averaged quadratic operator T in the Minkowski vacuum state |0 a . To do so, we need to solve for the eigenvalues of T and corresponding probabilities of measuring those eigenvalues in the vacuum state, which amounts to a diagonalization problem.
Recall that a generic spacetime-averaged operator can be expanded in the form given in Eq. (7). Here we will assume that F and G are real and symmetric matrices, so that In general, the vacuum state will not be an eigenstate of T . We perform a Bogoliubov transformation [29] to convert T into a diagonal form, with creation and annihilation oper- acting on a different set of particle number states labeled by the subscript b, Such a transformation is done assuming {a k , a † k } can be written as a linear combination of {b k , b † k }, which obey their own sets of commutation relations (see Sec. III in Ref. [17] for further details). Requiring the diagonal where 1 is the identity operator, the conditions outlined above give expressions for λ k and C shift , constants that depend on F kk and G kk . In analogy with ladder operators in quantum mechanics, we can express the a vacuum state |0 a in terms of the b particle number states |{n k } b through clever use of the creation and annihilation operators. Recall that the a vacuum is the physical state in which we wish to study the fluctuations. It can be shown Here T refers to the matrix transpose and b denotes a column matrix composed of b k 's.
The matrix M is derived by noting that a k , which can be written as a linear combination annihilates the a vacuum state |0 a , which itself is a linear combination of the b number states |{n k } b . The constant N emerges from the usual normalization a 0|0 a = 1.
Doing the calculation in full, both M and N can be derived from F kk and G kk .
Thus, for any vector in the eigenbasis |{n k } b , the eigenvalue and corresponding measurement probability in the original a vacuum state |0 a is given by B. Specific case: square of the time derivative of a massless scalar field We are interested in the specific case discussed earlier, with the spacetime-averaged operator T given in Eq. (22) and the sampling functions behaving as discussed in Eqs. (15) and (24). The Klein-Gordon equation for a massless scalar field is Although Sec. II is done in rectangular coordinates, here we work in spherical coordinates to take advantage of the spherical symmetry of the spatial sampling function g(r). The solution for the positive frequency mode function is given by where ω = k, j l (r) are the spherical Bessel functions, P m l (x) are the associated Legendre functions, ξ ωlm is some phase factor, and A ωlm is some constant to be determined. A convenient choice for the phase factor is Applying vanishing boundary conditions at the surface of a sphere of radius R, we find where z nl is the nth zero of the spherical Bessel function j l (kr). Since the frequencies ω depend on n and l, it is more convenient to label the solutions with {n, l, m} instead of {ω, l, m}. Requiring that the commutation relations hold in second quantization, we find For the Condon-Shortley phase convention, there is an extra factor of (−1) m . Expanding ϕ(t, r) in terms of creation and annihilation operators, Differentiating in time, squaring the result, and ordering normally gives where H.c. refers to the Hermitian conjugate. The spacetime average of :φ 2 (t, r) : can be done by recalling the definition of the Fourier transform, Eq. (14), and making use of the orthonormality conditions of the spherical harmonics. We find, identically to Eq. (29), 2F nlm,n l m a † nlm a n l m + G nlm,n l m (a † nlm a † n l m + a nlm a n l m ) , where Here we have assumed that g(r) has a compact support of [0, r 0 ]. With the F and G matrices known, the procedure in Sec. III can be performed to construct P (x), which is done in Sec.
IV for the case n = 1 ∼ 600, l = m = 0. In the limit of no spatial averaging, we expect to recover the purely time-averaged result. Indeed, if we let we find The integral vanishes for l = 0 because j l (0) = 0 in those cases. For l = 0 we have Y 00 (θ, φ) = Y * 00 (θ, φ), so the result holds for all four terms in Eq. (41). We then get F n00,n 00 = (nn ) 3/2 τ 4 π 2 2(−1) n+n R 4f and G n00,n 00 = −(nn ) 3/2 τ 4 π 2 2(−1) n+n R 4f where it is understood that the matrix elements with l = 0 are zero. Equations (47) and (48) are identical to the results in Ref. [17], except for the extra factor of 1/(−1) n+n , which arises due to a different choice of the phase factor ξ nlm .
As discussed in Sec. II, the asymptotic behavior of P (x) is primarily determined by the Fourier transforms of the sampling functions, so it can be useful to rewrite Eqs. (43) and (44) in terms ofĝ(k) instead of g(r). This can be done by calculating the Fourier transform using the plane wave expansion, giving where P l (x) are the Legendre polynomials. Equation (49) can be substituted into Eqs. (43) and (44)  A compactly supported time sampling function with a Fourier transform that asymptotically approaches the limit in Eq. (15) can be constructed following Sec. IIB of Ref. [16]. Let where it can be shown that specification of the parameters {α, δ, τ } thus generates a particular time sampling function and corresponding Fourier transform.
Although we could perform the above procedure for the full set of {ω} in the computation, in practice there are a number of complications that motivate an approximation. For a computation with many frequency modes, calculatingf (ω) point-by-point is time consuming and susceptible to numerical error at large ω. We also need to differentiate and integratê f (ω), a task made easier with an analytic form. We choose to approximatef (ω) in the following way:f The spline interpolation is performed on a sample dataset with ω ∈ [0, ω c ]. When ω > ω c , we directly evaluate the theoretically expected form, Eq. (15). Here ω c is chosen to be a point where the numerically computedf (ω) approaches the theoretically expected limit but before any severe numerical error sets in.

Construction of g(r) orĝ(k)
A compactly supported, spherically symmetric spatial sampling function g 1 (r) that has a Fourier transform asymptotically approachingĝ 1 (k) ∼ e −(k ) λ can be found using the method in Sec. IID of Ref. [16]. Although this method was originally used to construct a onedimensional temporal sampling function in Ref. [16], the argument holds for constructions of spherically symmetric spatial sampling functions, which take only a single argument r.
The asymptotic behavior ofĝ 1 (k) strongly depends on the properties of g 1 (r) near the end points, where g 1 (r) switches on and off. Suppose we want g 1 (r) compactly supported in r ∈ [0, r 0 ]. Because we assume a spherically symmetric g 1 (r), the relevant behavior is the switch-on and switch-off as r → r + 0 and r → r − 0 , respectively. Note that the function does not switch on or off at r = 0, which is in the interior of the sampling region.
Crudely, to getĝ 1 (k) ∼ e −(k ) λ , direct application of the method in Sec. IID of Ref. [16] requires g 1 (r) to switch on as e −r λ/(λ−1) as r → 0 + . In our case, since there is no switch-on at r = 0, we instead want e −(r 0 −r) λ/(λ−1) as r → r − 0 , which is just a reflection and translation to convert the switch-on at r = 0 to a switch-off at r = r 0 . For the case λ = 0.5, one option is then Here A = (3e)/{2πr 3 0 [8 + 13e Ei(−1)]} is a normalization factor such that d 3 r g(r) = 1. While this construction gives us a Fourier transform that asymptotically approachesĝ 1 (k) ∼ e −(k ) λ , there are limitations to this method. We do not have fine control over the Fourier transform itself, which is the most important function in the context of the high moments, as discussed in Sec. II. In particular, the characteristic spatial sampling scale is unknown, complicating efforts to calculate the transition location x * . We are also not able to guarantee that the asymptotic limit goes exactly as Eq. (24), which is the assumed behavior in Ref. [18].
For the functions used in this paper, the characteristic decay length ofĝ is of the same order as the spatial sampling length.
For this reason, there are benefits to constructing the Fourier transform directly. Here we consider the construction discussed in Ref. [18]. Suppose we have a one-dimensional, compactly supported time sampling function h(t) with a Fourier transform that asymptotically approachesĥ (ω) ∼ C h e −η|ωτ | λ , ωτ 1.
We now define a different spatial sampling function for which the Fourier transform iŝ Hereĥ (ω) ≡ d dωĥ (ω), and is now an input parameter we control. Equation (58) can be found by explicitly taking the derivative ofĥ(ω) in Eq. (57), using the asymptotic limit in Eq. (55). The construction given in Eq. (57) has the same asymptotic behavior as Eq. (24) once we identify respectively. The most straightforward choice forĥ(ω) is then to chooseĥ(ω) =f (ω). In this case, the compact support of g 2 (r) is [0, 2δ /τ ].

Particle sectors and n,l,m
In principle, the diagonalization of a quadratic field operator calls for multiple infinite sums. We may readily see this from Eq. 2F nn a † n a n + G nn (a † n a † n + a n a n ) , where, recalling that the nth zero of j 0 (kr) is z n0 =nπ, and The equation relating the integral with g(r) to the integral withĝ(k), Eq. (49), can now be written in the simpler form The other upper bound to consider emerges from Eqs. (32) and (33), which ask for b particle number states |{n k } b to calculate the eigenvalues and probabilities. As there are infinitely many b number states, we need to choose which states to include in the computation. Let us first consider the a vacuum state |0 a , which can be written in terms of |{n k } b using Eq. (31) and a Taylor expansion of the exponential: The second sum in Eq. (65) is bounded above at i, j = 600 due to our choice of n = 1 . . . 600, l = m = 0. The first sum over ρ controls the number of b and b † 's in the expansion of the a vacuum state. Note that the b and b † 's always come in pairs, so the expansion of |0 a in the basis of b number states is such that only the even particle sectors of |{n k } b contribute: ρ = 0 is the b zero-particle sector, ρ = 1 is the b two-particle sector, and so on.
In the expansion for |0 a , Eq. (65), the only term that matters is the one that is proportional to |2 p b , because all other terms vanish when taking the inner product.
Now observe that the 2ρth particle sector contains a factor of (M ij ) ρ in Eq. (65). The probabilities in the 2ρth sector will then go as |M ij | 2ρ . Since |M ij | < 1, for fixed i and j the higher particle sectors contribute smaller and smaller probabilities. In light of the averaging process discussed later in section IV A 4, this suggests that the higher particle sectors can be ignored without affecting the probability distribution significantly. In our setup we will calculate the eigenvalues and probabilities of only the 2-particle sector, corresponding to ρ = 1. However, note that in a typical computation |M ij | for different i and j can span many orders of magnitude, so it is not true that the higher particle sectors always contribute negligible probabilities. However, given the necessity for limiting the particle sectors in the computation, the most straightforward choice is to focus on the two-particle sector, which contributes significantly to the probability distribution across the board. Here it may be worth recalling the results in the worldline case. In Ref. [17], Table I, some results for cumulative probability distributions are given. In all the cases studied there, the fourparticle sector gives a contribution of the order of 4% or less of that of the two-particle sector.
The two-particle sector also offers the advantage of being simpler to manage, as there are only two possible configurations of the momentum states. If both particles are in the same momentum state, the eigenvalues and probabilities are given by If the two particles are in different momentum states, we instead get As we go to higher particle sectors, the number of configurations rises quickly, raising the additional question of which configurations to include in the computation, a problem we avoid with the two-particle sector.

Averaging P (x)
The probability distribution P (x) can be constructed by calculating the eigenvalues from Eq. (32) and probabilities from Eq. (33). Note that these equations do not provide any a priori reason to expect that P (x) is a smooth distribution. The arguments in Sec. II that lead to a smooth distribution rely on the moments of a quadratic operator, which do not necessarily encode the finer details of the distribution. The analytical treatment assumes that the moments of a quadratic operator can be related to the moments of a smooth probability distribution function, which is a sensible conjecture but not proven.
Indeed, for a generic computation following Sec. III, the constructed probability distribution is highly degenerate: for eigenvalues that are close together, the probabilities of measuring those eigenvalues can vary significantly. The simplest way to see this is an extension of the argument in Sec. IV A 3. Higher b number states are more likely to have smaller probabilities of being measured, but the outcomes of these measurements are not guaranteed to be much different from lower particle number states. As an example, let us compare a measurement with two particles in the same momentum state with a measurement with four particles in the same momentum state. The former case is given in Eqs. (68) and (69). The analogous expressions for four particles are One can imagine finding states p and p for which the outcomes are similar, but the probabilities are vastly different, In fact, to avoid the described situation across the infinitely many particle sectors would require much coincidence on the part of nature. Perhaps such fine tuning in fact occurs in nature, but for a finite mode computation we need to deal with a degenerate P (x).
In experimental and observational settings, we may not need to worry about these degeneracies. No realistic physical measurement probes a single eigenvalue x of an operator, so the probability distribution P (x) will always be integrated over some finite region of x. We claim that the theoretical predictions in Sec. II are describing this physically observable probability distribution instead, one that has already been coarse-grained by the measurement process. This coarse-grained distribution is the one we expect to be smooth and asymptotically approach the theoretical predictions.
Under this view, repeated measurements probing some range ∆x = x f − x i will find, on or, for the discrete probability distributions we are dealing with, Here the sum over j is understood to be over all eigenvalues x j in the measurement domain, where the analogous expression for the discrete case is Numerically, our prescription is to average the raw data using bins of width ∆x, where the values ofx andP (x) in each bin are given by Eqs. (77) and (79), respectively. 2 Doing so allows us to coarse-grain the raw data in a manner consistent with what we might expect from physical measurements.
The physical meaning of the bin sizes ∆x is not so clear. One could argue that the sizes of these bins correspond to limiting factors in an experiment, such as the resolution of a detector or the uncertainty in the momentum of a photon probe, but such claims are purely speculative and would need to address the distinction between the averaging via binning discussed here and the spacetime averaging of the operators discussed in Section II. As the averaged data is fairly independent of the bin size ∆x, our numerical simulations do not rely on particular choices of ∆x. In this sense we can also view the binning as a mathematical tool used to better analyze the numerical results and leave the trickier question of physical meaning for future investigation.

B. Results
Because the full asymptotic form of the probability distribution is not well-predicted, as shown in Eq. (27) compared with Eq. (1), our analysis needs to accommodate the undetermined parameters. We would like to verify the unambiguous theoretical predictions: asymptotically, the averaged probability distribution will first decay as before transitioning to a form that decays as where the location of this transition is expected to occur at x ∼ x * , given in Eq. (28). In Eqs. (80) and (81), it is understood that Eq. (1) gives the full asymptotic form, of which we are concerned with the parameter c that governs the exponential fractional decay rate. Here we have written the equations using the averaged quantitiesP (x) andx to emphasize that we expect the theoretical predictions to hold for the coarse-grained probability distributions.
From Eq. (1) we have assumingx is sufficiently large. Note that Eq. (1) itself already assumesx 1, but here we requirex to be even greater so that the third term in Eq. (82) is dominant. We thus find A plot of ln[− lnP (x)] against lnx then gives a slope of c, which is α/3 in the worldline limit, Eq. (80), and α in the spacetime-averaged limit, Eq. (81).
As discussed in Sec. IV A 2, it is more advantageous to work withĝ 2 (k) instead of g 1 (r) to allow for finer control over the Fourier transform and better compatibility with Ref. [18]. However, in principle the numerical diagonalization can be performed using either function, as shown in Fig. 1, where the averaged probability distributions are remarkably similar despite g 1 (r) = g 2 (r). In particular, we see that two sampling functions with similar asymptotic behavior in Fourier space do indeed produce similar probability distributions, lending credence to the claim that the Fourier transforms govern the tail region of the probability distributions. We work with datasets that perform the numerical diagonalization using the Fourier transforms of the sampling functions,f (ω) andĝ 2 (k), and we consider the cases α = λ = 0.5; α = 0.7, λ = 0.5; and α = λ = 0.7. The parameters used in the construction off (ω) and g 2 (k) are shown in Table I, where we explicitly work in τ = 1 units. In principle, we would like to take the limit where the boundary of the sphere is infinite, i.e., R → ∞. Because we are limited to a finite number of modes, we are forced to consider a finite boundary. A compromise is to take 2R or 4πR 3 /3 to be greater than the sampling times and volumes, respectively, so that the presence of the boundary would not be observable in a measurement.
We satisfy this requirement for our sampling volumes, but numerical instabilities do not allow us to choose small enough sampling times without sacrificing the high frequency modes. The latter problem is also noted in Ref. [17], and as in that publication we acknowledge that the total sampling duration, 4δ, is greater than 2R, the distance to travel from the origin to the boundary and back, or from one end of the boundary to the opposite end. The simulations may thus be an imperfect approximation of Minkowski space, but we expect the setup to be approximate enough to compare with the theoretical calculations.
For each of the three cases, the averaging bin size ∆x is chosen to reduce scatter in the spacetime-averaged regime without smearing out the worldline behavior completely. For Table I this reason, datasets for which the worldline behavior ends earlier require smaller ∆x. As we do not expect the measurement outcomes x, and thus the averaged outcomesx, to be correlated, a least-squares linear fit is sufficient for our purposes. However, in order to account for the different number of data points in each bin, we perform weighted leastsquares linear fits to find the slopes. As a shorthand, let us write and We want a linear fit to Eq. (84), which can now be written as where the expected values for γ i are given in Eq. (84): γ 1 = d, γ 2 = ln a. We estimate γ i by finding the parametersγ i such that the weighted squared residuals are minimized, wherez(y;γ 1 ,γ 2 ) is a linear fit function. Here w i is the weight associated with the ith squared residual; for non-weighted least-squares fits the conventional choice is to take w i = 1 for all i. In our numerical computation, w i is taken to be the number of points in the ith bin divided by the total number of points.
The fit results are compiled in Table II and Fig. 2. The numerical results for the exponential decay rates in the spacetime-averaged limits match the predictions exceedingly well.
In contrast, the decay rates in the worldline limits are lower than expected in all cases. One explanation is the largex approximation in Eq. (83) may not hold for our worldline data in the region lnx = 4 ∼ 9. Moreover, as shown in Fig. 3, datasets with greater have less data in the worldline region, limiting the reliability of our fits. Nonetheless, the worldline limit has been extensively investigated in Ref. [17], where numerical simulations verified the the two limits, ln x * , from Eq. (28). The parameters used for each dataset are given in Table II. full asymptotic form in Eqs. (1) and (21), so the poorer worldline results here are likely due to limitations in our analysis or data sets.
As shown Fig. 3, the predicted transition locations from worldline to spacetime-averaged behavior are underestimates, suggesting the worldline behavior is more robust than predicted. We see that the transition behavior qualitatively holds, as datasets with smaller do transition at greaterx, though the calculations leading to Eq. (28) may be too rough to accurately pinpoint the locations of these transitions.
We can consider a more detailed analysis of the transition behavior by noting that which can be numerically probed via linear fits. The actual transition locations can be numerically estimated through a number of methods, one of which we discuss here. Theoretically, x * is predicted to be the point at which the worldline behavior ends. At this point, by definition the data will deviate from the linear fits in the worldline region. Following the notation in Eqs. (85)-(88), we can estimate these points by taking ln x * to be the smallest provided (y i , z i ) is not a data point used in the linear fit, which has been assumed to be part of the worldline regime. Here is a constant that sets the threshold for how much the data need to deviate from the worldline behavior, so > 1. Some estimates of ln x * are given in Fig. 4, where we observe that larger values of put the transition further from the worldline region. Note that Eq. (90) is equivalent to the fractional residual form Linear fits to Eq. (89) with fixed {α, λ} and varying do not conclusively confirm the theoretical predictions for the choices of attempted, which is not surprising given the discrepancies in Fig. 3. Early results suggest where B 1 can vary by a few orders of magnitude and B 2 is roughly in the ballpark of 3.
However, for each case of {α, λ} we only have four or five different datasets, limiting the reliability of the fit results. Furthermore, as shown in Fig. 3, for large values of there is very little data in the worldline regime to begin with, heavily skewing the estimates of ln x * . The omission of the higher l modes could exacerbate this problem. One trial dataset using l = 0, 1 shows more data at smaller values of x, which could shift the location of the transition upon coarse-graining. We leave for the future more careful analysis of the transition behavior.  Table I.

V. OUTLOOK AND DISCUSSION
Large fluctuations of stress-tensor-like operators have a number of potentially observable and interesting effects, including fluctuating gravity waves [19], increased barrier penetration probabilities for charged particles [20], alternative processes for false vacuum decay [21], and greater variance of scattered photons in a low-temperature light scattering experiments [23].
The extent to which these effects have physically observable manifestations depends on the likelihood of these large fluctuations, which can be investigated through the underlying probability distributions. The asymptotic behavior of these distributions can be deduced from the high moments of the quadratic operators in question, although the operators need to be averaged in time alone or space and time for the moments to be finite. Purely timeaveraged operators have been discussed in two dimensions with Gaussian time sampling functions [13], in four dimensions with Lorentzian time sampling functions [15], and more recently in four dimensions with compactly supported time sampling functions [16]. The latter two scenarios were numerically verified in Ref. [17], providing confirmation of the high moments method. However, a physical experiment takes place not only in a finite duration but also in a finite volume, motivating work on stress tensor operators averaged by compactly supported space and time sampling functions [18]. In four dimensions, the probability distributions of spacetime-averaged operators asymptotically approach the worldline limit described by a purely time-averaged operator, P (x) ∼ e −x −α/3 , before transitioning to a form that decays faster, P (x) ∼ e −x α , where x is a dimensionless quantity proportional to the eigenvalues of the operator and the value of α determines the switch-on/off behavior of f (t).
In this paper, we adapt the method developed in Ref. [17] for the case of spacetimeaveraged operators. The Minkowski vacuum state is generally not an eigenstate of an arbitrary normal-ordered quadratic operator averaged in space and time, so repeated measurements of the operator lead to different outcomes with different probabilities of occurence. To construct the associated probability distribution, a Bogoliubov transformation is performed to find the eigenvalues and eigenkets. The probabilities of measuring these eigenvalues in the vacuum state is then given by the squared inner products of the eigenkets with the vacuum state. Choosing a suitable set of eigenkets, we can numerically construct the probability distributions by finding the eigenvalues and corresponding probabilities.
Numerically constructing the distribution for the spacetime average of :φ 2 (t, r) :, where ϕ(t, r) is the massless scalar field, we find that similar outcomes can have wildly different probabilities of being measured. As physical measurements are not precise enough to probe single eigenvalues, we argue that binning the data is a plausible resolution that produces smoother, more well-behaved data sets. Alternative, rigorously developed coarse-graining methods may perhaps already exist in other fields, but in any case our procedure should capture the key qualities we expect from these averaged distributions. Whether the bin sizes carry any physical meaning is speculated but better left for future investigation.
Our results show clear worldline and spacetime-averaged behavior with obvious transitions between the two limits, allowing analysis of the asymptotic behavior of the probability distributions. Fitting to the asymptotic regions of the averaged probability distributions for the cases α = λ = 0.5; α = 0.5, λ = 0.7; and α = λ = 0.5, we find that the decay rates in spacetime-averaged limits are consistent with prediction, whereas those of worldline limits are slightly lower than expected. The latter inconsistency may result from limitations of the datasets and the approximations used to analyze the transition to and behavior in the spacetime-averaged limit. When the full asymptotic form in the spacetime-averaged limit is predicted theoretically, we expect that a more careful analysis of the numerical data will find good consistency in both the worldline and spacetime-averaged limits.
In contrast, the predicted locations x * for the transitions between these two limits are somewhat inconsistent with data, though the qualitative behavior holds. We do not have enough datasets to numerically analyze the transition behavior in detail, a problem compounded by the lack of data in the worldline region in some computations. Preliminary analysis suggests the predicted power law behavior for x * , Eq. (28), may hold, but we lack sufficient data to conclusively show this. Though the effects of higher l and m modes are not well understood, we may expect nontrivial contributions because the spherical Bessel functions are nonzero as one moves away from the origin. Such effects could shift the transition location, a speculation that appears plausible from some early trial runs. Further exploration of the consequences of the l and m modes, perhaps in conjunction with more n modes, may be worth pursuing in the future.
In addition to the transition behavior, another work in progress is the generalization to rectangular coordinates. Although in this paper we work in spherical coordinates to take advantage of the spherical symmetry present in setup, spherical coordinate systems can be unwieldy in many contexts. Computations in rectangular coordinates could allow for more straightforward applications in a variety of scenarios and perhaps even facilitate more sophisticated simulations, but the selection of modes can be tricky and more work remains to be done. In spherical coordinates, a sufficient choice for mode selections is to set l = m = 0, but the analogous choice in rectangular coordinates is not so clear.

VI. ACKNOWLEDGMENTS
We thank Chris Fewster for valuable discussions. This work was supported by the

Academy of Finland Grant 318319 and by the National Science Foundation under Grant
No. PHY-1912545.