Vacuum Quantum Stress Tensor Fluctuations: A Diagonalization Approach

Large vacuum fluctuations of a quantum stress tensor operator can be described by the asymptotic behavior of the probability distribution of the time or spacetime averaged operator. Here we focus on the case of stress tensor operators averaged with a sampling function in time. The Minkowski vacuum state is not an eigenstate of the time-averaged operator, but can be expanded in terms of its eigenstates. We calculate the probability distribution and the cumulative probability distribution for obtaining a given value in a measurement of the time-averaged operator taken in the vacuum state. In these calculations, we use the normal ordered square of the time derivative of a massless scalar field in Minkowski spacetime as an example of a stress tensor operator. We analyze the rate of decrease of the tail of the probability distribution for different temporal sampling functions, such as compactly supported functions and the Lorentzian function. We find that the tails decrease relatively slowly, as exponentials of fractional powers, in agreement with previous work using the moments of the distribution. Our results lead additional support to the conclusion that large vacuum stress tensor fluctuations are more probable than large thermal fluctuations, and may have observable effects.


I. INTRODUCTION
The definition and the use of the expectation value of a quantum stress tensor operator have been a topic of intense study in recent decades. The semiclassical theory for gravity uses the renormalized expectation value of the quantum matter stress tensor to give an approximate description of the effects of quantum matter fields on the gravitational field.
As in the semiclassical theory of electromagnetic radiation, it is expected that this theory is a reasonable approximation to a more complete quantum theory of gravity coupled to matter fields. It is known that a renormalized stress energy operator for quantum fields in curved spacetime is associated with quantum corrections to Einstein's equations, via higher order derivative terms [1]. These corrections lead to physical effects, such as small scale factor oscillations around an expanding background universe and quantum particle creation [2].
Moreover, this theory has been successful about giving a plausible description of the back reaction to black hole evaporation through Hawking radiation [3]. However, the semiclassical theory does not consider the quantum fluctuations of the stress tensor around its expectation value and their possible effects. Several authors have studied a variety of physical effects associated with quantum stress tensor fluctuations [4]. These effects include, for example, potentially observable gravity waves from quantum stress tensor fluctuations in inflationary models [5], effects of vacuum electric field fluctuations on light propagation in nonlinear materials [6,7], and barrier penetration of charged or polarizable particles through large vacuum radiation pressure fluctuations [8,9].
In general, the physical effects of large fluctuations of a quantum stress tensor operator can be studied through the analysis of the probability distribution for the time or spacetime averaged operator. This probability distribution can be inferred (at least qualitatively) from the moments of the averaged operator, and the exact distribution was found in a twodimensional model in Ref. [10]. The moments method was used in Ref. [11] to infer the probability distribution for several normal-ordered quadratic operators in four dimensional Minkowski spacetime with Lorentzian time averaging. These included the square of the electric field and the energy densities of a massless scalar field and of the electromagnetic field.
This idea was extended in Ref. [12]  can dominate over thermal fluctuations at large energies. However, the moments of a quantum stress tensor operator grow very rapidly, to the extent that they might not uniquely determine the probability distribution, so it is desirable to seek alternative methods.
In the present paper, we develop such an independent test of the moments approach for the study the probability distribution of time-averaged quantum stress tensor operators.
The main idea is to diagonalize the time-averaged operator through a change of basis and calculate the cumulative probability distribution function of their quantum fluctuations in the vacuum state. We are interested in checking the behavior predicted by the high moments approach, and in determining which modes and particle numbers give the dominant contribution to the large fluctuations. Unlike the moments approach, which primarily gives information about the asymptotic behavior of the probability distribution for large vacuum stress tensor fluctuations, the diagonalization approach in principle gives a unique probability distribution for a broad range of fluctuations x. We take the normal ordered square of the time derivative of a massless scalar field in Minkowski spacetime as our stress tensor operator, and find the tail of the probability distribution for different temporal sampling functions, specifically a class of compactly supported functions and the Lorentzian function.
The tails decrease relatively slowly, as exponentials of fractional powers, in agreement with previous results using the moments of the distribution.
The paper is organized as follows: In Sec. II, we review the main results of Ref. [12] on the high moments approach to the analysis of the probability distribution for quantum stress tensor operators. In Sec. III, we develop an independent approach to the study of probability distributions based on the diagonalization of the operator. In Sec. IV, we show the numerical results obtained for different time sampling functions. In Sec. V, we summarize and discuss the main results of the paper.

II. MOMENT-BASED APPROACH TO THE PROBABILITY DISTRIBUTION
Here we review the main results of Ref. [12]. Working in 4-dimensional Minkowski spacetime, let T (t, x) be a operator which is a quadratic function of a free field operator and define its time average with a real-valued sampling function f (t) by We will consider measurements of the time average T rather than T . The sampling function has a characteristic width τ and should decay quickly as |t| τ . One example is a Lorentzian function, used in Ref. [11], whose mathematical expression and Fourier transform are given by where the Fourier transform of f L (t) and its normalization are given bŷ However, if the measurement of the operator occurs in a finite interval of time, the sampling function is better described by a smooth and compactly supported function. This kind of sampling function is strictly zero outside a finite region, avoiding the long temporal tails of functions like the Lorentzian. It therefore gives a better description of a measurement which begins and ends at finite times. We will be interested in compactly supported nonnegative functions whose Fourier transform has the following asymptotic form when ωτ 1: where α, γ, and β are constants. Here α ∈ (0, 1) is a decay parameter which defines the rate of decrease off (ω) (values α ≥ 1 are incompatible with f having compact support). It is worth emphasising that τ does not directly measure the support of f , but rather indicates the shortest characteristic timescale associated with f ; in our examples, this will characterise the switch-on and switch-off regions.
For any given f (compactly supported or not) define the n-th moment of the normalordered time-averaged quadratic operator T , Eq. (2.1), as where |0 is the Minkowski vacuum vector of the theory. As we will now see, the form of the Fourier transformf defines the rate of growth of the moments µ n and, as a result, the probability for large fluctuations.
In the first instance, we work in a box of finite volume and express T in a mode sum of creation and annihilation bosonic operators as whereÃ ij andB ij are components of symmetric matricesÃ andB, which have the functional where ω i are the mode frequencies. Precise forms ofÃ andB will be given when we come to specific examples in Section IV. The moment µ n can be expressed as an n-th degree polynomial in these components. As n increases, the number of terms in the expression for the n-th moment grows rapidly. Fortunately, only one term gives the dominant contribution for n 1: First,B j 1 j 2 andB * jnj 1 have to begin and end, respectively, the expression for M n becausẽ B * ij a † i a † j andB ij a i a j in Eq. (2.6) are the only terms which do not annihilate the vacuum from the left and right, respectively. Second, all the remaining coefficients in M n areÃ ij 's, which fall slower thanB ij 's when ω i becomes large. This arises because theÃ ij involve a difference in frequencies, as opposed to the sum in theB ij . Provided thatf ≥ 0, all the terms contributing to the n-th moment are nonnegative, so M n is actually a lower bound on µ n , which will gives us a lower bound on the probability distribution for large vacuum fluctuations.
To be more specific, now consider the time average of :φ 2 :, where φ is a massless scalar field in four-dimensional Minkowski spacetime. Then, passing to a continuous mode sum, the dominant term takes the form (2.10) Iff has the asymptotic form (2.4), then the dominant term has the asymptotic form, in units in which τ = 1, [12]). The most important part of this expression is the gamma function factor, which leads a rapid rate of growth of the high moments, M n ∝ (3n/α)!. Thus, the parameter α is crucial in determining the rate of growth of the moments when n 1.
The goal is to use the asymptotic form for the moments, Eq. between this feature of the stress tensor probability distribution and quantum inequality bounds, which is explained in detail in Refs. [10,11]. We define the tail distribution (also called the complementary cumulative distribution function), P > (x), as the probability of finding any value y ≥ x in a measurement x P (y)dy (2.12) and of course P is normalized so that P > (x) = 1 for x ≤ −x 0 . The n-th moment of T can be written in terms of P as and this can be compared with the the asymptotic form of the dominant contribution M n , Eq. (2.11), to infer information about P (x) and P > (x). In this way, we are led to consider the asymptotic forms for large vacuum fluctuations, x 1, where c 0 , a, b, and c are constants to be determined, and for which the corresponding moments obey when n becomes large. The similarity between this expression and the asymptotic form for M n , Eq. (2.11), is evident, and leads to the identifications (2.16) However, the situation is a little bit more subtle, because it is not guaranteed that a set of moments growing as fast as (3n/α)! (for α < 1) determines a unique probability distribution [14]. Fortunately, the difference between two probability distributions with the same moments is just an oscillatory function, which does not add any interesting feature to the general form of P (x) for our purposes. Therefore the parameters in Eq. (2.16) should provide a good approximation to the asymptotic behavior of P (x) and P > (x). Rigorous arguments to this effect are given in Sec. VI of [11].
The argument just given applies to the case of a compactly supported function with asymptotics given by Eq. (2.4). For the case of a non-compactly supported sampling function such as a Lorentzian, Eq. (2.2), a slightly different argument is needed to compute the asymptotic form of the dominant contribution M n , as is explained in detail in Ref. [11].
However, the analysis of high moments still leads to an asymptotic form for P (x) given by Eq. (2.14) with c = 1/3. This is consistent with the α → 1 limit of the relation c = α/3 derived for compactly supported functions, in which limit the asymptotic form (2.4) agrees with that of the Lorentzian (2.2), with γ = β = 1.
In general, we see that the decay parameter α in the asymptotic form of the sampling function's Fourier transform determines the rate of decay in P (x) for large x, and hence the probability of large vacuum fluctuations. The smaller α is, the more slowly the tail decreases and the greater the probability of large fluctuations becomes. For compactly supported functions, the value of α is related to the rate of switch-on and switch-off of f (t).
[See Eqs. (51) and (52) in Ref. [12].] and occupation numbers to the probability distribution, in addition to providing a uniquely defined probability distribution.

A. Bogoliubov diagonalization
We express a general quadratic operator H as a mode sum involving bosonic creation and annihilation operators for N modes as and 1 1 is the identity operator. Here the coefficients of Eq. (3.1) correspond to elements of N -square matrices {D r } 4 r=1 which form the so-called dynamical matrix Here we follow an approach developed by Colpa [15] for the diagonalization of D. This approach was previously applied to stress tensor operators by Dawson [16], who was primarily concerned with quantum inequality bounds on expectation values. The diagonalization of the quadratic operator H implies a homogeneous linear transformation (Bogoliubov transformation [17]) to go from the original set of bosonic operators, in which H takes a diagonal form. For our purposes, we consider the case D 1 = D 4 = F and D 2 = D 3 = G with F and G real and symmetric matrices. Under these conditions, we may normal order the operator H in Eq. (3.1) to obtain Now we apply a Bogoliubov transformation where A and B are real N ×N matrices, and the new set of bosonic operators satisfy the usual Note that the commutation relations for the a and a † operators and the Bogoliubov transformation, Eq. (3.6), impose conditions upon A and B matrices of the form where I and 0 are the identity and null N × N matrices, respectively. A consequence of Substituting Eq. (3.6) into Eq. (3.5), we obtain Now we impose a diagonalization condition in Eq. (3.8), where Λ = diag(λ 1 , . . . , λ N ). Using the canonical commutation relations for the b i , we obtain It is clear that :H: is diagonal in the orthonormal basis formed by vectors where n = (n 1 , . . . , n N ) with each n i a nonnegative occupation number, so that b † i b i |n b = n i |n b and |0 b is annihilated by all the b i . The eigenvalues are easily read off from where the i-index runs from 1 to N , and a sum on repeated indices is understood. The operator :H: is bounded from below provided that λ 1 , . . . , λ N are all nonnegative, in which case C shift is the lowest eigenvalue. This gives a quantum inequality bound ψ|:H:|ψ ≥ C shift (3.14) for all physical normalized states ψ. Note that C shift is both the lowest eigenvalue of the time-averaged stress tensor operator, and the lower bound on its probability distribution, Let us return to the problem of achieving the diagonalization in practice. Noting that Eq. (3.7) can be written in matrix notation as we use the diagonalization condition, Eq. (3.9), to obtain which is equivalent to a set of 2N -equations to be solved for A, B, and Λ, given F and G: and as we are interested in the case where Λ is positive definite, it follows that a solution is only possible if both F + G and F − G are also positive definite. In this case, the equations can be solved as follows. First, because F −G is positive, we may use the Cholesky decomposition [18] to find a real and invertible matrix K such that K † K = F − G. The matrix K(F + G)K † is real, symmetric and positive definite and therefore can be brought where all the diagonal entries are strictly positive and U is a real orthogonal matrix. We then define It may be verified (see Appendix A) that the solution to (3.17) and (3.18) is given by Λ together with Express where C n are coefficients to be determined. Apply the a-annihilation operator from the left and use the Bogoliubov transformation for the single mode case, Eq. (3.6), according to As the |n b form an orthonormal basis, we deduce to obtain It can be easily proved by induction that the general term in Eq. (3.28) has the form (3.29) where N ≡ C 0 . We apply the normalization condition to obtain N as As a result, the probability P 2n of finding the original a-vacuum state in a specific b-state, |2n b , and the corresponding eigenvalue of T , Eq. (3.10), are given by where n = 0, 1, 2, 3 . . . . From these equations, we see that the lowest possible outcome in a measurement of T is just C shift , the a-vacuum state is only connected with 2n-particle sectors of the b-state, and the probability of finding the a-vacuum state in a specific b-state is concentrated in the lower particle number sectors. Indeed, the asymptotic expression for P 2n decreases rapidly with n according to P 2n ∼ |N | 2 M 2n √ πn , for large n and |M| < 1. (3.34) These three features of the single mode case hold for the general case that we now proceed to develop in the next subsection.

C. Probabilities for particle sectors and outcomes for the general case
Express the a-vacuum state as a linear combination of ψ n , where each ψ n belongs to the n-particle subspace for b-states, as follows where we have used b † k b k ψ n = nψ n in the last line. The expression inside the bracket in Eq. (3.38) consists of n-particle terms with n ≥ 2, thus we can only have a solution with ψ 1 = 0. That means that and ψ 1 = ψ 3 = · · · = ψ 2n+1 = 0. We can rewrite Eq. (3.39) by relabeling n by 2n and expressing ψ 2n in terms of ψ 0 as where N is a normalization constant to be determined. Now, we diagonalize M such that M = S T ΞS with S a real and orthogonal matrix and Ξ = diag(µ i ). Set c i = S ij b j and They satisfy the bosonic commutation relations, because Here we have used the commutation relation of b-operators and the orthogonality of S. Now note that we can rewrite the exponent in Eq. (3.42) using The normalization constant N is calculated using the a-vacuum normalization and the definition for c † i 's. For a single mode we have (3.53) From this expression, we can determine, for example, the probability of finding the system in the b-vacuum state, P {0} , or in some configuration in the two-particle sector such as P {2 i } or P {1 i 1 j } . These probabilities and the corresponding outcomes associated with a measurement of T , Eq. (3.10), are specifically given by Here it is understood that i, j run from 1 to N . Now, we can re-express the normalization constant, N , to obtain information about the total probability for each particle sector. We take Eq.  Then, the contributions to the total probability of the b-vacuum and the two-particle sector, for instance, are |N | 2 and (1/2)|N | 2 Tr(M 2 ), respectively. Each 2n-particle sector contributes with terms having 2n-factors of M.

IV. MASSLESS SCALAR FIELD IN MINKOWSKI SPACETIME
A. The square of the time derivative of the field We consider a minimally coupled massless scalar field, φ(x), in a four-dimensional Minkowski spacetime in spherical coordinates (t, r, θ, ϕ), with the origin of the spherical polar coordinates placed at the fixed spatial point at which :φ 2 : will be evaluated. The equation of motion is given by the usual wave equation Solutions of this equation take the form [19] where g ωl (r) = ω 2 R j l (ωr) , Here z nl is the n-th zero of the spherical Bessel function, j l .
We expand the quantized field in terms of creation and annihilation operators, a ωlm and a † ωlm , as where a sum on ω is abbreviated notation for the sum on n = 1, 2, . . .
where a ω ≡ a ω00 , the sums run over the range given in Eq. (4.10), and H.c. means hermitian conjugate. Convergence here should be understood in a distributional sense, so that when we now let T = :φ 2 : in Eq. (2.1), we find wheref is the Fourier transform of the sampling function f (t), Eq. (2.3).
, (4.14) β = 2 cos πα 2 . We define dimensionless variables x 1 = T (τ 2 ) 2 and x 2 = T (4πτ 2 ) 2 for the compactly supported functions and the Lorentzian function, respectively. (The difference in the numerical factors is to facilitate comparison with the results of Refs. [11] and [12], which used slightly different conventions.) Using the expression for ω, Eq. (4.10), these variables become τ 4 0 2π 2 (rs) 3/2 a † r a sf (|r − s|τ 0 ) − a r a sf ((r + s)τ 0 ) + H.c. , (4.16) and 8τ 4 0 (rs) 3/2 a † r a s e −|r−s|τ 0 − a r a s e −(r+s)τ 0 + H.c. , (4.17) where we have defined The F and G matrices are all that we need to calculate, for a given number of modes, the probability distribution and the cumulative probability distribution associated with a measurement of x 1 or x 2 .

B. Numerical results for the cumulative probability distribution function and tail for large fluctuations
Here we explain the general features of the numerical calculation that we carry out to calculate the probability, P (x), and cumulative probability distribution function, P > (x), for the two cases mentioned above. Here x denotes either x 1 or x 2 , defined in Eqs. (4.16) and (4.17). For a given number of modes, we calculate all possible outcomes in a measurement of x up to and including the 6-particle sector, except for the following outcomes which have been omitted: and (4.15). For a fixed characteristic timescale, τ , the radius of the sphere is inversely proportional to dimensionless variable τ 0 . For a given number of modes, if the size of the sphere is too large, there will not be enough data in the tail (x 1) of the P > (x)-curve to perform a reliable fit. By contrast, if the size of the sphere is too small, the P > (x)-curve will not be smooth, showing a step-like behavior. For the compactly supported functions, we also have to determine values for δ, which defines the support of the sampling function f (t), i.e., the duration of the sampling. We choose these values to be slightly larger than the first maximum of the corresponding ϕ(t), and the results are given in Table I, working in units where τ = 1. Then the sphere radius R = τ 0 /π gives values 1.14, 0.64 and 0.64 for α = 0.5, 0.6, 0.7, respectively, for the values of τ 0 considered. Note that for α = 0.5 we have R > 2δ, which means that the total sampling time is less than the time taken for light to travel to the boundary and back. Accordingly, the numerics ought to give a good approximation to sampling in Minkowski space; this is an instance of local covariance, which has a number of applications to quantum inequalities [20]. By contrast, in the other two remaining cases we have R < 2δ, so the sampling process can be sensitive to the presence of the bounding sphere. The reduced values of δ used for α = 0.6, 0.7 were required to obtain numerical stability.
We build P > (x)-curves for compactly supported functions whose Fourier transform is given by Eq. x max is the maximum value obtained in a measurement of x for a given number of modes  Fig. 2, for the case of compactly supported functions with different values of α and for the Lorentzian function. Units in which τ = 1 have been adopted. Here values of P > (x) for different particle sectors are calculated adding all probabilities for all possible outcomes for the given sector as is indicated in Table III. Since x max is the maximum value obtained in a measurement of x for a given number of modes and size of the sphere, the expression [1 − P > (x max )] gives us the loss of probability. and size of the sphere. All analyzed cases show a small loss of probability of the order of 10 −5 or less. This small loss of probability indicates that the outcomes which have been included provide a reasonable approximation for P > (x).
Our calculated values of the lower bound C shift can be compared with results from other approaches. In the case of the Lorentzian function, our calculated value C shift = −0.0593338, is of the order of the predicted value from the analysis using high moments, x 0 = −0.0236 [11] and well within the (non-optimal) theoretical bound C shift ≥ −27/128 = −0.211 given by the method of [21]. For a general compactly supported test function f , the theoretical bound which can be obtained by setting p = √ ω in Eq. (3.11) of [21]. For the case α = 0.5, the integral on the right-hand side of (4.23) can be evaluated numerically and yields the bound C shift ≥ −0.3592. Our calculated value C shift = −0.0781613 is therefore consistent with the theoretical bound and indicates that the latter bound is weaker than the sharpest  Table I. possible bound by a factor of approximately 4.6. This result is broadly in line with Dawson's computations [16], where a ratio of about 3 was found. Note that Dawson used a toroidal spatial geometry rather than a ball and a squared Lorentzian sampling function of infinite duration, so one would not expect an exact match with our results.
Since we want to test the predicted behavior of the cumulative probability distribution for large fluctuations in vacuum, we focus on the tail of each P > (x)-curve and propose a trial function inspired by Eq. (2.14). Specifically,  The fitting procedure is based on the least-squares method to find the specific set of values of parameters which minimize the error variance. We name this specific set as θ * = (p * 1 , a * , b * , c * , c * 0 ). The estimation of the error variance, s 2 , is given by   Figure 3 shows the P > (x)-curves with their respective best fits to the trial function, Eq. (4.24). In the case of the Lorentzian function, the P > (x)curve and its respective fit are indistinguishable on the scale shown. Figure 3 shows that the diagonalization procedure is able to reproduce smooth tails for all the cases considered, which are well fitted by the trial function given by an incomplete gamma function, Eq. (4.24). The variance of the fits are small in comparison to the variation of the P > (x)-curves. For example, for the case of the compactly supported function with α = 0.7, we have s 2 ∼ O(10 −22 ) but the change of the P > (x)-curve over the range plotted in Fig. 3 is the order of 10 −9 . All values of parameters obtained through the best-fitting procedure agree reasonably well with the predicted ones from the high moments approach [11,12], In complete agreement with previous results based on the high moments analysis [11,12], our results confirm the fact that averaging over a finite time interval compactly supported functions results in a probability distribution which falls more slowly than for the case of the Lorentzian function, and both fall more slowly than exponentially.

V. SUMMARY AND DISCUSSION
Large vacuum fluctuations of quantum stress tensor operators can have a variety of physical effects such as production of gravity waves in inflationary models [5], fluctuations of the light propagation speed in nonlinear materials [6,7], and enhancing barrier penetration of charged or polarizable particles [8,9]. These quantum fluctuations can be studied through the analysis of the probability distribution for the time or spacetime averaged operator in Minkowski spacetime. The asymptotic behavior of the probability distribution can be inferred by studying the moments of the normal ordered operator. The study of several normal-ordered quadratic operators time averaged with a Lorentzian function [11] or compactly supported functions [12] predict an asymptotic form of the probability distri-  The value of α is related to the rate of switch-on and switch-off of compactly supported functions.
In the present paper, we have developed a method which is independent of the moments approach for the study of the probability distribution for quantum vacuum fluctuations of a time averaged quantum stress tensor operator, T , in Eq. (2.1). Since the vacuum state is not usually an eigenstate of T , we diagonalize this operator through a change of basis.
Expressing the vacuum state in terms of the new basis in which T is diagonal, we are able to calculate the probability distribution, P (x) and the cumulative probability distribution function, P > (x) for obtaining a specific result in a measurement of T . Specifically, we work with the time averaged quadratic operator T = ∞ −∞ :φ 2 (t, 0):f (t)dt, where φ is a massless minimally coupled scalar field and f (t) is the sampling function . We use a dimensionless variable x ∝ T τ 4 , where τ is a characteristic timescale of the sampling function. Numerical results for both Lorentzian and compactly supported functions show that the probability distribution of vacuum quantum fluctuations is bounded below at x = −x 0 < 0, and that the tail of the probability distribution varies as an incomplete gamma function in agreement with the previous studies [11,12]. We apply a best-fit procedure through a least-squares method to the tail of the P > (x)-curves in order to determine values for parameters in Eq. (4.24). The results for p 1 , a, b, and c parameters show good agreement with the predictions of the high moments approach. (See Table II.) The diagonalization procedure is able to reproduce with great accuracy the rate of decrease of the tail of the cumulative probability distribution.
We reproduce the relation c = α/3 for α = (0.5, 0.6, 0.7, 1), where α = 1 corresponds to the case of the Lorentzian function, with percentage errors less than 6% compared to the theoretical values predicted by the high moments approach [11,12]. Our results confirm that averaging over a finite time interval, with compactly supported functions, results in a probability distribution which falls more slowly than for the case of Lorentzian averaging, and both fall more slowly than exponentially.
Recall that we have quantized the scalar field in a sphere with finite radius R, so the probability distribution which we calculate could differ from that of empty Minkowski spacetime.
As was noted in Sect. IV B, there should be no difference for the α = 1/2 case, as the duration of the sampling is less than the light travel time to the boundary and back. In the other cases, there could in principle be an effect of the boundary. However, this is likely only to alter the lower frequency modes, which are not expected to give a large contribution to the tail of the distribution.
The diagonalization method is free of the ambiguity potentially present in the high moments approach, and leads to a unique result for the probability distribution. It also has the potential to determine the entire distribution, including its lower bound, which is also the optimum quantum inequality bound on expectation values of the averaged operators.
In addition, it can provide information about the particle content of the eigenstates of the averaged stress tensor which are associated with the large fluctuations.

VI. ACKNOWLEDGMENTS
We would like to thank Tom Roman for valuable discussions in the early stage of this project. CJF thanks Simon Eveson for a useful discussion concerning the numerical evaluation of (4.23). This work was supported in part by the National Science Foundation under Grant PHY-1607118.